Finite-Scale Emergence of 2+1D Supersymmetry at First-Order Quantum Phase Transition

TL;DR

This paper demonstrates that 2+1D supersymmetry can emerge at a first-order quantum phase transition in an Ising magnet, with potential experimental realization using topological superconductor surfaces, revealing finite-scale supersymmetry phenomena.

Contribution

It shows for the first time that supersymmetry can emerge at a first-order quantum phase transition in 2+1D systems, expanding the understanding of SUSY in condensed matter.

Findings

Supersymmetry emerges at the first-order transition of the Ising magnet.

Emergence of SUSY is linked to a topological phase transition of Majorana fields.

The scale of SUSY emergence can be large with small initial masses.

Abstract

Supersymmetry, a symmetry between fermions and bosons, provides a promising extension of the standard model but is still lack of experimental evidence. Recently, the interest in supersymmetry arises in the condensed matter community owing to its potential emergence at the continuous quantum phase transition. In this work, we demonstrate that 2+1D supersymmetry, relating massive Majorana and Ising fields, might emerge at the first-order quantum phase transition of the Ising magnetization by tuning a single parameter. Although the emergence of the SUSY is only allowed in a finite range of scales due to the existence of relevant masses, the scale range can be large when the masses before scaling are small. We show that the emergence of supersymmetry is accompanied by a topological phase transition for the Majorana field, where its non-zero mass changes the sign but keeps the magnitude. An experimental realization of this scenario is proposed using the surface state of a 3+1D time-reversal invariant topological superconductor with surface magnetic doping.

Figures

Tables

Table 1. Table 1: Classification of partial derivatives and Pauli matrices according to the angular momentum L z subscript 𝐿 𝑧 L_{z} .

$L_z$	Expressions
0	${\partial }_{\tau },{\partial }_{\tau }^2,{\sigma }_0,{\sigma }_y,{\nabla }^2$
-1	${\partial }_x-i{\partial }_y,{\sigma }_z-i{\sigma }_x$
1	${\partial }_x+i{\partial }_y,{\sigma }_z+i{\sigma }_x$
-2	${\left({\partial }_x-i{\partial }_y\right)}^2$
2	${\left({\partial }_x+i{\partial }_y\right)}^2$

Equations

$$S=\int d^{d}x[\frac{1}{2}{\gamma }^T({\partial }_{\tau }-i{v}_f\mathbf{\alpha }\cdot \nabla +m{\sigma }_y)\gamma +\frac{1}{2}g\varphi {\gamma }^T{\sigma }_y\gamma$$

$$+\frac{1}{2}\varphi (-{\partial }_{\tau }^2-v_b^2{\nabla }^2+r)\varphi +a\varphi +\frac{1}{3!}b{\varphi }^3+\frac{1}{4!}u{\varphi }^4],$$

$$v_f=v_b=1,b=3mg,u=3g^2,$$

$$a=\frac{m}{g}(r-m^2),r<\frac{3m^2}{2}.$$

$$S=\int d^{d}x[\frac{1}{2}{\gamma }^T({\partial }_{\tau }-i{v}_f\mathbf{\alpha }\cdot \nabla +m^{\prime }{\sigma }_y)\gamma +\frac{1}{2}g\bar{\varphi }{\gamma }^T{\sigma }_y\gamma$$

$$+\frac{1}{2}\bar{\varphi }(-{\partial }_{\tau }^2-v_b^2{\nabla }^2+r^{\prime })\bar{\varphi }+\frac{1}{3!}{b}^{\prime }{\bar{\varphi }}^3+\frac{1}{4!}u{\bar{\varphi }}^4],$$

$$\frac{d{v}_b}{dl}=\frac{{\tilde{g}}^2\left(v_f^2-v_b^2\right)}{32{\pi }^2{v}_b{v}_f^3},\frac{d{v}_f}{dl}=\frac{{\tilde{g}}^2(v_b-v_f)}{6{\pi }^2{v}_b(v_b+v_f{)}^2},$$

$$\frac{d\tilde{g}}{dl}=\tilde{g}\left(\frac{ϵ}{2}-\frac{7{\tilde{g}}^2}{32{\pi }^2}\right)$$

$$\frac{d\tilde{u}}{dl}=\frac{12{\tilde{g}}^4-2{\tilde{g}}^2\tilde{u}-3{\tilde{u}}^2}{16{\pi }^2}+\tilde{u}ϵ,$$

$$\frac{dm}{dl}=m-\frac{3{\tilde{g}}^{2}m}{16{\pi }^2}$$

$$\frac{d\tilde{b}}{dl}=\tilde{b}(1+\frac{ϵ}{2})+\frac{\tilde{b}\left(-3{\tilde{g}}^2-6\tilde{u}\right)+24{\tilde{g}}^{3}m}{32{\pi }^2}.$$

$$\frac{d}{dl}\left(\frac{b}{mg}\right)=-\frac{{\tilde{g}}^2}{4{\pi }^2}(\frac{b}{mg}-3).$$

$$\frac{dr}{dl}=2r-\frac{{\tilde{b}}^2+{\tilde{g}}^2\left(r-6m^2\right)+r\tilde{u}}{16{\pi }^2}$$

$$\frac{d\tilde{a}}{dl}=(3-\frac{ϵ}{2})\tilde{a}-\frac{\tilde{a}{\tilde{g}}^2+2\tilde{b}r-4\tilde{g}{m}^3}{32{\pi }^2},$$

$$\frac{d}{dl}\left(\frac{r}{m^2}\right)=\frac{{\tilde{g}}^2}{8{\pi }^2}\left[\left(\frac{r}{m^2}\right)-\frac{3}{2}\right]$$

$$\frac{d}{dl}\left(\frac{ag}{m^3}-\frac{r}{m^2}+1\right)=\frac{5{\tilde{g}}^2}{16{\pi }^2}\left(\frac{ag}{m^3}-\frac{r}{m^2}+1\right),$$

$$S^{\ast }=\int d^{d}x[\frac{1}{2}{\gamma }^T(i{\partial }_{\mu }{\alpha }^{\mu }+m^{\ast }{\sigma }_y)\gamma +\frac{1}{2}{g}^{\ast }\varphi {\gamma }^T{\sigma }_y\gamma$$

$$+\frac{1}{2}\varphi (-{\partial }^2+r^{\ast })\varphi +a^{\ast }\varphi +\frac{1}{2}{g}^{\ast }{m}^{\ast }{\varphi }^3+\frac{1}{8}(g^{\ast }{)}^2{\varphi }^4],$$

$$S_{MGNY}=\int d^{d}x[{\bar{\Psi }}_j({\partial }_{\tau }{\bar{\gamma }}^0+v_f\bar{\mathbf{\gamma }}\cdot \nabla +m){\Psi }_j+g\varphi {\bar{\Psi }}_j{\Psi }_j$$

$$+\frac{1}{2}\varphi (-{\partial }_{\tau }^2-v_b^2{\nabla }^2+r)\varphi +a\varphi +\frac{1}{3!}b{\varphi }^3+\frac{u}{4!}{\varphi }^4],$$

$$\frac{d{v}_b}{dl}=\frac{N{\tilde{g}}^2\left(v_f^2-v_b^2\right)}{32{\pi }^2{v}_b{v}_f^3},\frac{d{v}_f}{dl}=\frac{{\tilde{g}}^2(v_b-v_f)}{6{\pi }^2{v}_b(v_b+v_f{)}^2}.$$

$$\frac{d\tilde{g}}{dl}=\tilde{g}\left(\frac{ϵ}{2}-\frac{(6+N){\tilde{g}}^2}{32{\pi }^2}\right)$$

$$\frac{d\tilde{u}}{dl}=\frac{12N{\tilde{g}}^4-2N{\tilde{g}}^2\tilde{u}-3{\tilde{u}}^2}{16{\pi }^2}+\tilde{u}ϵ$$

$$\frac{dm}{dl}=m-\frac{3{\tilde{g}}^{2}m}{16{\pi }^2}$$

$$\frac{d\tilde{b}}{dl}=\tilde{b}(1+\frac{ϵ}{2})+\frac{\tilde{b}\left(-3N{\tilde{g}}^2-6\tilde{u}\right)+24N{\tilde{g}}^{3}m}{32{\pi }^2}$$

$$\frac{dr}{dl}=2r-\frac{{\tilde{b}}^2+N{\tilde{g}}^2\left(r-6m^2\right)+r\tilde{u}}{16{\pi }^2}$$

$$\frac{d\tilde{a}}{dl}=(3-\frac{ϵ}{2})\tilde{a}-\frac{N\tilde{a}{\tilde{g}}^2+2\tilde{b}r-4N\tilde{g}{m}^3}{32{\pi }^2}.$$

$$\frac{d{v}_b}{dl}=(2{\pi }^2{K}_d)\frac{{\tilde{g}}^2\left(v_f^2-v_b^2\right)}{32{\pi }^2{v}_b{v}_f^3}$$

$$\frac{d{v}_f}{dl}=(2{\pi }^2{K}_d)\frac{{\tilde{g}}^2(v_b-v_f)}{6{\pi }^2{v}_b(v_b+v_f{)}^2},$$

$$\frac{d\tilde{g}}{dl}=\tilde{g}\left(\frac{ϵ}{2}-(2{\pi }^2{K}_d)\frac{7{\tilde{g}}^2}{32{\pi }^2}\right)$$

$$\frac{d\tilde{u}}{dl}=(2{\pi }^2{K}_d)\frac{12{\tilde{g}}^4-2{\tilde{g}}^2\tilde{u}-3{\tilde{u}}^2}{16{\pi }^2}+\tilde{u}ϵ$$

$$\frac{dm}{dl}=m-(2{\pi }^2{K}_d)\frac{3{\tilde{g}}^{2}m}{16{\pi }^2}$$

$$\frac{d\tilde{b}}{dl}=(1+\frac{ϵ}{2})\tilde{b}+(2{\pi }^2{K}_d)\frac{\tilde{b}\left(-3{\tilde{g}}^2-6\tilde{u}\right)+24{\tilde{g}}^{3}m}{32{\pi }^2}$$

$$\frac{dr}{dl}=2r-(2{\pi }^2{K}_d)\frac{{\tilde{b}}^2+{\tilde{g}}^2\left(r-6m^2\right)+r\tilde{u}}{16{\pi }^2}+{\Delta }_r$$

$$\frac{d\tilde{a}}{dl}=(3-\frac{ϵ}{2})\tilde{a}-(2{\pi }^2{K}_d)\frac{\tilde{a}{\tilde{g}}^2+2\tilde{b}r-4\tilde{g}{m}^3}{32{\pi }^2}+{\Delta }_a,$$

$$\frac{d}{dl}\left(\frac{ag}{m^3}-\frac{r}{m^2}+1\right)=\frac{5}{8}{\tilde{g}}^2{K}_d\left(\frac{ag}{m^3}-\frac{r}{m^2}+1\right),$$

$$\text{I}:r^{\ast }>(m^{\ast }{)}^2,⟨\varphi ⟩=0$$

$$\text{II}:r^{\ast }<(m^{\ast }{)}^2,m^{\ast }>0,⟨\varphi ⟩=\frac{-3m^{\ast }-\sqrt{9(m^{\ast }{)}^2-8r^{\ast }}}{2g^{\ast }}$$

$$\text{III}:r^{\ast }<(m^{\ast }{)}^2,m^{\ast }<0,⟨\varphi ⟩=\frac{-3m^{\ast }+\sqrt{9(m^{\ast }{)}^2-8r^{\ast }}}{2g^{\ast }}.$$

$$S_0=\int \frac{d\tau d{k}^3}{(2\pi {)}^3}[{\bar{\psi }}_{\tau ,\mathbf{k}}({\partial }_{\tau }+h(\mathbf{k})){\psi }_{\tau ,\mathbf{k}}+\frac{1}{2}{\bar{\psi }}_{\tau ,\mathbf{k}}{\Delta }_{\mathbf{k}}({\bar{\psi }}_{\tau ,-\mathbf{k}}{)}^T$$

$$+\frac{1}{2}{\psi }_{\tau ,-\mathbf{k}}^T{\Delta }_{\mathbf{k}}^{†}{\psi }_{\tau ,\mathbf{k}}]$$

$$S_E=\int d^{d}x[\frac{1}{2}{\gamma }^T({\partial }_{\tau }-i{v}_f\mathbf{\alpha }\cdot \nabla )\gamma +\frac{1}{2}g\varphi {\gamma }^T{\sigma }_y\gamma +$$

Full text

I Introduction

Originally proposed as a means to evade the Coleman-Mandula theorem and unify internal and space-time symmetries, supersymmetry (SUSY) has been an active area of research due to its potential in solving the hierarchy and cosmological constant problems and other puzzles in high-energy physics Gervais and Sakita (1971); Wess and Zumino (1974); Zumino (1975); Dimopoulos and Georgi (1981); Wess and Bagger (1992). Bosons and fermions related by SUSY transformations, referred to as superpartners, have the same mass Wess and Bagger (1992); SUSY breaking lifts this degeneracy proportionally to the breaking scale. Despite extensive searches at high energies, conclusive experimental evidence for SUSY is yet to be found.

The last 30 years have witnessed active SUSY-related research in solid-state systems, including the tricritical Ising model Friedan et al. (1984), the boundary of topological insulators and superconductors Hasan and Kane (2010); Qi and Zhang (2011); Ponte and Lee (2014); Grover et al. (2014); Zerf et al. (2016); Witczak-Krempa and Maciejko (2016); Li et al. (2017, 2018); Jian et al. (2017), the bulk of semimetals Lee (2007); Jian et al. (2015), high-$T_{c}$ superconductors Balents et al. (1998), the Josephson-junction array Foda (1988) and various other model systems Thomas (2005); Fendley et al. (2003); Huijse et al. (2008); Bauer et al. (2013); Rahmani et al. (2015); Huijse et al. (2015), as well as in the cold atom system Yu and Yang (2010). It has been argued that, even though the corresponding microscopic models do not exhibit it, SUSY emerges macroscopically at continuous quantum phase transitions of such solid-state systems. The gapless nature of the continuous phase transition implies that the resulting superpartners are massless.

In this work, we present an example of emergent SUSY at the first-order quantum phase transition (FOQPT) of a solid-state system. The FOQPT can be achieved by tuning only one parameter and the corresponding superpartners are massive. Although the emergent SUSY is only valid in a finite range of scales owing to the gapped nature of the FOQPT, the scale range can be large if the inital mass (mass before scaling) is small. Specifically, we consider a 2+1D Majorana field coupled to an Ising field, and perform the one-loop renormalization group (RG) analysis in three different schemes. Within the finite range of scales allowed by the gapped theory, all three schemes show that the FOQPTs can be reached by tuning one parameter, and have emergent SUSY with gapped Majorana and Ising fields serving as the massive superpartners. Interestingly, the emergent SUSY is always accompanied by a topological phase transition of the Majorana field even though its mass does not vanish. Finally, we propose an experimental realization of the emergent SUSY based on the time-reversal (TR) invariant topological superconductor (TSC).

II SUSY in a Massive Majorana-Ising System

We consider a 2+1D action that describes a Majorana fermion $\gamma$ interacting with an Ising field $\phi$ and discuss its SUSY. The 2+1D action is

[TABLE]

where $x=(\tau,\boldsymbol{x})$ with the imaginary time $\tau$, and $\boldsymbol{\alpha}=(\sigma_{z},\sigma_{x})$ with the Pauli matrices $\sigma_{i}$. Without loss of generality, we choose $v_{f}>0$ by rotating the index of the Majorana field, $g\geq 0$ by flipping the sign of $\phi$, and $u>0$ to make the bosonic potential bounded below. Eq. (II) has rotational invariance along $z$ and, for a uniform vacuum expectation value (VEV) of $\phi$, it is the most general rotationally invariant action to $\phi\gamma^{T}\sigma_{y}\gamma$ and $\phi^{4}$ order. (See more details in Appendix. A) The action (II) is not invariant under the TR transformations, $\gamma_{\tau,\boldsymbol{x}}^{T}\rightarrow\gamma_{\tau,\boldsymbol{x}}^{T}(i\sigma_{y})$ and $\phi_{\tau,\boldsymbol{x}}\rightarrow-\phi_{\tau,\boldsymbol{x}}$, unless $m,b,a=0$. The TR-invariant case was analyzed in Ref. [Grover et al., 2014], where it was shown that SUSY with massless superpartners emerges by tuning the parameter $r$ to the continuous phase transition point $r=0$.

Eq. (II) has SUSY when

[TABLE]

To see this, let us perform the vacuum shift $\phi=\bar{\phi}+\phi_{0}$ to Eq. (II) with $\phi_{0}$ satisfying $a+r\phi_{0}+\frac{b\phi_{0}^{2}}{2}+\frac{u\phi_{0}^{3}}{3!}=0$, resulting in

[TABLE]

where $m^{\prime}=m+g\phi_{0}$, $r^{\prime}=r+b\phi_{0}+\frac{1}{2}u\phi_{0}^{2}$, and $b^{\prime}=b+u\phi_{0}$. When the SUSY condition (II) holds, $\phi_{0}$ has three different choices: $\phi_{0}=(-m\pm\sqrt{3m^{2}-2r})/g,-m/g$, while only the first two are vacua of the Eq. (II). Around either of the two vacua, the SUSY condition (II) further leads to $r^{\prime}=m^{\prime 2}$ and $b^{\prime}=3m^{\prime}g$ in Eq. (3), making Eq. (3) the “real” version of the 2+1D Wess-Zumino SUSY model. Indeed, it is invaraint under infinitesimal SUSY transformation: $\delta_{\xi}\bar{\phi}=\xi^{T}\sigma_{y}\gamma$ and $\delta_{\xi}\gamma=\sigma_{y}\alpha^{\mu}\xi(i\partial_{\mu}\bar{\phi})+\xi(-m^{\prime}\bar{\phi}-g\bar{\phi}^{2}/2)$, where $\xi=(\xi_{1},\xi_{2})^{T}$ are the constant Grassmann-valued parameters of the SUSY transformation Wess and Bagger (1992), $\mu=0,1,2$, $\partial_{\mu}=(i\partial_{\tau},\boldsymbol{\nabla})$ and $\alpha^{\mu}=(-\mathds{1},-\boldsymbol{\alpha})$. (See Appendix. B for more details.) Therefore, Eq. (II) with Eq. (II) has SUSY with $\gamma$ and $\phi$ serves as superpartners.

Now the question becomes how to realize the SUSY in Eq. (II) since Eq. (II) is not typically satisfied. Naively, one may think 4 parameters need to be fine-tuned to realize Eq. (II), i.e. tuning $(\frac{v_{f}}{v_{b}},\frac{b}{mg},\frac{u}{g^{2}},\frac{ag}{m^{3}}-\frac{r}{m^{2}})$ to $(1,3,3,-1)$, as the velocities $v_{f},v_{b}$ can always be chosen as 1 by rescaling the spacial coordinate once they are equal and $r<3m^{2}/2$ is a parameter region instead of a critical condition. In the following, we show through one-loop RG analysis that the first three of the above parameters can naturally flow to the SUSY-required values as the scale increases, resulting in the emergence of SUSY achievable by finely tuning only one parameter $(\frac{ag}{m^{3}}-\frac{r}{m^{2}})$.

III One-loop RG Equations and Emergent SUSY at Finite Scales

The one-loop RG analysis is performed in three schemes in $d=4-\epsilon$ dimensionss Lee (2007); Fei et al. (2016): (i) dimensional regularization (DR) for Eq. (II), (ii) DR from the so-called “massive Gross-Neveu-Yukawa (GNY) model” Fei et al. (2016), and (iii) Wilson RG scheme with spacial momentum cutoff Shankar (1994) for Eq. (II). All three schemes show that the SUSY might emerge as the scale increases by tuning only one parameter. Furthermore, we discuss the finite-scale nature of the emergent SUSY.

III.1 DR and Emergent SUSY at finite Scales

In this part, we apply DR to Eq. (II) and discuss the finite-scale nature of the emergent SUSY. The RG equations decouple in sectors which can be studied sequentially, starting with the completely decoupled one of the bosonic and fermionic velocities:

[TABLE]

where $l$ parametrizes the scaling $(\tau,\boldsymbol{x})\rightarrow e^{l}(\tau,\boldsymbol{x})$, $\widetilde{g}=g\widetilde{\mu}^{-\epsilon/2}$ with $\widetilde{\mu}$ independent of $l$ and having the energy unit, see e.g. Ref. [Srednicki, 2007]. Structure of perturbation theory implies that the dimensionful parameters $m$, $b$, $r$ and $a$ cannot appear in Eq. (4) in DR; thus the velocities flow stably toward $v_{f}=v_{b}\equiv v$ for any non-zero $g$, as in the TR invariant case with $r=0$. Grover et al. (2014) By rescaling the spacial coordinate as $\boldsymbol{x}\rightarrow v\boldsymbol{x}$, we can choose $v_{f}=v_{b}=1$ to study the RG flows of other parameters.Nielsen and Ninomiya (1978); Kane and Fisher (1995)

We next consider the RG equations of $g$ and $u$,

[TABLE]

where $\widetilde{u}=u\widetilde{\mu}^{-\epsilon}$. Similar to the velocities, the RG equations of $g$ and $u$ are the same as the TR invariant case at $r=0$, since they are dimensionless in the absence of the dimensional regulator Grover et al. (2014). Thereby, $g$ stably flows toward a non-zero value $g^{*}=\sqrt{16\pi^{2}\epsilon\widetilde{\mu}^{\epsilon}/7}$, while $u$ stably flows toward $u^{*}=3(g^{*})^{2}$. The RG equations for $m$ and $\widetilde{b}=b\widetilde{\mu}^{-\epsilon/2}$ are

[TABLE]

For $g=g^{*}$, the mass anomalous dimension is $3\epsilon/7<1$. Thus, $m$ is relevant and flows toward a value $m^{*}$ determined by the system scale. In contrast to the conventional Ising model without TR-breaking term Goldenfeld (2018), the $\phi^{3}$ term is included in Eq. (II) since it can be generated by the TR-breaking $m$ term as shown below. For a non-zero $m$, it is useful to consider the flow of $b/(mg)$ at the fixed point of the flow of $u$, i.e. $u=3g^{2}$. It is given by

[TABLE]

This indicates that $b$ flows stably toward $b^{*}=3m^{*}g^{*}$, and thus $b$ can be driven away from zero by a non-zero $m$ as long as $g\neq 0$. The RG flow of $(u,b)$, stably flowing toward $(u^{*},b^{*})=(3(g^{*})^{2},3m^{*}g^{*})$, is also verified by numerically plotting the RG flows of $u/g^{2}$ and $b/(gm)$ in Fig. 1(a).

The RG equations of $r$ and $a$ read

[TABLE]

where $\widetilde{a}=a\widetilde{\mu}^{\epsilon/2}$. Though $r$ and $a$ do not have any stable flow as suggested by the above equations, they have unstable fixed points at $r/m^{2}=3/2$ and $ag/m^{3}-r/m^{2}+1=0$, shown by the following RG equations with $(u/g^{2},b/(mg))=(3,3)$:

[TABLE]

as well as Fig. 1(b). In summary, one-loop RG analysis shows that as the scale $l$ increases, the action (II) flows toward

[TABLE]

where the summation over repeated index is implied, $r^{*}$ and $a^{*}$ are the macroscopic values of $r$ and $a$, and $\partial^{2}=\partial_{\tau}^{2}+\boldsymbol{\nabla}^{2}$. Compared with Eq. (II), only $(\frac{a^{*}g^{*}}{(m^{*})^{3}}-\frac{r^{*}}{(m^{*})^{2}})$ needs to be finely tuned to $-1$ in order to achieve SUSY, and $r^{*}<3(m^{*})^{2}/2$ is guaranteed if $r<3(m)^{2}/2$ holds at $l=0$ as suggested by Eq. (III.1), indicating that the blue line in Fig. 1(b) has SUSY. Therefore, the SUSY can emerge as the scale increases after tuning one parameter. (See more details on RG equations in Appendix. C.)

However, unlike quantum critical points, the emergent SUSY of Eq. (II) is only valid in a finite range of scales due to the existence of relevant parameters and ultraviolet cutoff $\Lambda$ of the effective actions. Let us first focus on the SUSY-invaraint version of Eq. (3), since it is the expansion of SUSY-invariant version of Eq. (II) around the SUSY vacua. In $d=4-\epsilon$ dimensions, Eq. (3) has one relevant parameter $m^{\prime}$ when SUSY exists, which exponentially grows as $l$ increases from [math], as shown in Fig. 1(c). When $l$ reaches a critical scale $l_{c}$, $m^{\prime}$ becomes comparable with $\Lambda$; the effective theory fails as $l$ goes beyond $l_{c}$. As a result, the SUSY of the theory is only meaningful when $0<l\leq l_{c}$. However, in experiments and numerical simulations, the system always has finite size, or equivalently finite scale $l$, and thus it is possible to make a sample with scale smaller than $l_{c}$ in order to probe the SUSY signature in this case.

Now we return to Eq. (II) and consider the case where the system is initially away from the SUSY hypersurface (e.g. $b/(mg)$ and $u/g^{2}$ are not exactly 3 at $l=0$) while keeping $\frac{ag}{m^{3}}-\frac{r}{m^{2}}=-1$ and $r<3m^{2}/2$ by parameter tuning. The system still flows to the SUSY hypersurface as the scale $l$ increases before reaching $l_{c}$, as suggested by the RG equations Eq. (4),(III.1) and (III.1). (See Fig. 1(d).) The validity of one-loop RG equations also requires $l<l_{c}$ as discussed in Appendix. C. Reflected in experiments and numerical simulations, the SUSY signature becomes better and better as the sample size increases before reaching the critical scale $l_{c}$. Moreover, the SUSY signature of the sample at the size of $l_{c}$ becomes better if the initial deviation from the SUSY hypersurface decreases. Therefore, we can define a region called “SUSY-emergent region” such that if a system is initially in this region, the low energy physics of the system at the scale $l=l_{c}$ is controlled by the SUSY fixed point suggested by the one-loop RG equations, possibly leading to “exact” SUSY signature within the numerical and experimental errors. (Fig. 1(d)). In this way, Eq. (III.1) can be interpreted as the action of a system in the SUSY-emergent region at a scale $l^{*}$ that is close to the critical scale $l_{c}$. Since $l_{c}$ increases as the initial mass (e.g. $m^{\prime}$ in Eq. (3) at $l=0$) decreases, the SUSY-emergent region expands as the initial mass decreases, and covers the whole parameter space if the initial mass approaches zero, restoring the massless limit. In sum, despite of the existence of the relevant parameters, the signatures of the emergent SUSY can be found for a wide range of scales in a large parameter region (though not completely generic as the massless case).

III.2 Massive GNY model

A concern of the RG analysis done in the last part is that the validity of RG in $d=4-\epsilon$ dimensions might be undermined by the fact that the action (II) does not have 3+1D correspondence since all Pauli matrices have been used for the fermion in 2+1D. Fei et al. (2016) To resolve this issue, let us first introduce the so-called “massive GNY model”. By adding mass-related terms and breaking the Lorentz invairance of the originally massless GNY model Fei et al. (2016), we arrive at the massive GNY model that reads

[TABLE]

where $j=1,...,N_{f}$ is summed over, $\Psi_{j}$ is the four-component Dirac spinor, $N_{f}$ is the number of Dirac fermions, and $\bar{\gamma}^{\mu}$’s are the $4\times 4$ matrices that satisfy $\{\bar{\gamma}^{\mu},\bar{\gamma}^{\nu}\}=2\delta^{\mu\nu}$. Eq. (III.2) has similar form as Eq. (II) if defining $\bar{\psi}=\gamma^{T}\sigma_{y}$, $\psi=\gamma$, $\bar{\gamma}^{0}=\sigma_{y}$ and $\boldsymbol{\bar{\gamma}}=-\mathrm{i}\sigma_{y}\boldsymbol{\alpha}$ for Eq. (II). There are two key differences between Eq. (II) and Eq. (III.2): (i) that Eq. (II) is in 2+1D while Eq. (III.2) is in 3+1D, and (ii) that Eq. (II) only has one complex fermionic degree of freedom while Eq. (III.2) has $N=4N_{f}$ ones. The first difference does not matter since the RG equations of both actions are derived in $d=4-\epsilon$ dimensions, and the second difference can be resolved by formally limiting $N$ to $1$ in the RG equations of Eq. (III.2) as proposed in Ref. [Fei et al., 2016].

The RG equations of velocities in Eq. (III.2) are still decoupled from others and read

[TABLE]

The above equation is exactly the same as Eq. (4) after taking $N\rightarrow 1$, implying the emergence of Lorentz invariance $v_{f}=v_{b}=1$. Similar as the last part, we study the RG flows of other parameters for $v_{f}=v_{b}=1$, and thereby the RG equations of $g$, $u$, $m$, $b$, $r$ and $a$ in Eq. (III.2) take the form

[TABLE]

The above RG equations of $g$ and $u$ coincide with those in Ref. [Fei et al., 2016] at the one-loop level since they are not influenced by the relevant parameters. Formally taking the $N=1$ limit in Eq. (III.2) renders the exact match to (III.1), (III.1) and (III.1), verifying the RG scheme used in the last part.

III.3 Momentum Cutoff

In this part, we revisit the emergent SUSY with the momentum cutoff regularization, i.e. integrating the high-energy modes with spacial momentum $\boldsymbol{k}$ satisfying $\Lambda(1-dl)<|\boldsymbol{k}|<\Lambda$. The finite-scale nature of the emergent SUSY is also reflected in this scheme.

Since this RG scheme does not need the counter-terms, the one-loop Feynman diagrams for this RG scheme are the same as those without counter-terms in Fig. 4. Before deriving the RG equations according to the diagrams, let us first redefine the following quantites: $\widetilde{g}=g\Lambda^{-\epsilon/2}$, $\widetilde{u}=u\Lambda^{-\epsilon}$, $\widetilde{b}=b\Lambda^{-\epsilon/2}$ and $\widetilde{a}=a\Lambda^{\epsilon/2}$. We assume that the relevant parameters are much smaller than $\Lambda$ (with proper power according to the dimension of the quantity), and only keep the zeroth order of $1/\Lambda$ in the RG equations. As a result, the RG equations of the velocity reads

[TABLE]

where $K_{d}=\Omega_{d-1}/(2\pi)^{d-1}$ with $\Omega_{d-1}$ the solid angle in $d-1$ dimensions. Since the above equations are exactly the same as Eq. (4) if setting $K_{d}=1/(2\pi^{2})$, the velocities stably flow to the Lorentz invariant condition $v_{f}=v_{b}=1$. With this condition, the RG equations of $\tilde{g},\tilde{u},m,\tilde{b}$, $r$ and $\tilde{a}$ are

[TABLE]

where $\Delta_{r}=\frac{K_{d}}{8}(2\widetilde{u}-4\widetilde{g}^{2})\Lambda^{2}$ and $\Delta_{a}=\frac{K_{d}}{16}\left(4\widetilde{b}-8\widetilde{g}m\right)\Lambda^{2}$. There are two differences between the above equations and those abtained from DR: (i) the $(2\pi^{2}K_{d})$ factor, and (ii) the extra terms $\Delta_{r}$ and $\Delta_{a}$ in the RG equations of $r$ and $a$. Since the RG equations of $g,u,b,m$ in Eq. (15) have the same form as Eq. (III.1) and (III.1), the stable flow to $(u/g^{2},b/(mg))=(3,3)$ is independent of $r$ and $a$. Despite the extra terms $\Delta_{r}$ and $\Delta_{a}$ in the RG equations of $r$ and $a$, the flow of $(\frac{ag}{m^{3}}-\frac{r}{m^{2}}+1)$ for $(u/g^{2},b/(mg))=(3,3)$ is similar to that in Eq. (III.1), which reads

[TABLE]

indicating the SUSY point is still a fixed point for the RG flow and only needs to finely tune one parameter. As discussed at the beginning of this part, the RG equations (III.3) and (15) are obtained in the condition that the relevant parameters $r,b,m,a$ are small compared with the cutoff $\Lambda$. This condition again reflects the finite-scale nature of the emergent SUSY at the FOQPT discussed in the last section: the emergent SUSY can only be observed before the relevant parameters become comparable with the ultraviolet cutoff as the scale increases.

IV Topological FOQPT

In this section, we demonstrate that SUSY-invariant action (Eq. (II) with Eq. (II)) features a topological FOQPT in the sense that the Majorana field undergoes an unusual topological phase transition, which leads to experimentally testable phenomena. We first describe it in general and then focus on the $a^{*}=0$ case to elaborate the phenomenon.

As discussed in Sec. II, there are two vacua for Eq. (II) when SUSY exists. The two vacua must have the same energy since both vacua have SUSY and SUSY requires the ground state energy to be zero (after removing the constant that is decoupled to the superparnters). Wess and Bagger (1992) This indicates that the SUSY-invariant version of Eq. (II) describes a system right at the FOQPT with the bosonic vacuum expectation value (VEV) $\langle\phi\rangle$ jumping between $(-m\pm\sqrt{3m^{2}-2r})/g$. Here we neglect the possibility that a new vacuum with lower energy appears after including all orders of quantum correction Srednicki (2007). As the fermion mass has the expression $m_{f}=m+g\langle\phi\rangle$, $m_{f}$ takes the values $\pm\sqrt{3m^{2}-2r}$ across the transition. The sign flip of the fermion mass signals a topological phase transition, as the fermion mass domain wall can trap a 1+1D chiral Majorana mode.Wang et al. (2011) More importantly, the non-zero $m_{f}$ across the transition indicates the fermion gap does not close, thus representing a unique topological phase transition without gap-closing owing to its first-order transition nature. Although a similar scenario has been discussed in the literature Amaricci et al. (2015); Roy et al. (2016); Zhu et al. (2018), our case is special because the unchanged mass amplitude $|m_{f}|$ across the transition is required by the emergent SUSY. This feature is better exposed by Eq. (3), which is equivalent to Eq. (II). A suitable rewrite of the bosonic potential $V(\bar{\phi})=\frac{1}{4!}\bar{\phi}^{2}u[\frac{12}{u}(r^{\prime}-\frac{1}{3}\frac{b^{\prime 2}}{u})+(\bar{\phi}+\frac{2b^{\prime}}{u})^{2}]$ shows that the FOQPT now can occur at $r^{\prime}=b^{\prime 2}/(3u)$ and $b^{\prime}\neq 0$, where $\langle\bar{\phi}\rangle$ changes between [math] and $-2b^{\prime}/u$, and $m_{f}$ jumps between $m^{\prime}$ and $m^{\prime}-2b^{\prime}g/u$. Therefore, the relation $b^{\prime}g=um^{\prime}$ imposed by SUSY is essential to maintain the unchanged magnitude of $m_{f}$ and flip its sign across the transition.

To better illustrate the topological FOQPT, we next discuss the action Eq. (III.1) with $a^{*}=0$ as an example. This case is quite general since $a^{*}=0$ can always be achieved by a vacuum shift like that for Eq. (3). We first derive the bosonic VEV $\langle\phi\rangle$ in Eq. (III.1) with $a^{*}=0$, which stands for the macroscopic magnetic ordering, by searching for the global minimum of the bosonic part of the action. Owing to the negative sign in front of $\partial^{2}$, $\langle\phi\rangle$ must be uniform in $(\tau,\boldsymbol{x})$. By minimizing the boson potential $V(\phi)=r^{*}\phi^{2}/2+g^{*}m^{*}\phi^{3}/2+(g^{*})^{2}\phi^{4}/8$, we found a non-magnetic phase (I) and two magnetic phases (II,III) with opposite values of $\langle\phi\rangle$ (See Fig. 2(a)):

[TABLE]

The phases I and II (III) are separated by the transition line $r^{*}=(m^{*})^{2}$ with $m^{*}>0$ ($m^{*}<0$) as depicted by the line 1 (2) in Fig. 2(a), while the Phases II and III are separated by the line $m^{*}=0$ and $r^{*}<0$ (the line 3 in Fig. 2(a)).

According to the above discussion, all the three transition lines are FOQPT lines with emergent SUSY. To show this, we consider a path across the phase transition line 1 or 2, such as the path $\eta_{1}$ in Fig. 2(a). As shown by the blue line in Fig. 2(b), $\left\langle\phi\right\rangle$ vanishes on the phase I side, but approaches $\langle\phi\rangle=-2m^{*}/g^{*}$ as $r^{*}=(m^{*})^{2}+0^{-}$ (on the phase II or III side). Therefore, lines 1 and 2 are FOQPT lines with two degenerate vacua $\langle\phi\rangle=0$ and $-2m^{*}/g^{*}$, where the emergent SUSY was demonstrated above for the former vacuum. Around the latter vacuum, the action for the boson fluctuation $\delta\phi=\phi+2m^{*}/g^{*}$ and the fermion $\gamma$ has the same form as Eq. (III.1) with $a^{*}=0$ and the replacement $m^{*}\rightarrow-m^{*}$ and $r^{*}\rightarrow(m^{*})^{2}$, which also exhibits SUSY. Then, SUSY emerges along the FOQPT lines 1 and 2 in either of the two vacua. Similarly, the line 3 is also a FOQPT line with two degenerate vacua $\langle\phi\rangle=\pm\sqrt{-2r^{*}}/g^{*}$ as shown by the blue line in Fig. 2(c). Around either of the two vacua, the corresponding action for the boson fluctuation $\delta\phi=\phi-(\pm\sqrt{-2r^{*}}/g^{*})$ and the fermion $\gamma$ can be obtained from Eq. (III.1) with $a^{*}=0$ by replacements $m^{*}\rightarrow\pm\sqrt{-2r^{*}}$ and $r^{*}\rightarrow-2r^{*}$, and thus possesses SUSY. We conclude that SUSY exists for all three FOQPT lines in Fig. 2(a) around any of the degenerate vacua. Moreover, the topological feature of the three FOQPT lines is shown by the red lines in Fig. 2(b) and (c), where the fermion mass $m_{f}=m^{*}+g^{*}\langle\phi\rangle$ changes suddenly between $\pm m^{*}$ across the lines 1 and 2, and between $\pm\sqrt{-2r^{*}}$ across the line 3.

V Experimental Setup

In this section, we demonstrate that, by tuning an external magnetic field $B$, the emergent SUSY at FOQPT might be realized on the surface of a TR-invariant TSC with surface magnetic doping, as shown in Fig. 3(a). The action of the TR-invariant TSC reads Roy (2008); Qi et al. (2009); Chung and Zhang (2009)

[TABLE]

where $\bar{\psi}_{\tau,\boldsymbol{k}}=(\bar{\psi}_{\tau,\boldsymbol{k},\uparrow},\bar{\psi}_{\tau,\boldsymbol{k},\downarrow})$ is the Grassman field for the electron with the momentum $\boldsymbol{k}=(k_{x},k_{y},k_{z})$, $h(\boldsymbol{k})=\frac{\boldsymbol{k}^{2}}{2m_{0}}-\mu$ with the chemical potential $\mu$, and $\Delta_{\boldsymbol{k}}=\Delta_{p}(\boldsymbol{k}\cdot\boldsymbol{s})is_{y}$ with the p-wave pairing $\Delta_{p}>0$ and the Pauli matrices $s_{i}$ for spin. Eq. (V) may be used to describe the Ce-based heavy fermion SCs and half-Heusler SCs, and its superfluid version has been realized in B phase of He-3.  Bauer and Sigrist (2012); Savary et al. (2017); Leggett (1975) For $\mu m_{0}>0$, one can solve Eq. (V) with an open boundary condition at $z=0$ and obtain gapless Majorana modes $\gamma$ at the surface Roy (2008); Qi et al. (2009); Chung and Zhang (2009). (See Appendix. E for details.)

The surface magnetic doping of TSC can be phenomenologically described by the standard Ginzburg-Landau free energy of Ising magnetism, $S_{M}=\int d^{d}x[\frac{1}{2}\phi(-\partial_{\tau}^{2}-v_{b}^{2}\boldsymbol{\nabla}^{2}+r_{0})\phi+\frac{1}{4!}u\phi^{4}]$, where $\phi_{\tau,x,y}$ is the order parameter of surface Ising magnetism along $z$. $\phi$ is coupled to electrons at the surface through the exchange interaction, which takes the form $\frac{1}{2}g\int d^{d}x\phi\gamma^{T}\sigma_{y}\gamma$ after the surface projection. Furthermore, a magnetic field $B$ along $z$ is applied and coupled to both electron spin and Ising magnetism on the surface through the Zeeman-type action, and with that, we arrive at the total action:

[TABLE]

(See more details in Appendix. E.) Here, we neglect the orbital effect as all fields are charge neutral, and add an $\phi^{3}$ term since it is allowed by symmetry and can be generated at the quantum level. Since Eq. (V) has exactly the same form as Eq. (II), its RG equations are the same as Eq. (4),(III.1), (III.1) and (III.1), resulting in $a_{2}^{*}=3\mu_{B}^{*}g^{*}$ in addition to $v_{f}=v_{b}=1$ and $u^{*}=3(g^{*})^{2}$. For simplicity, we neglect the $B$-dependence of $\mu_{B}^{*}$, $r_{0}^{*}$ and $a_{1}^{*}$. In this case, as long as $r_{0}^{*}>\frac{3a_{1}^{*}g^{*}}{\mu_{B}^{*}}$, there exists a critical magnatic field $B_{c}=\pm\frac{1}{\mu_{B}^{*}}\sqrt{r_{0}^{*}-\frac{a_{1}^{*}g^{*}}{\mu_{B}^{*}}}$ such that FOQPT with emergent SUSY happens at $B=B_{c}$. To demonstrate this possibility, we choose the values of parameters as $r_{0}^{*}/(g^{*})^{4}=0.2$, $a_{1}^{*}/g^{*}=-1$, $\mu_{B}^{*}(g^{*})^{2}=1$ and $B\geq 0$. As $B$ increases to the critical value $B_{c}\approx 1.1(g^{*})^{4}$, the FOQPT is reached. A signature of emergent SUSY at FOQPT is that the fermion mass $m_{f}=m^{*}+g^{*}\langle\phi\rangle$ should have unchanged magnitude and flip sign across the transition, verified by the red line in Fig. 3(b). The unchanged amplitude of $|m_{f}|$ can be confirmed by local density of states measurement with scanning tunneling microscopy, and the sign flip can be tested by checking the resulting topological phase transition, i.e. the appearance or disappearance of chiral 1+1D domain-wall fermion. Moreover, the surface magnetization $\langle\phi\rangle$ has a sudden change (see the blue line in Fig. 3(b)) as an evidence of FOQPT, which can be measured by superconducting quantum interference devices.

VI Conclusion and Discussion

In conclusion, SUSY with massive superpartners can emerge at the topological FOQPT occurring on the surface of a TR-invariant TSC in a tunable external magnetic field. Although the emergence of the SUSY only happens in a finite range of scales owing to the gapped nature of FOQPT, the scale range can be large when the initial mass is small. Although the emergent SUSY with massive superpartners was proposed in a cold-atom system with spontaneous symmetry breaking Yu and Yang (2010), similar to that at line 3 in Fig. 2(a), that proposal requires to tune more than one parameter and does not discuss the relation between SUSY and the topological FOQPT. Our work also helps shed light upon other emergent symmetries of a FOQPT.Zhao et al. (2019)

VII Acknowledgments

J.Y. thanks Zhen Bi, Shinsei Ryu, Ashvin Vishwanath, Juven Wang and Igor Klebanov for helpful discussions. J.Y. and C.-X.L. acknowledges the support of the Office of Naval Research (Grant No. N00014-18-1-2793), the U.S. Department of Energy (Grant No. DESC0019064) and Kaufman New Initiative research grant KA2018-98553 of the Pittsburgh Foundation. R.R. is supported by the U.S. Department of Energy (Grant No. DE-SC0013699). S.-K. J. is supported by the NSFC under grant 11825404.

Appendix A Rotational Invariance of Eq. (II)

In this section, we show that Eq. (II) of the main text is the most general rotationally invariant action to $\phi\gamma^{T}\sigma_{y}\gamma$ and $\phi^{4}$ order if assuming $\phi$ has uniform classical vacuum.

The rotation transformation along $z$ is defined as $\gamma^{T}_{\tau,\boldsymbol{x}}\rightarrow\gamma^{T}_{\tau,R_{\theta}\boldsymbol{x}}e^{-i\sigma_{y}\theta/2}$ and $\phi_{\tau,\boldsymbol{x}}\rightarrow\phi_{\tau,R_{\theta}\boldsymbol{x}}$, where $R_{\theta}$ is the rotational matrix along $z$ for angle $\theta$ counterclockwise. According to the transformation, the Pauli matrices and the derivatives in the action can be classified according to their angular momentum $L_{z}$, as shown in Tab.1. Based on Tab.1, the most general action without derivatives is

[TABLE]

where Hermiticity requires $g_{n},u_{n}\in\mathds{R}$, and the Hermitian conjugate transformations of the fields are $\gamma_{\tau,\boldsymbol{x}}^{\dagger}=\gamma_{-\tau,\boldsymbol{x}}^{T}$ and $\phi_{\tau,\boldsymbol{x}}^{\dagger}=\phi_{-\tau,\boldsymbol{x}}$ since $\tau$ is imaginary time. In $S_{0}$, $u_{1}$ is the linear $\phi$ term, $g_{0},u_{2}$ are the mass terms, and all other terms are the on-site interaction terms. The most general kinetic term of $\gamma$ with leading order derivatives is

[TABLE]

where $A\in\mathds{R}$. The most general kinetic term of $\phi$ with leading order derivatives is

[TABLE]

where $C_{0},C_{1}\in\mathds{R}$ and there are no first order derivatives since $[\phi_{\tau,\boldsymbol{x}},\phi_{\tau^{\prime},\boldsymbol{x}^{\prime}}]=0$. Then, $S^{\prime}=S_{0}+S_{1}+S_{2}$ is the most general rotationally invariant action if the kinetic term only contains the leading order derivatives and the on-site interaction terms have no derivatives.

Next we show how to derive Eq. (II) of the main text from $S^{\prime}$. Since we want the classical vacuum of $\phi$ to be uniform, we assume $C_{0}>0,C_{1}\geq 0$. We also assume $A\neq 0$ which is typically true for a legitimate fermion action. By changing integration variable $\tau\rightarrow\text{sgn}(A)\tau$ and defining $B_{1}=-\frac{1}{2}|A|v_{f}e^{i\theta_{B}}$, $g_{n}=\widetilde{g}_{n}|A|(C_{0})^{n/2}$, $u_{n}=\widetilde{u}_{n}C_{0}^{n/2}$, $C_{1}=v_{b}^{2}C_{0}$, $\gamma^{T}_{\text{sgn}(A)\tau,\boldsymbol{x}}=\frac{1}{\sqrt{2|A|}}\gamma^{\prime T}_{\tau,\boldsymbol{x}}e^{-i\sigma_{y}\theta_{B}/2}$ and $\phi_{\text{sgn}(A)\tau,\boldsymbol{x}}=\phi^{\prime}_{\tau,\boldsymbol{x}}/\sqrt{2C_{0}}$, $S^{\prime}$ becomes

[TABLE]

where $v_{f}\geq 0$ can be chosen by adjusting $\theta_{B}$, and $\gamma^{\prime}$ and $\phi^{\prime}$ are renamed as $\gamma$ and $\phi$ since they are integrated over in the partition function. The renaming of $\gamma^{\prime}$ and $\phi^{\prime}$ would not cause any physical confusion since $\gamma^{\prime}$ and $\phi^{\prime}$ behave the same as $\gamma$ and $\phi$ under rotation, TR and hermitian conjugate:

[TABLE]

and obvious for $\phi^{\prime}$ and $\phi$ due to their simple relation $\phi^{\prime}_{\tau,\boldsymbol{x}}=\sqrt{2C_{0}}\phi_{\text{sgn}(A)\tau,\boldsymbol{x}}$. If we only keep the on-site interaction terms up to the $\phi\gamma^{T}\sigma_{y}\gamma$ and $\phi^{4}$ order and define $\widetilde{g}_{0}=m$,$\widetilde{g}_{1}=\sqrt{2}g$, $\widetilde{u}_{1}=\sqrt{2}a$, $\widetilde{u}_{2}=r$, $\widetilde{u}_{3}=\sqrt{2}b/3$ and $\widetilde{u}_{4}=u/3!$, Eq. (A) is the same as Eq. (II) . Therefore, Eq. (II) is the most general rotational invariant action if (i) only keeping leading order derivatives in the kinetic terms, (ii) neglecting derivatives in the on-site interaction, (iii) only keeping terms to $\phi\gamma^{T}\sigma_{y}\gamma$ and $\phi^{4}$ order for the on-site interaction, and (iv) assuming the classical vacuum of $\phi$ is uniform.

At last, we show Eq. (II) is in the most general rotationally invariant form if only keeping terms to $\phi\gamma^{T}\sigma_{y}\gamma$ and $\phi^{4}$ order for the on-site interaction and assuming the classical vacuum of $\phi$ is uniform, which is the statement at the beginning of this section. It means that we need to argue why higher-order derivatives in the kinetic energy terms and derivatives in the on-site interaction can be neglected. The argument will be done by dimension analysis. The dimensions of fields are $[\gamma]=(d-1)/2$ and $[\phi]=(d-2)/2$. As a result, we have $[v_{f}]=[v_{b}]=0$, $[\widetilde{u}_{n}]=(2-n)d/2+n$ and $[\widetilde{g}_{n}]=-nd/2+n+1$. The zero dimension of $v_{f}$ and $v_{b}$ means any higher-order derivatives in the kinetic terms of $\gamma$ and $\phi$ are irrelevant, and thus can be neglected. Now we discuss the on-site interaction. Since we perform RG analysis in $d=4-\epsilon$, we may analyze the dimension of the interaction for $d=4$, resulting that $[\widetilde{u}_{n}]=4-n$ and $[\widetilde{g}_{n}]=1-n$. It means, keeping $\phi\gamma^{T}\sigma_{y}\gamma$ and $\phi^{4}$ order is equivalent to neglect all the irrelevant on-site interaction terms at $d=4$. The remaining on-site interaction terms include $b$, $g$ and $u$. $g$ and $u$ are marginal, and thus adding derivatives to the two terms would make them irrelevant. $[b]=1$ for $d=4$ and thus allows one derivative. However, this derivative would be a total derivative since $(\partial_{x}\phi)\phi^{2}=(\partial_{x}\phi^{3})/3$ and can be neglected. Therefore, the derivatives in the on-site interaction can be neglected if only keeping terms to the $\phi\gamma^{T}\sigma_{y}\gamma$ and $\phi^{4}$ order. The statement at the beginning of this section is proven.

Appendix B SUSY

In this section, we discuss the SUSY algebra and transformation. Since the fields considered here are one two-component Majorana field and one Ising field, there should be two supercharges $Q_{a}$ with $a=1,2$ instead of four for $N=1$ Wess-Zumino model in 3+1D. Wess and Bagger (1992) The supercharges satisfy

[TABLE]

where $\bar{\alpha}^{\mu}=\sigma_{y}(\alpha^{\mu})^{T}\sigma_{y}$ and $P_{\mu}$ is the energy-momentum operator that gives $[P_{\mu},\varphi]=i\partial_{\mu}\varphi$ for any field operator $\varphi$. Here the metric is chosen as $(-,+,+)$. The infinitesimal SUSY transformation is defined as

[TABLE]

which gives

[TABLE]

Here $\xi,\eta$ are two two-component Grassmann numbers. Eq. (26) shows the closure of the SUSY algebra.

Now, we show the SUSY of SUSY-invaraint version of Eq. (3). For simplicity, we replace $m^{\prime}$ and $\bar{\phi}$ by $m$ and $\phi$, respectively. The action can be re-written as

[TABLE]

where

[TABLE]

The SUSY transformation in this case can be defined as

[TABLE]

To demonstrate $S^{*}$ is SUSY invariant, we can first act $\delta_{\xi}$ on $S_{\gamma,0}$ and get

[TABLE]

Then act $\delta_{\xi}$ on $S_{\gamma\phi}$ and get

[TABLE]

As a result, we have

[TABLE]

The defined SUSY transformation must be close, i.e. satisfying Eq. (26). It is true for $\phi$, i.e.

[TABLE]

For $\gamma$, we have

[TABLE]

where $\eta\xi^{T}-\xi\eta^{T}=\sum_{\mu}\eta^{T}\alpha^{\mu}\xi\alpha^{\mu}$ and $\sigma_{y}\alpha^{\mu}\alpha^{\nu}+\alpha^{\nu}\alpha^{\mu}\sigma_{y}=2\sigma_{y}\delta^{\mu\nu}$ are used. Clearly, the closure of the algebra for $\gamma$ requires the equatin of motion of $\gamma$, which is $(i\alpha^{\mu}\partial_{\mu}+m\sigma_{y}+g\phi\sigma_{y})\gamma=0$. We call the algebra is close through the equation of motion. The requirement of the equation of motion is because we integrate out the auxiliary field.Wess and Bagger (1992)

Appendix C Details For RG Equations

In this section, we derive Eq. (4),(III.1), (III.1) and (III.1). We first derive the Callan-Symanzik equation for $N$-point function, and then show the RG equations to the one-loop order.

C.1 Callan-Symanzik Equation

The regularization scheme chosen here is the dimensional regularization with $d=4-\epsilon$. For a generic dimension $d$, the dimensions of the fields, velocities, masses and interaction couplings are $[\gamma]=(d-1)/2$, $[\phi]=(d-2)/2$, $[v_{b}]=[v_{f}]=0$, $2[m]=[r]=2$, $[a]=1+d/2$, $[g]=(4-d)/2$, $[b]=(6-d)/2$ and $[u]=4-d$. For $d=4$, $[a]=3$, $[g]=[u]=0$ and $[b]=1$. In order to keep the dimensions of the couplings the same as $d=4$ in the $d=4-\epsilon$ scheme, we introduce a parameter $\widetilde{\mu}$ with $[\widetilde{\mu}]=1$ by doing the transformation

[TABLE]

In addition, in order to include the quantum corrections, we should introduce $Z$ factors to the action, and the action becomes

[TABLE]

where the partition function is $Z=\int D\gamma D\phi e^{-S}$ and $X=(\tau,\boldsymbol{x})$. Here, $Z$ factors are chosen to only cancel the divergent part of the quantum corrections ($MS$ schemeSrednicki (2007)). Since the divergence of the quantum corrections is given by $1/\epsilon^{n}$ with $n$ positive integer, the $Z$ factors must have the form

[TABLE]

with $i=\gamma,\phi,v_{f},v_{b},m,r,a,g,b,u$ and $``\ln"$ function formally defined as $\ln(1+x)=\sum_{n=1}^{+\infty}\frac{(-1)^{n-1}x^{n}}{n}$. Now, we derive the expressions $\beta$ and $\Gamma$ functions, where we use $\Gamma$ instead of commonly used $\gamma$ to label the Gamma functions since the latter is reserved for the Majorana field. Since $\widetilde{\mu}$ is not physical, the physical action as well as any physical observable should not depend on $\widetilde{\mu}$. By defining

[TABLE]

the partition function becomes $Z=C_{0}\int D\gamma^{(0)}D\phi^{(0)}e^{-S}$ with $S$ becomes

[TABLE]

In this way, we obtain a partition function with an action that is independent of $\widetilde{\mu}$. Since $\widetilde{\mu}$ is not physical, $S$ in the above equation should be the physical action that gives the physical observables. It means the fields $\phi^{(0)}$ and $\gamma^{(0)}$ and the parameters $w^{(0)}$’s in the above equation are physical and independent of $\widetilde{\mu}$, where $w=v_{f},v_{b},r,m,a,g,b,u$. As a result, we have

[TABLE]

where $D=\frac{g}{2}\partial_{g}+\frac{b}{2}\partial_{b}+u\partial_{u}-\frac{a}{2}\partial_{a}$ and the condition that $\beta$ and $\Gamma$ functions are finite at $\epsilon\rightarrow 0$ is used. Now we derive the Callan-Symanzik equation. The physical $N$-point function is defined as

[TABLE]

with $X_{n}$ indicating all its coordinate dependence. The above equation is related with the $N$-point function of $\gamma$ and $\phi$ by Eq. (C.1):

[TABLE]

Combining the fact that $F^{(0)}_{N_{\phi},N_{\gamma}}$ is independent of $\widetilde{\mu}$ and Eq. (C.1), we arrive at the Callan-Symanzik equation with respect to $\widetilde{\mu}$:

[TABLE]

But we want the Callan-Symanzik equation with respect to the physical scale instead of the non-physical $\widetilde{\mu}$. To do so, consider the $N$-point function with scaled coordinates $F_{N_{\phi},N_{\gamma}}(tX_{n},w,\widetilde{\mu})$. By defining $\gamma_{X}^{\prime}=\gamma_{tX}t^{[\gamma]}$, $\phi_{X}^{\prime}=\phi_{tX}t^{[\phi]}$, $w^{\prime}=wt^{[w]}$ and $\widetilde{\mu}^{\prime}=\widetilde{\mu}t$ , we have

[TABLE]

Differentiating the above equation by $t$ at $t\rightarrow 1$ and using Eq. (C.1), we have

[TABLE]

The meaning of the above equation can be better illutrated in the integrated form:

[TABLE]

where

[TABLE]

and $\zeta_{\phi}(0)=\zeta_{\gamma}(0)=1$. Eq. (C.1) indicates that the $N$-point function at a larger scale $e^{l}X_{n}$ is the same as the $N$-point function at $X_{n}$ with scaled fields $\phi,\gamma$ and parameters $w$ according to Eq. (C.1). It means, we can define an action $S^{l}$ that is the same as $S$ except that the parameters $w$ in $S^{l}$ can scale with $l$ according to Eq. (C.1). In this case, the N-point function generated by $S^{l}$ at $X_{n}$ and the N-point function given by $S$ at $e^{l}X_{n}$ just deviate from each other by a factor $[\zeta_{\phi}(l)]^{N_{\phi}}[\zeta_{\gamma}(l)]^{N_{\gamma}}$. Then, $S_{l}$ can be viewed as an effective action of $S$ at a larger scale $e^{l}$. And Eq. (C.1) is called the RG equations. In the following, we will derive the RG equations to the one-loop order.

C.2 One-Loop RG Equations

In order to obtain Eq. (4),(III.1), (III.1) and (III.1), we need to find the corresponding $\beta$ and $\Gamma$ functions in Eq. (C.1). According to Eq. (C.1) and Eq. (37), it is equivalent to deriving the expression of the $Z$ factors in Eq. (C.1). Since the $Z$ factors are given by the quantum corrections, we need to evaluate the loop diagrams. For convenience, we seperate Eq. (C.1) into two parts $S=S_{c}+S_{ct}$ with

[TABLE]

and

[TABLE]

Here $S_{c}$ is Eq. (C.1) without quantum corrections and $S_{ct}$ is called the counter-terms. According to Eq. (C.2), the fermion propagator is

[TABLE]

and the boson propagator reads

[TABLE]

Here $k=(\omega,\boldsymbol{k})$ and $q=(\nu,\boldsymbol{q})$. Equiped with Eq. (C.2)-(51), we can evaluate the loop diagrams. For simplicity, we only consider the one-loop diagrams, which together with the corresponding counter-terms are shown in Fig. 4.

The one-loop contribution to the fermion self-energy $\Sigma(k)$ is given by Fig. 4 a and b, which reads

[TABLE]

where we only keep the divergent terms since $Z$ factors are chosen to only cancel the divergent part.

The one-loop contribution to the boson self-energy $\Pi(q)$ is given by Fig. 4 c, d, e and f, which reads

[TABLE]

The fermion-boson coupling term $\Gamma^{\phi\gamma\gamma}$ has the one-loop correction given by Fig. 4 g,h and i, which reads

[TABLE]

where the contribution of Fig. 4h is not explicitly included since it is not divergent for $d=4-\epsilon$.

For the $\phi^{3}$ term, the one-loop contribution is given by Fig. 4 j, k, l and m, and reads

[TABLE]

where Fig. 4l is not divergent and included in $O(\epsilon^{0})$.

Fig. 4 n, o, p, q and r give the one-loop contribution of the $\phi^{4}$ term. However, only Fig. 4 n, o and r are divergent and contribute to the RG equations, which gives

[TABLE]

Fig. 4 s, t and u give the one-loop contribution of the $\phi$ term, which gives

[TABLE]

According to Eq. (C.2)-(C.2) and Eq. (37), we have the expressions of $A^{(1)}_{a}$’s as

[TABLE]

Plug the above equation into Eq. (C.1), we can get the $\beta$ and $\Gamma$ functions and then we can derive the RG equations according to Eq. (C.1). Due to the transformation Eq. (35), $a$, $g$, $u$ and $b$ in the obtained RG equations should be replaced by $\widetilde{a}$, $\widetilde{g}$, $\widetilde{u}$ and $\widetilde{b}$ according to the convention defined in the main text. As a result, we can get Eq. (4),(III.1), (III.1) and (III.1) in the main text.

From Eq. (4),(III.1), (III.1) and (III.1) in the main text, we discuss the RG equations of $w=v_{f}/v_{b}$ and the RG equation of $u/g^{2}$. The RG equation of $w$ reads

[TABLE]

which can be re-written as

[TABLE]

with $\widetilde{w}=w-1$. The RG equation of $\widetilde{w}$ shows that the velocity ratio $v_{f}/v_{b}$ flows exponentially to 1 when it is close to 1. On the other hand, the RG equation of $u/g^{2}$ for $v_{f}=v_{b}=1$ reads

[TABLE]

Define $y=u/g^{2}-3$, then the above equation can be re-written as

[TABLE]

where the existence of y-linear term indicates the exponential flow of $u/g^{2}$ to $3$. Moreover, according to Eq. (5) of the main text, the RG equation of $z=b/(mg)-3$ for $v_{f}=v_{b}=1$ and $u/g^{2}=3$ is linear in $z$:

[TABLE]

indicating the exponential flow of $b/(mg)$ to $3$. Therefore, although the classical dimensions of $v_{f}/v_{b}$, $u/g^{2}$ and $b/(mg)$ are zero, their stable flows to the SUSY hypersurface are exponential instead of logarithmic after including the quantum correction.

C.3 Higher-loop Contribution

Unlike the continuous phase transition, the validity of perturbation theory in series of loops is not obvious in our case due to the existence of relevant $m$, $r$ and $b$, which might make the higher-loop terms not small. In this part, we discuss the validity of neglecting the higher-loop terms. Here we still replace the $\tilde{a},\tilde{g},\tilde{u},\tilde{b}$ in the main text by $a,g,u,b$ as above, and we only consider the 1PI connected graphs with loops, which do not contain $a$.

As shown above, the RG equations are obtained from the counter-terms. The counter-terms are determined by the graphs with non-negative superficial degree of divergence $D=4L-2I_{\phi}-I_{\gamma}$, where $L$ is the number of loops in the graph, and $I_{\phi}$ and $I_{\gamma}$ are the numbers of internal bosonic and fermionic propagators, respectively.Srednicki (2007) To estimate the order of higher-loop terms, we need to know the structure of a generic connected graph. As one end of the external propagator and the two ends of the internal propagator are connected to vertexes, we have $E_{\phi}+2I_{\phi}=V_{g}+3V_{b}+4V_{u}$ and $E_{\gamma}+2I_{\gamma}=2V_{g}$, where $E_{\phi}$ and $E_{\gamma}$ are the numbers of external bosonic and fermionic propagators, respectively, and $V_{g}$, $V_{b}$ and $V_{u}$ are number of $g$, $b$ and $u$ vertexes, respectively. Furthermore, the momentum conservation gives $L=I_{\phi}+I_{\gamma}-(V_{g}+V_{b}+V_{u})+1$. From the above relations, we have

[TABLE]

Equipped with those relations, we next address the possible issues brought by the three relevant parameters by estimating the order of the $L$-loop contribution when $a\sim m^{3}/g$, $u\sim g^{2}$, $b\sim gm$, $r\sim m^{2}$ and $v_{f}\sim v_{b}\sim 1$, which is near the SUSY hypersurface. With the assumption, we also have $m^{\prime}\sim m$.

Since only $D\geq 0$ graphs contribute to counter-terms, high powers of $b$ do not exist in the RG equations as they can make $D$ negative according to Eq. (C.3). Therefore, the relevant $b$ cannot cause any divergence in summing the series. On the other hand, $D\geq 0$ requires $2E_{\phi}+3E_{\gamma}\leq 8$, leading to the following six combinations: $(E_{\phi},E_{\gamma})=(1,0)$, $(2,0)$, $(3,0)$, $(4,0)$, $(0,2)$ or $(1,2)$. Here we use the fact that $E_{\gamma}$ can only be even. Clearly, those six combinations correspond exactly to those counter-terms that we include in Eq. (C.1), verifying the renormalizability of our theory. In the following, we estimate the order of the counter-terms for the six combinations, namely the $Z$ factors. A generic connected graph with loop has the form $g^{V_{g}}b^{V_{b}}u^{V_{u}}F(k)$ with the loop-integral part $F(k)$ having dimension $[F(k)]=D$, where $\widetilde{\mu}$ is omitted since $\epsilon\rightarrow 0$ and will be taken into account in the $1/\epsilon$ expansion in the following. For the condition that we choose, only $m$ and $k$ can carry dimensions, and thus the graph is approximately $g^{2L-2+E_{\phi}+E_{\gamma}}\sum_{n}C_{n}m^{D+V_{b}-n}k^{n}$, where $C_{n}$ is dimensionless and given by the loop integral. Based on this estimation, all the $Z$ factors approximately have the form $1+\sum_{L=1}^{\infty}C_{L}(g^{2})^{L}$, where the dimensionless $C_{L}$ has the form $\sum_{i=1}^{L}C^{(i)}_{L}/\epsilon^{i}$ and only depends on $\ln(|m|/\tilde{\mu})$. Since only $C^{(1)}_{L}$ contributes to the RG equations, we only need to care about whether $\sum_{L=1}^{\infty}C_{L}^{(1)}(g^{2})^{L}$ is well defined. As $m$ is relevant, $|m|\gg\tilde{\mu}$ is commonly true at a relatively large scale, and thus the $C^{(1)}_{L}\sim[\ln(|m|/\tilde{\mu})]^{L-1}$. Then, the $L$-loop term in the series of interest is of the order $\sim[\ln(|m|/\tilde{\mu})]^{L-1}(g^{2})^{L}$. Therefore, the higher-loop contributions to the RG equations are neglectable if $\ln(|m|/\tilde{\mu})<1/g^{2}\sim 1/\epsilon\sim\ln(\Lambda/\tilde{\mu})$, where $\Lambda$ is the ultraviolet energy cut-off of the model. Furthermore, it means $\ln(|m^{\prime}|/\tilde{\mu})<\ln(\Lambda/\tilde{\mu})$. Clearly, in the assumption that the ultraviolet energy cut-off is much larger than any dimensionful parameter of the model (or similarly $l<l_{c}$), the higher-loop terms can be neglected and the one-loop result is trustworthy.

C.4 Emergent-SUSY Region

In this part, we discuss the example of paths 1 and 2 of Fig. 1d, as shown in Fig. 5(a) and (b), respectively. Fig. 5 is ploted according to Eq. (III.1), Eq. (III.1) and Eq. (III.1) with $a=\frac{m}{g}(r-m^{2})$ fixed. $m^{\prime}$ is determined by the relation for Eq. (3), and the critical scale is defined as $\ln(|m^{\prime}/\widetilde{\mu}|)=\frac{1}{5}\frac{1}{\epsilon}$ since $1/\epsilon\sim\log(\Lambda/\widetilde{\mu})$. The values of parameters are $\widetilde{g}_{i}=0.8,\epsilon=10^{-3},m_{i}/\widetilde{\mu}=0.1,r_{i}/\widetilde{\mu}^{2}=0.1$ for both subgraphs of Fig. 5, while $(u_{i}/g_{i}^{2},b_{i}/(m_{i}g_{i}))=(0.4,0.3)$ for (a) and $(u_{i}/g_{i}^{2},b_{i}/(m_{i}g_{i}))=(2.4,2.3)$ for (b). Here the subscript “$i$” means they are initial (at $l=0$) values. According to Fig. 5, path 1 cannot reach the SUSY point before $l$ reaches $l_{c}$, while path 2 can be very close to the SUSY point.

Appendix D First-Order Phase Transitions

In this section, we show details on determining the phases of Eq. (III.1) with $a^{*}=0$. For convenience, we drop the “$*$” of all parameters in Eq. (III.1).

$\langle\phi\rangle$ is given by the global minimum of bosonic part of Eq. (III.1) with $a^{*}=0$, which reads

[TABLE]

Since $\langle\phi\rangle$ indicates the macroscopic magnetic ordering and must be real, we should impose $\phi^{\dagger}_{\tau,\boldsymbol{x}}=\phi_{\tau,\boldsymbol{x}}$ when solving for the global minimum of $S^{*}_{b}$. $S_{b}^{*}$ can be splitted into two parts $S_{b}^{*}=S_{T}+S_{V}$ with

[TABLE]

Here $V(\phi)=\frac{1}{2}r\phi^{2}+\frac{1}{2}gm\phi^{3}+\frac{1}{8}g^{2}\phi^{4}$. By defining $\phi_{\tau,\boldsymbol{x}}=\frac{1}{\beta\mathcal{V}}\sum_{\omega,\boldsymbol{k}}e^{-i\omega\tau+i\boldsymbol{k}\cdot\boldsymbol{x}}\phi_{\omega,\boldsymbol{k}}$ and $\phi_{\omega,\boldsymbol{k}}^{*}=\phi_{-\omega,-\boldsymbol{k}}$, $S_{T}$ can be re-written as

[TABLE]

It means $S_{T}\geq 0$ and $S_{T}=0$ holds only if $\phi$ is uniform in $(\tau,\boldsymbol{x})$. On the other hand, for any non-uniform field $\phi_{\tau,\boldsymbol{x}}$, we can always find a uniform field $\phi^{(0)}$ such that $S_{V}[\phi^{(0)}]\leq S_{V}[\phi]$. (Simply, one can pick one position $(\tau_{0},\boldsymbol{x}_{0})$ such that $V(\phi_{\tau,\boldsymbol{x}})\geq V(\phi_{\tau_{0},\boldsymbol{x}_{0}})$ holds for any $(\tau,\boldsymbol{x})$ and define the uniform field as $\phi^{(0)}=\phi_{\tau_{0},\boldsymbol{x}_{0}}$.) Therefore, the global minimum of $S^{*}_{b}$ must be uniform given $\phi^{\dagger}_{\tau,\boldsymbol{x}}=\phi_{\tau,\boldsymbol{x}}$. Therefore, we only need to minimize the bosonic potential $V(\phi)$ to obtain the minimum field. The extrema of $V(\phi)$ are given by

[TABLE]

In this case, if $9m^{2}<8r$, there is only one extremum that is $\phi_{0}=0$, and if $9m^{2}\geq 8r$, the extrema are

[TABLE]

Clearly, $\phi_{\pm}\neq 0$ as long as $r\neq 0$ and they exist. The values of $V(\phi)$ at the three extrema are

[TABLE]

Now we discuss the global minimum of $V(\phi)$. As mentioned before, if $9m^{2}<8r$, $\phi_{0}$ is only one extremum and thereby is the global minimum. If $9m^{2}/8\geq r>m^{2}$, we have $4r-3m^{2}>m^{2}>|m\sqrt{9m^{2}-8r}|$ , which gives $2gm\phi_{\pm}+4r=4r-3m^{2}\pm m\sqrt{9m^{2}-8r}>0$ and thus $V(\phi_{\pm})>0$. Therefore, $\phi_{0}=0$ is the global minimum if $r>m^{2}$. For $r<m^{2}$ and $m>0$, we have

[TABLE]

which gives $V(\phi_{-})<0$. $m>0$ results in $2gm\phi_{-}+4r<2gm\phi_{+}+4r$ and $\phi_{-}^{2}>\phi_{+}^{2}$. In this case, if $2gm\phi_{+}+4r\geq 0$, $V(\phi_{+})\geq 0>V(\phi_{-})$. And if $2gm\phi_{+}-4r<0$, $\phi_{-}^{2}>\phi_{+}^{2}$ gives $V(\phi_{+})>V(\phi_{-})$. Then, we have $V(\phi_{-})<V(\phi_{+})$. Therefore, $\phi_{-}$ is the global minimum if $r<m^{2}$ and $m<0$. Since $V(\phi)$ is invariant under $m\rightarrow-m$ and $\phi\rightarrow-\phi$, the global minimum is at $\phi_{+}$ if $r<m^{2}$ and $m>0$. Therefore, $r=m^{2},m\neq 0$ and $r<0,m=0$ are where the first order phase transition happens. In the following, we show the form of the action at the first order phase transition.

At $r=m^{2}$ with $m\neq 0$, there are two degenerate vacua: (i) $\langle\phi\rangle=0$, and (ii) $\langle\phi\rangle=-2m/g$. Around the first vacuum, the form of the action is just $S^{*}$ with $r=m^{2}$ and $a=0$. Around the second vacuum, we need to use the fluctuation $\delta\phi=\phi-(-2m/g)$ in $S^{*}$ with $a=0$, and the resulted action for $\gamma$ and $\delta\phi$ reads

[TABLE]

At $r<0$ and $m=0$, $S^{*}$ reads

[TABLE]

Since $r<0$, the boson part has spontaneous symmetry breaking(SSB) and the new vacua are at $\phi=\pm\sqrt{-2r}/g$. Around the new vacua $\phi=\delta\phi\pm\sqrt{-2r}/g$, we have

[TABLE]

Appendix E Details on the Experimental Setup

In this section, we derive Eq. (V) from Eq. (V).

To solve for the surface modes, we consider a semi-infinite configuration ($z<0$) with open boundary condition at $z=0$. Then, we address the $z$ direction in the real space with $\bar{\psi}_{\tau,\boldsymbol{k}}=\int dze^{ik_{z}z}\bar{\psi}_{\tau,\boldsymbol{k}_{\shortparallel},z}$, and Eq. (V) becomes

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

$\bar{\mu}=\mu-\boldsymbol{k}_{\shortparallel}^{2}/(2m_{0})$, $\tau_{i}$’s are Pauli matrices for the particle-hole index and $\bar{\Psi}_{\tau,\boldsymbol{k}_{\shortparallel},z}=(\bar{\psi}_{\tau,\boldsymbol{k}_{\shortparallel},z},\psi^{T}_{\tau,-\boldsymbol{k}_{\shortparallel},z})$. In the following, we first solve for the zero modes of $h_{0}$ and then add $h_{1}$ as a perturbation since $\boldsymbol{k}$ is small. The zero mode equation for $h_{0}$ reads

[TABLE]

Define four orthonormal vectors $\xi_{a,j}$’s as

[TABLE]

where $a=\pm$, $j=1,2$ and $\tau_{y}s_{x}\xi_{a,j}=a\xi_{a,j}$. The wave function can be re-expressed as $\Phi_{z}=\sum_{a,j}f_{a,j}(z)\xi_{a,j}$, and Eq. (79) is equivalent to

[TABLE]

for all $a=\pm$ and $j=1,2$ with boundary condition $f_{a,j}(0)=f_{a,j}(-\infty)=0$. Without loss of generality, we choose $\Delta_{p}m_{0}>0$. In this case, we only have solution for $a=+$ and $2m_{0}\mu-\boldsymbol{k}_{\shortparallel}^{2}>0$, which is $f_{+,j}(\boldsymbol{k}_{\shortparallel}^{2},z)=f_{+}(\boldsymbol{k}_{\shortparallel}^{2},z)=C_{0}e^{m_{0}\Delta_{p}z}\sinh(z\sqrt{\Delta_{p}^{2}m_{0}^{2}-2\bar{\mu}m_{0}})$ with $C_{0}$ the normalization constant that makes $f_{+}$ real. The wave function of the zero modes in general has the form

[TABLE]

Clearly, there are two independent zero modes, of which the wavefunction can be chosen as

[TABLE]

To get the low-energy effetive model, we define $\bar{\gamma}_{\tau,\boldsymbol{k}_{\shortparallel},i}=\int^{0}_{-\infty}dz\bar{\Psi}_{\boldsymbol{k}_{\shortparallel},z}\Phi_{i,\boldsymbol{k}_{\shortparallel}^{2},z}$. Using $\bar{\Psi}_{\tau,\boldsymbol{k}_{\shortparallel},z}=\sum_{i}\bar{\gamma}_{\tau,\boldsymbol{k}_{\shortparallel},i}\Phi^{\dagger}_{i,\boldsymbol{k}_{\shortparallel}^{2},z}+...$ with “$...$” the high energy contribution, we can project $h_{1}$ to the zero modes and get

[TABLE]

where $\bar{\gamma}_{\tau,\boldsymbol{k}_{\shortparallel}}=(\bar{\gamma}_{\tau,\boldsymbol{k}_{\shortparallel},1},\bar{\gamma}_{\tau,\boldsymbol{k}_{\shortparallel},2})$ and the terms related with the high-energy modes are neglected here. Define $\gamma_{\tau,\boldsymbol{k}_{\shortparallel}}=\int d^{2}xe^{-i\boldsymbol{k}_{\shortparallel}\cdot\boldsymbol{x}}\gamma_{\tau,\boldsymbol{x}}$. Due to $\tau_{x}\Phi_{i,\boldsymbol{k}_{\shortparallel}^{2},z}^{*}=\Phi_{i,\boldsymbol{k}_{\shortparallel}^{2},z}$, $i\tau_{0}s_{y}\Phi_{i,\boldsymbol{k}_{\shortparallel}^{2},z}^{*}=\sum_{j}\Phi_{j,\boldsymbol{k}_{\shortparallel}^{2},z}(i\sigma_{y})_{ji}$, $\Psi^{T}_{\tau,-\boldsymbol{k}_{\shortparallel},z}=\bar{\Psi}_{\tau,\boldsymbol{k}_{\shortparallel},z}\tau_{x}$ and $\bar{\Psi}_{\tau,\boldsymbol{k}_{\shortparallel},z}\xrightarrow{\text{TR}}\Psi_{\tau,-\boldsymbol{k}_{\shortparallel},z}(i\tau_{0}s_{y})$, we have $\bar{\gamma}_{\tau,\boldsymbol{x}}=\gamma_{\tau,\boldsymbol{x}}^{T}$ and $\gamma_{\tau,\boldsymbol{x}}^{T}\xrightarrow{\text{TR}}\gamma_{\tau,\boldsymbol{x}}^{T}(i\sigma_{y})$. As a result, the above action becomes

[TABLE]

Replacing $\Delta_{p}$ by $v_{f}$ in the above action, we can get fermionic part of the Eq. (V). Note that we neglect the high-energy modes, and thus $\Delta_{p}=v_{f}$ holds only to the leading order.

The surface magnetic doping of TSC can be phenomenologically described by the standard Ginzburg-Landau free energy of Ising magnetism, $S_{M}=\int d^{3}x[\frac{1}{2}\phi(-\partial_{\tau}^{2}-v_{b}^{2}\boldsymbol{\nabla}^{2}+r_{0})\phi+\frac{1}{4!}u\phi^{4}]$, where $\phi_{\tau,x,y}$ is the order parameter of surface Ising magnetism along $z$. $\phi$ is coupled to electrons at the surface through the exchange interaction $S_{ex}=\int d\tau\int d^{3}r\phi_{\tau,x,y}g_{M}(z)\bar{\psi}_{\tau,\boldsymbol{r}}s_{z}\psi_{\tau,\boldsymbol{r}}$ with $g_{M}(z)$ localized at the surface. The matrix form of the exchange interaction in the BdG bases is $\tau_{z}s_{z}/2$, of which the projection to the Majorana modes is $\sigma_{y}/2$. And thus the Ising coupling after the surface projection takes the form $\frac{1}{2}g\int d^{3}x\phi\gamma^{T}\sigma_{y}\gamma$ with $g=\int dzf_{+}(0,z)^{2}g_{M}(z)$. Here we neglect the momentum dependence in $f_{+}(\boldsymbol{k}_{\shortparallel}^{2},z)$ as $\boldsymbol{k}_{\shortparallel}$ is small. Now the total action should read

[TABLE]

which has TR symmetry. Furthermore, a magnetic field $B$ along $z$ is applied and coupled to both electron spin and Ising magnetism on the surface through the Zeeman-type action: $S_{B}=\int d\tau d^{3}rg_{s}(z)B\bar{\psi}_{\tau,\boldsymbol{r}}s_{z}\psi_{\tau,\boldsymbol{r}}+\int d^{3}x(a_{1}B\phi+\frac{a_{2}}{3!}B\phi^{3})$ with $g_{s}(z)$ localized at the surface. Here, we neglect the orbital effect as all fields are charge neutral, and add an $\phi^{3}$ term since it is allowed by symmetry and can be generated at the quantum level. After projecting the $g_{s}$ term to the surface in the same way as the Ising coupling, we can include $S_{B}$ into the total action, and get Eq. (V) with $\mu_{B}=\int dzf_{+}(0,z)^{2}g_{s}(z)$.