Interacting Majorana fermions in strained nodal superconductors

TL;DR

This paper explores how non-uniform strain induces Landau levels and Majorana fermions in 2D and 3D superconductors with Dirac nodal points, revealing topological phase transitions and potential for exotic quasiparticles.

Contribution

It demonstrates the spinful Majorana nature of zeroth Landau levels in strained nodal superconductors and identifies topological phases with localized Majorana modes at interfaces.

Findings

Zeroth Landau levels exhibit Majorana fermion characteristics.

Strain induces topological phase transitions in superconductors.

Localized Majorana modes appear at interfaces under certain strains.

Abstract

Landau levels (LL) have been predicted to emerge in systems with Dirac nodal points under applied non-uniform strain. We consider 2D, $d_{xy}$ singlet (2D-S) and 3D $p \pm i p$ equal-spin triplet (3D-T) superconductors (SCs). We demonstrate the spinful Majorana nature of the bulk gapless zeroth-LLs. Strain along certain directions can induce two topologically distinct phases in the bulk, with zeroth LLs localized at the the interface. These modes are unstable toward ferromagnetism for 2D-S cases. Emergent real-space Majorana fermions in 3D-T allow for more exotic possibilities.

Equations

$${\gamma }_{s0}(\mathbf{k})=\sum _{s^{\prime }}{M}_{s{s}^{\prime }}{\gamma }_{s^{\prime }0}^{†}(-\mathbf{k}),$$

$$H=\sum _{\mathbf{k}}{\mathbf{\Psi }}^{†}(\mathbf{k})\left[\hat{h}(\mathbf{k})+\hat{\Delta }(\mathbf{k})\right]\mathbf{\Psi }(\mathbf{k}),$$

$$c_{P,\sigma }(\mathbf{r})=$$

$$=$$

$$H_{\text{Pairing}}=$$

$$\mathbf{c}\left(\mathbf{k}\right)=$$

$$H\left(\mathbf{k}\right)=$$

$$H(\mathbf{k})=\left(\begin{array}{cc}{H}_{s=\uparrow }(\mathbf{k})& 0\\ 0& H_{s=\downarrow }(\mathbf{k})\end{array}\right),$$

$${\gamma }_{s,n}(\mathbf{k})=$$

$$c_{\sigma }\left(x_i,\mathbf{k}\right)=\sum _{n,s}{u}_{s\sigma ,n}(x_i,\mathbf{k}){\gamma }_{s,n}(k)+v_{s\sigma ,n}^{\ast }(x_i,-\mathbf{k}){\gamma }_{s,n}^{†}(-\mathbf{k}).$$

$$H_s\left(\mathbf{k}\right){\Phi }_{s,n}(\mathbf{k})=$$

$${\Phi }_{s,n}(\mathbf{k})=$$

$$C_{\text{Singlet}}=(-i{\tau }_y)1_{N\times N}K,$$

$$C_{\text{Triplet}}={\tau }_x{1}_{N\times N}K,$$

$$C{H}_s\left(\mathbf{k}\right)C^{-1}=-H_s\left(-\mathbf{k}\right).$$

$$H_s\left(\mathbf{k}\right)\left[C{\Phi }_{s,n}(-\mathbf{k})\right]=$$

$$C{\Phi }_{s,0}(-\mathbf{k})=$$

$$=$$

$$H_s(k)=H_s(-k).$$

$$-H_s\left(-k\right)\left[C{\Phi }_{s,n}(k)\right]=$$

$$H_s\left(k\right)\left[C{\Phi }_{s,n}(k)\right]=$$

$$H=$$

$$H_{TB}=$$

$$H_{\text{Pair}}=$$

$$t(x_i)=$$

$$\delta t(x_i)=$$

$$c_{\sigma }(x_i,y_i)=$$

$$H_{TB}=$$

$$+$$

$$H_{\text{Pair}}=$$

$$t(x_i)\left[u_{s,n}(x_{i+1},k)+u_{s,n}(x_{i-1},k)\right]-\mu u_{s,n}(x_i,k)+2t(x_i)\cos \left(ka\right)u_{s,n}(x_i,k)$$

$$+$$

$$-i\Delta \sin \left(ka\right)\left[u_{s,n}(x_{i+1},k)-u_{s,n}(x_{i-1},k)\right]$$

$$-$$

$$c_{\sigma }\left(x_i,y_i\right)=$$

Full text

Strain engineering has emerged in recent years as a promising way to access topologically non-trivial phases. The initial proposal of Landau levels (LLs) in strained graphene Guinea et al. (2009) was followed by a number of extensions, most notably to nodal Dirac superconductors (SCs) Nica and Franz (2018); Massarelli et al. (2017) and Weyl semi-metals and SCs Grushin et al. (2016); Liu et al. (2017). These intriguing proposals raise questions regarding the nature and the stability of the strain-induced gapless modes. Such issues have not been addressed in detail except in the case of strained grapheneGhaemi et al. (2012). Here, we study the bulk gapless LLs which emerge under applied non-uniform strain in prototypical models which describe a large class of well-studied 2D d-wave SCs, as exemplified by high-$T_{c}$ cuprates Hashimoto et al. (2014). Furthermore, we consider possible instabilities driven by small residual interactions in the low-energy sector. We likewise analyze 3D time-reversal symmetric equal-spin triplet SCs.

We first elucidate the general spinful Majorana nature of the bulk zeroth LLs which satisfy

[TABLE]

where $s=\uparrow/\downarrow$ is a spin index with $M=-\sigma_{y}$ for 2D singlet (2D-S) SCs with $d_{xy}$ pairing. Similarly, $M=-i\sigma_{z}$ for 3D triplet (3D-T) with $p\pm ip$ pairing. Subsequently, we show that the topological properties of the strained bulk depend on the direction of the uniaxial strain. When non-uniform strain is applied along an axis of the Brillouin Zone (BZ), a single zeroth LL per spin $s$ naturally emerges as a consequence of a phase transition in the bulk which separates two topologically distinct gapped regions. Whenever two independent zeroth LLs per spin are present, as in the case of 2D-S SCs with strain along a diagonal of the BZ, the bulk is topologically trivial. Furthermore, we show that the instabilities of these strained systems depend on the Majorana properties of the zeroth LLs. In 2D-S systems, the halving of the degrees-of-freedom (DOF) associated with the zeroth LLs effectively excludes any valley-polarization instabilities, unlike the case in graphene Ghaemi et al. (2012). Instead, the zeroth LLs in flat bands are generically unstable to ferromagnetism induced by small residual short-range interactions. By contrast, we show that the strained 3D-T SCs exhibit emergent real-space Majorana gapless modes in the bulk. Generic density-density interactions are suppressed in these cases, but non-local spin-exchange four-fermion terms of the form $U(\bm{R},\bm{R^{\prime}})c_{P,\sigma}(\bm{R})c_{P,\bar{\sigma}}(\bm{R})c_{P,\sigma}(\bm{R}^{\prime})c_{P,\bar{\sigma}}(\bm{R}^{\prime})$ are allowed and can lead to more exotic instabilities. The real-space Majorana fermions $c_{P,\sigma}(\bm{R})$ arise upon projection onto the zeroth LL sector. Our conclusions are supported by detailed analytical calculations, which are available in the Supplementary Materials (SM), as well as by numerical results presented further down.

Spinful Majorana fermions - The pairing Hamiltonian under applied uniaxial strain along the $x$ direction can be written as

[TABLE]

where $\Psi^{T}(\bm{k})=\left(c_{\sigma}(x_{i},\bm{k}),(i\sigma_{y})_{\sigma\sigma^{\prime}}c^{{\dagger}}_{\sigma}(x_{i},-\bm{k})\right)$. The index $i$ runs over all $N_{x}$ lattice sites along the $x$ direction. The conserved momentum $\bm{k}$ is perpendicular to the direction of the applied strain. The normal part given by the spin-independent $\hat{h}\sim\sigma_{0}\tau_{z}$ contains the effects of strain. The pairing part is given by $\hat{\Delta}=\hat{\Delta}_{S}\sigma_{0}\tau_{x},\hat{\Delta}_{T,x}\sigma_{x}\tau_{x}+\hat{\Delta}_{T,y}\sigma_{y}\tau_{x}$, for 2D-S and 3D-T cases, respectively. Here, the $\sigma$ and $\tau$ represent Pauli matrices in spin and Nambu spaces, respectively. The Hamiltonian represents a set of effective 1D chains defined for each conserved momentum $\bm{k}$. For more details, we refer the reader to the SM.

In the absence of any strain, the low-energy solutions of the Hamiltonian determine four nodal points on a generic Fermi surface (FS) for the 2D-S and two nodal points for 3D-T cases. The resulting spectra are illustrated in Figs. 1 (a) and (b), respectively.

By analogy with graphene Guinea et al. (2009), weak and non-uniform strain plays the role of an effective vector potential in the low-energy theory and the spectrum exhibits LLs. As pointed out in Ref. Nica and Franz, 2018, the strained system is time-reversal symmetric. Consequently, the pseudo-magnetic fields are not subject to Meissner screening, in contrast to the case of genuine magnetic fields where LLs are suppressed Franz and Tes̆anović (2000).

In either 2D-S or 3D-T cases, the Bogoliubov-de Gennes (BdG) Hamiltonian can be decoupled into two independent spin sectors $H_{s}(\bm{k})$. To illustrate, in the 2D-S case the two sectors involve $c_{\uparrow}(x_{i},\bm{k}),c^{{\dagger}}_{\downarrow}(x_{i},-\bm{k})$ and their spin-flipped counterparts, respectively. In the 3D-T case, the two sectors involve fermions of equal spin. The BdG equations for each spin sector are $H_{s}(\bm{k})\Psi_{ns}(\bm{k})=E_{ns}\Psi_{ns}(\bm{k}),$ where $\Psi^{T}_{ns}=\left(u_{ns}(x_{i},\bm{k}),v_{ns}(x_{i},\bm{k})\right)$ is a BdG spinor, $s$ is the spin of the BdG quasiparticles, and $n$ is a LL index. We define particle-hole (p-h) transformations Sato and Ando (2017) for singlet and triplet pairing as $C_{S}=(i\tau_{y})\otimes(\bm{1})_{N\times N}K$ and $C_{T}=\tau_{x}\otimes(\bm{1})_{N\times N}K$, respectively, where $K$ represents complex-conjugation. Since $CH_{s}(\bm{k})C^{-1}=-H_{s}(-\bm{k})$, it follows that $H_{s}(\bm{k})$ has eigenstates $\Psi_{sn}(\bm{k})$ and $C\Psi_{sn}(-\bm{k})$ at energies $E$ and $-E$, respectively.

Whenever $H_{s}(\bm{k})$ exhibits a single zeroth LL, only one of the two states $\Psi_{s0}(\bm{k})$ and $C\Psi_{s0}(-\bm{k})$ is a non-trivial solution. In practice, only half of all momenta for each spin sector are associated with independent DOF. We take this as a definition of spinful Majorana modes. Equivalently, the analytical solutions for all cases (Sec. II C, III B, and IV B of the SM) indicate that the BdG quasiparticles in the zeroth LLs are pair-wise related via Eq. 1. In order to reconcile this redundancy with the BdG transformation, we constrain the zeroth LL sector to half of the allowed momenta. This procedure preserves the correct normalization of the $c_{\sigma}(x_{i},\bm{k})$ operators, as shown in Sec. II D of the SM.

*Topological origin of zeroth LLs - * Whenever strain induces a single zeroth LL of spin $s$ and momentum $k_{y}$ in the bulk, it is located at the interface of two fully-gapped and topologically distinct phases. To illustrate, we consider an effective 1D chain which is gapped in the pristine 2D-S case. Non-uniform strain along one of the axes locally closes the gap in the bulk of the 1D chain, and induces a transition between two topologically distinct phases which gives rise to the lowest LL. This mechanism is illustrated in Figs. 1 (c) and by our numerical results presented further down. A similar picture emerges in the 3D T case. Whenever two independent zeroth LLs per spin $s$ are induced in the bulk by applied strain, the gapped regions on either side are topologically trivial. This occurs in the 2D-S case with strain along a diagonal (Fig.1(b)).

Thus the zeroth LLs are closely related to topologically nontrivial edge states. Indeed, under applied strain, our effective 1D and 2D models belong to the DIII class of the Altland-Zirnbauer classification Teo and Kane (2010). The latter predicts Majorana Kramers doublets Teo and Kane (2010) for point defects in 1D, while line defects in 2D lead to helical Majorana fermions, both associated with a non-trivial bulk $Z_{2}$ invariant. This is precisely what we obtain in the 2D-S case with strain along an axis of the BZ and in the 3D-T case, respectively. Similarly, the presence of two zeroth LLs in the 2D-S case with strain along a diagonal of the BZ signals topologically trivial bulk phases.

We stress that the topological transition in the bulk via applied strain does not imply destruction of the pairing ground-state. We argue by analogy to vortex states in standard Type-II SCs, where the condensate likewise survives de Gennes (1999).

*Results - * We first consider the case of a 2D-S cases. Referring to Fig. 1 (a), the unstrained system exhibits Dirac points on the axes of the BZ, located at $\bm{k}=\left(0,\pm K_{F}\right)$ and labeled by valley indices $(1)$ and $(1^{\prime})$. The pair of valleys located at $\bm{k}=\left(\pm K_{F},0\right)$ are labeled by $(2)$ and $(2^{\prime})$, respectively. We introduce uniaxial strain along the $x$ direction (green arrow in Fig. 1 (a)) via a slowly-varying correction to the hopping term $\delta t(x_{i})=t\epsilon x_{i}$, where $t$ is the pristine hopping coefficient. As discussed in Sec. V, the strain parameter $\epsilon$ is proportional to the gradient of the applied strain. The detailed form of the lattice Hamiltonian is given in Sec. II A of the SM. Next we solve this system analytically in the continuum limit about each of the four valleys. Details of the calculation are presented in Sec. II B and C of the SM. We find that strain induces LLs at valleys $(1)$ and $(1^{\prime})$. Furthermore, the zeroth LL states obey Eq. 1. Therefore, the gapless modes of opposite valley and spin indices are not independent.

Likewise, the numerical solutions for momenta approaching the Dirac points are consistent with the analytical results. To illustrate, in Fig. 2 (a) we show the spectrum of $H_{s=\uparrow}$ about valley $(1)$. The results were determined using a chain of 500 sites with lattice spacing $a=1$, pairing amplitude $\Delta=0.1$, and strain parameter $\epsilon=9\times 10^{-4}$. We can distinguish the two zero-energy states via their respective spectral weights. In Fig. 2 (b) we show the spectral weight $|u(x_{i},k_{y})|^{2}+|v(x_{i},k_{y})|^{2}$ as a function of position along the effective 1D chain, for the several decreasing values of $k_{y}$ indicated in Fig. 2 (a). As we approach the center of the flat zeroth LL band, the spectral weight shows well-defined peaks in the bulk and at the edge, corresponding to a zeroth LL and an edge state, respectively. This confirms the topological origin of the zeroth LL modes as pointed out in Fig. 1 (c). Indeed, the bulk of the strained sample now contains two topologically-distinct phases, separated from each other by the zeroth LL. Likewise, an edge state emerges at the boundary between the vacuum and the non-trivial sector. As we move away from the node, the zeroth LL merges with an edge state (orange points), separating the topologically trivial vacuum from the non-trivial bulk. These results also illustrate that there is only one bulk, zeroth LL per spin sector, confirming it’s Majorana nature. A similar result is obtained for the opposite spin sector and about valley $(1^{\prime})$.

Our hypothesis can be further supported by considering the evolution of the gapped system at $k_{y}=1.84$, as a function of applied strain, as shown in Fig. 3 (a). Here $E_{0}\approx O(10^{-2})$ is the value of the gap in the absence of strain. As the strain parameter increases, the spectral weight is increasingly localized at the edge, as shown in Fig. 3 (b). Upon reaching a strain parameter $\epsilon\approx 1.2\times 10^{-4}$, the gap first decreases dramatically then stabilizes close to zero. The corresponding spectral weights in Fig. 3 (c) exhibit two sharp peaks corresponding to zeroth LL and edge states, respectively.

Next, we discuss the many-body instabilities of the spinful Majorana fermions. We consider residual short range interactions of the Hubbard type in the paired state Potter and Lee (2014) (see SM Sec. II E). Provided that the cyclotron frequency $\omega_{c}\sim\sqrt{B}$ associated with the pseudo-magnetic field $B$ is much smaller than the interaction strength U’, we project the lattice operators onto the zeroth LL sector. Since only half of the zero modes are independent (Eq.1), we can formally eliminate the valley $(1^{\prime})$ DOF (See SM Sec. II D). Consequently, we obtain an effective Hubbard interaction involving quasiparticles of both spins at valley $(1)$. At mean-field level, this model predicts an instability toward ferromagnetism with a flat gap of magnitude $n_{0}U^{\prime}/2$. Here, $n_{0}\sim l^{-2}_{B}$ is the density of zeroth LLs and $l_{B}$ is a pseudo-magnetic length. The remaining possible instabilities are either energetically unfavored, as in the case of singlet pairing, or are excluded due to the halving of the DOF in the case of a spin-density wave.

Note that LLs do not emerge at valleys $(2)$ and $(2^{\prime})$. As shown in Sec. II C of the SM, in the continuum limit, the effects of the strain can be incorporated via a global phase. The resulting states remain gapless and exhibit an approximate linear-in-$q$ dispersion. Using standard RPA arguments, we see that coupling the projected zeroth LLs to the gapless modes at valleys $(2)$ and $(2^{\prime})$ does not qualitatively change our conclusion.

We now consider 2D-S SCs with strain applied along the diagonal of the BZ (red arrow in Fig. 1 (a)). These differ from the analogous cases of strain along an axis of the BZ in several ways. First, there are only two valleys at $k\approx\pm K_{F}$, where the conserved momentum is perpendicular to the diagonal of the BZ. Secondly, the 1D effective models are defined on two sublattices. Moreover, the Hamiltonians are invariant under a reflection $H_{s}(k)=H_{s}(-k)$. Together with the p-h symmetry discussed previously, it ensures that zero-energy modes of the same spin and at the same valley emerge in pairs with BdG coefficients $\Psi^{T}_{s}(k)=\left(u_{s}(x_{i},k),v_{s}(x_{i},k)\right)$ and $C\Psi_{s}(k)=(-v_{s}(x_{i},k),u_{s}(x_{i},k))$. Note the absence of the minus sign for the momentum of the second spinor. For more details, please consult Sec. III A of the SM.

A consequence of the combined mirror and p-h symmetry is that there are two independent zeroth LLs per valley and spin, as illustrated by analytical solutions of the BdG equations in the continuum limit (Sec. III B of the SM). The numerical solutions of the effective 1D lattice model likewise support this claim. In Fig. 4, we present our numerical results for a chain of 500 sites with pairing amplitude $\Delta=10^{-1}$ and strain parameter $\epsilon=5\times 10^{-3}$. The BdG coefficients of the two degenerate zero-modes in panel (a) are precisely of the form imposed by the enhanced symmetry, as indicated in panels (b) and (c). This figure also illustrates the absence of any zero-energy edge states, thereby confirming our previous statement that the bulk is in a topologically trivial phase.

The pairs of independent zeroth LLs at opposite valleys are related via an analog of Eq. 1. In turn, each pair of independent zeroth LLs determine the projections of operators defined for the two sublattices. Upon inclusion of a repulsive Hubbard interaction, as in the case of strain applied along an axis of the BZ, we find that the leading instability is also toward ferromagnetism. A more detailed discussion is given in Sec. III C and D of the SM. The ferromagnetic state shows a full gap in the bulk. This state would exhibit exponential activation in temperature dependence of the specific heat or the penetration depth to name a few.

We discuss the case of 3D-T case next. In the absence of strain, we choose a $\Gamma^{-}_{1}$ representation of the tetragonal $D_{4h}$ group Sigrist and Ueda (1991) for the pairing. This leads to two nodal points at $k_{z}=\pm K_{F}$ for a spherical FS, as illustrated in Fig. 1 (b). This type of pairing can be thought of time reversal symmetric analog of the pairing in He-3A Leggett (1975); Volovik (2003). As discussed in Sec. IV A and B of the SM, under uniaxial non-uniform strain along the $x$-direction, and for fixed $k_{z}\approx\pm K_{z}$, the low-energy spectrum of $H_{s=\uparrow}$ along $k_{y}$ consists of single branches of right-moving bulk modes at each valley $(1)$ and $(1^{\prime})$. For fixed $k_{y}\approx 0$, a flat dispersion emerges along $k_{z}$. Likewise, the opposite-spin sector $H_{s=\downarrow}$ exhibits two left-moving modes at each of the two valleys. Moreover, we find that for each spin sector, the modes at opposite valleys are related via Eq. 1. These solutions are consistent with the general discussion of Ref. Teo and Kane, 2010, which predicts helical Majorana fermions associated with a line defect in a gapped 2D system belonging to the DIII class.

The analytical results are confirmed by our lattice calculations. In Fig. 5 (a), we show the spin up spectrum as a function of $k_{y}$ at valley $(1)$ corresponding to $k_{z}\approx K_{F}$. These results were obtained for a chain of 500 sites with a pairing amplitude of $\Delta=10^{-1}$ and strain parameter $\epsilon=9\times 10^{-3}$. As illustrated by Fig. 5 (b), only one of the two branches corresponds to bulk states, while the other is localized at the edge, confirming the scenario set out in Fig. 1 (c). Furthermore, we see that the two branches are interchanged for the spin-down sector at the same valley, as indicated by Figs. 5 (c) and (d).

The Fermi fields in the unpaired basis, projected onto the zeroth LL for each spin sector obey a real-space Majorana condition (see SM Sec. IV C)

[TABLE]

since the projected operators at each valley are related via $c^{(1^{\prime})}_{P,\sigma}=\sum_{\sigma^{\prime}}(i\sigma_{z})_{\sigma\sigma^{\prime}}c^{(1),{\dagger}}_{P,\sigma^{\prime}}$. Note that a similar relation does not hold in the singlet case, where zero-modes of opposite spin are related to each other.

In contrast to the 2D-S case, the real-space Majorana nature of the projected fields in the 3D-T case imposes strong constraints on possible residual interactions. It is well-known Elliott and Franz (2015) that the product of identical real-space Majorana operators reduces to a constant. This reflects the fact that ground-states of Majorana systems do not conserve particle number Elliott and Franz (2015). It follows that any local density $c^{{\dagger}}_{\sigma}(\bm{R}_{i})c_{\sigma}(\bm{R}_{i})\rightarrow c_{P,\sigma}(\bm{R}_{i})c_{P,\sigma}(\bm{R}_{i})=\text{Constant}$. As such, density-density interactions of any range between projected Majorana operators are dynamically trivial. However non-local, four-fermion exchange terms of the form $U(\bm{R},\bm{R^{\prime}})c_{P,\sigma}(\bm{R})c_{P,\bar{\sigma}}(\bm{R})c_{P,\sigma}(\bm{R}^{\prime})c_{P,\bar{\sigma}}(\bm{R}^{\prime})$ are allowed. These are unlike spinless Majorana fermions where a minimal interaction is defined on a plaquette Rahmani et al. (2015); Affleck et al. (2017); Wamer and Affleck (2018); Rahmani et al. (2019). To our knowledge, interactions between spinful Majorana fermions have not been proposed in a condensed matter setting. We anticipate that they can lead to more exotic phases which require further investigation.

*Discussion - * We showed that non-uniform uniaxial strain leads to the formation of spinful Majorana fermions in SCs with a Dirac nodal spectrum. We also demonstrated that the applied strain can induce topological transitions in the bulk of these gapless SCs. The resulting zeroth LLs are localized at the boundary between two topologically distinct phases within the strained sample. In many ways, these bulk states are analogous to the better-known cases of Majorana fermions localized in vortex cores of topological SCs Fu and Kane (2008). We note that a similar correspondance between edge states and the zeroth LLs in a strained Weyl semi-metal was pointed out in Ref. Grushin et al., 2016. To our knowledge, the more general mechanism involving a local closing of the gap has not been clearly identified as such.

We estimate a pseudo-magnetic length $l_{B}\approx 410~{}a$, where $a$ is the inter-atomic distance. We obtain a separation in energy of the LLs $E_{c}\approx 0.11~{}$meV, corresponding to a characteristic temperature scale of 1.2 K. These estimates were determined using parameters typical of the cuprate class of SCs and are discussed in detail in Sec. V of the SM.

As we have shown, the direct observation of these gapless states is precluded by various instabilities. In strained 2D d-wave SCs, the flat zeroth-LL band is unstable toward ferromagnetism. One immediate consequence is the coexistence of superconductivity and ferromagnetism in the bulk of a strained d-wave SC. The finite magnetic moment in the bulk is screened via the Meissner effect. However, signatures of the resulting screening currents could be detected by SQUID Frolov et al. (2008) or NMR experiments. For 3D triplet pairing, the low-energy gapless states are emergent real-space spinful Majorana fermions, which are likely to host more exotic phases.

Acknowledgements We thank Joel Moore for fruitful discussions. This work is supported by ASU startup grant.

I General Aspects of the Bogoliubov-de Gennes Hamiltonian

In this section, we discuss the general properties of spinful Majorana fermions which occur for time-reversal invariant singlet and equal spin-triplet pairing Hamiltonians under applied strain. The general discussion in this section is exemplified in subsequent sections which treat specific cases in detail.

We define a lattice Hamiltonian in $D=2,3$ dimensions under applier strain as

[TABLE]

where spinors are given by

[TABLE]

The position $x_{i},i\in\{1,\ldots,N\}$ indicates the position along the direction of the applied uniaxial strain, while $\bm{k}$ represents the $D-1$ component conserved momentum along the directions which are perpendicular to the applied strain. The Hamiltonian can be written as

[TABLE]

Due to the presence of the strain, each block is a matrix which depends on both position indices, as well as the momentum.

In either singlet or equal-spin triplet pairing cases, we can separate Hamiltonian into two independent spin sectors as

[TABLE]

corresponding to spin $(\uparrow\downarrow)$, $(\downarrow\uparrow)$ sectors in the singlet case, and spin $(\uparrow\uparrow)$ and $(\downarrow\downarrow)$ sectors in the equal-spin triplet cases.

The solutions of each $H_{s}$ are then obtained via a standard Bogoliubov-de Gennes (BdG) ansatz Elliott and Franz (2015)

[TABLE]

for energies $E\geq 0$. We likewise define the inverse transformation as

[TABLE]

The index $s$ also labels the conserved spin of the BdG quasiparticles. For the case of singlet pairing $u_{s\sigma}=\delta_{s\sigma}u_{s}$ and $v_{s\sigma}=(i\sigma_{y})_{s\sigma}v_{s}$. In the case of equal-spin triplet pairing $u$ has an identical form while $v_{s\sigma}=-(\sigma_{y})_{s\sigma}v_{s}$. The index $n$ represents all other labels. In the vicinity of the nodal points, it labels the bulk Landau Levels (LLs).

As discussed below and in the following sections, we find that only half of the degrees-of-freedom for the zeroth LL sector are independent since $\gamma_{s,0}(\bm{k})=M_{ss^{\prime}}\gamma^{{\dagger}}_{s^{\prime},0}(-\bm{k})$, where $M=\sigma_{y}$ for 2D-S, and $M=-i\sigma_{z}$ for singlet and triplet cases, respectively. Note that the inverse BdG transformation in Eq. S9 implicitly assumes that the two components are independent of each other. Therefore, in order to preserve a well-defined BdG transformation, we restrict the $u,v$ indices to half of the allowed momenta. This procedure is illustrated in Sec. II.4.

The Majorana nature of the zeroth LLs follows from the general p-h symmetry of the Hamiltonian. Indeed, each $H_{s}$ sector obeys the BdG eqs.

[TABLE]

where the BdG spinor is

[TABLE]

We define a particle-hole (p-h) transformation

[TABLE]

for singlet and

[TABLE]

for triplet pairing, respectively. The $N\times N$ identity matrix is defined for the position indices and $K$ represents complex conjugation. The two different forms are due to the distinct parities of the two types of pairing. The p-h transformation acts on the Hamiltonians for each of the two spin sectors as Sato and Ando (2017)

[TABLE]

As is well-know, this implies that the eigenstates of $H_{s}(\bm{k})$ come in pairs of opposite energy:

[TABLE]

When strain induces a single, zeroth LL mode, p-h symmetry imposes an additional constraint. Since only one of either $\Phi_{s,0}(\bm{k})$ or $C\Phi_{s,0}(-\bm{k})$ is a non-trivial solution of $H_{\sigma\sigma^{\prime}}(\bm{k})$, we require

[TABLE]

This necessarily implies that the Hilbert space of the zeroth LLs must be truncated to half of all allowed momenta $\bm{k}$.

The above argument excludes the case where the spinor is invariant under the p-h transformation: $C\Psi_{s}(-\bm{k})\propto\Psi_{s}(\bm{k})$. This is confirmed by explicit solutions presented in Sec. II.3.

Finally, we discuss the case where zero-energy LLs occur in pairs, as discussed in detail in Sec. III.2 and III.3. Consider the case where $H_{s}$ are invariant under a mirror plane involving the conserved momentum $k$:

[TABLE]

Together with the p-h transformation

[TABLE]

this implies

[TABLE]

Therefore, zero-energy solutions occur in pairs and are determined by $\Psi_{s,0}$ and $C\Psi_{s,0}(k)$.

II Two dimensional d-wave pairing with strain along one nodal axis

In this section, we discuss the case of a 2D $d_{xy}$ SC with strain along one of the nodal axes of the pairing. In Sec. II.1, we discuss the effective 1D lattice Hamiltonians which are appropriate under applied uniaxial strain. An effective low-energy Hamiltonian is derived in Sec. II.2 and it’s detailed LL solutions are presented in Sec. II.3. Sec. II.4 discuss the resulting projection of the lattice operators onto the zeroth LL sector. The last subsection is devoted to the discussion of a short-range Hubbard interaction projected onto the zeroth LLs and the ensuing ferromagnetic instability.

II.1 1D Lattice Hamiltonian in the presence of uniaxial strain

Using the convention of Eq. S5, we consider a simple lattice model for a 2D d-wave consisting of nearest-neighbor (NN) hopping together with next-nearest-neighbor (NNN) $d_{xy}$ pairing on a square lattice:

[TABLE]

Note that the pairing amplitude is independent of spin indices $\sigma,\sigma^{\prime}$. We define the NN and NNN vectors $\bm{\delta}_{1}=(a,0),\bm{\delta}_{2}=(0,a)$ and $\bm{\delta^{\prime}}_{1}=(a,a),\bm{\delta^{\prime}}_{2}=(-a,a)$, where $a$ is the NN distance. The sign change of the $d_{xy}$ pairing is induced via $\Delta(\bm{\delta^{\prime}}_{1/2})=\pm\Delta$. Without any strain, we have $t(x_{i})=t$ and the resulting dispersion shows gapless excitations around the four momenta $(\pm K_{F},0)$ and $(0,\pm K_{F})$, which denote the positions of four valleys denoted by $(2),(2^{\prime})$ and $(1),(1^{\prime})$, respectively.

We allow for the effects of strain by introducing a position dependent hopping along the $x$ direction

[TABLE]

where $\epsilon$ is approximately proportional to the gradient of the strain along $x$ (Sec. V.1). The modification to the unperturbed hopping $t$ is due to the non-uniform stretching of the NN bonds. For a microscopic derivation of this term, we refer the reader to the supplementary material of Ref. Nica and Franz, 2018. We can in principle allow an analogous variation in the pairing amplitudes. However, as shown in the supplementary material of Ref. Nica and Franz, 2018, in the low-energy limit, such a contribution can be eliminated via a gauge transformation and does not change our results qualitatively.

Our central assumption is that the contribution to the hopping due to strain is small compared to any of the parameters of the unperturbed system over the finite length of the sample i.e. $\text{max}(\delta t(x_{i}))\ll\Delta,t$. As a consequence, we assume that the applied strain does not destroy the d-wave superconductor.

We assume periodic boundary conditions along the $y$ direction and apply a Fourier transform:

[TABLE]

where the 1D Brillouin Zone (BZ) is defined by $k_{y}\in[-\pi/a,\pi/a]$. The Hamiltonian becomes

[TABLE]

Using the ansatz for singlet pairing in Eq. S8, we obtain the BdG Eqs for each $H_{s}$ sector as

[TABLE]

where we dropped the $y$ index.

II.2 Low-energy effective Hamiltonian under applied uniaxial strain

We consider the low-energy limit of the lattice BdG eqs. As we are interested in solutions in the bulk of a sample, we consider an infinite system and ignore the effects of the boundary along the direction of the strain. Our numerical results justify this approximation.

Note that the lattice and continuum fields are related via the Wannier states Auerbach (1994)

[TABLE]

A partial Fourier transform along the $y$ direction gives

[TABLE]

For sufficiently large systems, we allow the Fourier transform of the Wannier function to be sharply peaked

[TABLE]

Using these expression in Eq. S8, we obtain

[TABLE]

where we defined the continuum version of the BdG coefficients as

[TABLE]

We assume for simplicity that the Wannier states can be chosen to be real. We can formally switch to the continuum by multiplying the each one of eqs. S29 by the Wannier state $\tilde{\phi}$ and summing over all $i$. We subsequently drop the index $i$.

In the low-energy limit and for momenta in the vicinity of each valley s.t. $k\approx K^{(\alpha)}_{Fy}+q_{y}$, we employ the approximation

[TABLE]

where $\alpha\in\{1,1^{\prime},2,2^{\prime}\}$ is a valley index and $(u,v)^{(\alpha)}$ are envelope functions which are assumed to vary slowly on the scale of the lattice. Such an approximation is well-known in the context of graphene Castro et al. (2009). We introduce a cutoff scale $|q|\leq\Lambda$.

Our ansatz allows for a separation of the solutions in terms of valley index. We simplify common phase factors $e^{iK_{x}x}$ and expand the envelope functions s.t.

[TABLE]

We can further expand

[TABLE]

We keep only terms in Eqs. S29 which are first order in either of the three small quantities $\delta t(x),\partial_{x},q$. We note that the zeroth order terms vanish either due to presence of a Fermi surface or to the vanishing of the pairing at each valley. To first order, the first Eq. of S29 becomes

[TABLE]

We define

[TABLE]

where $K_{F}=\text{max}(|K^{(\alpha)}_{Fx}|,|K^{(\alpha)}_{Fy}|)$ is independent of $\alpha$ and positive. The BdG Eqs can be recast as

[TABLE]

II.3 Landau levels from the low-energy Hamiltonian

We first focus on the (1), (1’) valleys which correspond to $\bm{K_{F}}=(0,\pm K_{F})$. The strain-induced vector potential induces LLs about the nodal momenta. By contrast, LLs do not emerge in the vicinity of valleys (2), (2’), as shown further down below.

The BdG Eqs. in the vicinity of valley (1) are given by

[TABLE]

We apply a unitary transformation

[TABLE]

where $\sigma_{x}$ acts on the space of the BdG spinor. The new Bdg Eqs. for valley (1) read

[TABLE]

where we introduced the effective vector potential

[TABLE]

Recalling that $\tilde{A}(x)=Bx$ by construction, we define a pseudo-magnetic field as

[TABLE]

We now clearly see that the BdG equations correspond to Dirac fermions in the presence of a uniform magnetic field in the Landau gauge Castro et al. (2009). The well-known solutions for $n\neq 0$ are given by the Landau levels

[TABLE]

with energies

[TABLE]

The cyclotron frequencies and magnetic lengths are defined by

[TABLE]

$\Psi_{n}(x,q)$ represent harmonic oscillator wavefunctions Messiah (1999)

[TABLE]

where $H_{n}$ are Hermite polynomials and the dimensionless variable is $\xi=(x/l_{B})-(v_{F}/v_{\Delta})l_{B}q$. The $\psi_{n}$ functions obey the normalization condition

[TABLE]

Referring to the BdG Eqs. (S47), for valley (1’) we have

[TABLE]

These can be obtained from Eqs. S48 for valley (1) via complex conjugation and sending $q\rightarrow-q$. This procedure amounts to a time-reversal operation and illustrates the that this symmetry is preserved by the applied strain.

In the un-rotated basis, the solutions for the non-zeroth LLs are

[TABLE]

where $\text{sgn}(\alpha)=\pm 1$ for $\alpha\in\{1,1^{\prime}\}$.

The zeroth LL in the original basis is determined by

[TABLE]

The continuum version of the BdG ansatz in Eq. S8 is given by

[TABLE]

In the spin-singlet case $u_{s\sigma}=\delta_{s\sigma}$ and similarly for $v_{s\sigma}=(i\sigma_{y})_{s\sigma}v_{s}$. The explicit expressions are

[TABLE]

It follows that

[TABLE]

We now comment on the solutions at valleys (2), (2’). Based on Eqs. S47, the first of these is subject to the following eqs.

[TABLE]

These do not result in the emergence of LLs. The vector potential can be taken into account by incorporating a suitable global phase:

[TABLE]

The solutions correspond to linearly-dispersing, counter-propagating modes centered around $K_{Fy}=0$.

We finally comment on the topological nature of the zeroth LLs. To illustrate, we consider the zero-energy solutions of Eqs. S50. By canceling factors of $i$ these eqs. can be reduced to

[TABLE]

where $m(x)=v_{F}q-A(x)$. These are identical to the well-known Jackiw-Rebbi eqs. Jackiw and Rebbi (1976). The latter host a topologically-protected domain-wall solution centered about the point where $m(x)$ changes sign. Going beyond the low-energy approximation, the zeroth LL modes are expected to persist as topologically-protected modes centered around line-defects, in accordance with the findings of Ref. Teo and Kane, 2010.

II.4 Projection of lattice operators onto the zeroth LL modes

We consider solutions deep in the bulk of a strained sample and thus ignore contributions from any edge states. As before, we can pass onto the continuum via an expansion using the Wannier states. The continuum version of the inverse BdG transformation is

[TABLE]

Note that for the zeroth LL contribution, only two out of the four states in Eqs. S66-S69 are independent. This necessarily implies that the inverse BdG transformation must be restricted to half of the the zeroth LL states. Without loss of generality, we write the expansion of the fields as

[TABLE]

where we introduced the projection onto the zeroth LL

[TABLE]

We chose to keep both positive and negative values of $q$, but retained only those contributions from valley $(1)$. These expressions are consistent with Eqs. S66-S69. Using the completeness of the Hermite polynomials (Eq. S59) it is straightforward to verify that the fields in Eq. S77 satisfy canonical anti-commutation relations.

Via Eqs. S70-S71, and S76, one can verify that, without any restriction on the Hilbert space of the zeroth LLs, the inverse BdG transformation is indeed singular. The projections onto the zeroth LL are themselves not independent since

[TABLE]

We now illustrate that Eq. S82 is independent of the choice of Hilbert space for the zeroth LL. Consider retaining both valley indices, but restricting the momenta $q$ to positive values. The contribution of the zeroth LLs is given by

[TABLE]

Using Eqs. S70-S71 one can verify that these satisfy the relations in Eqs. S82.

The Fourier transforms of the projected operators, defined by

[TABLE]

likewise obey

[TABLE]

In addition, they are subject to the anti-commutation relations

[TABLE]

[TABLE]

where $\lambda=v_{F}/v_{\Delta}$. In the fourth line we approximated the integral by a Gaussian in the limit $\Lambda l_{B}\rightarrow\infty$. This approximation is justified for weak pseudo-magnetic fields considered here. We identify the coefficient as half of the density of LLs

[TABLE]

by using the well-known formula for the degeneracy of LLs

[TABLE]

where $L_{x},L_{y}$ are the lengths of the system in either direction and we accounted for the anisotropy between the Fermi velocities via the factor $\lambda$.

The field operator projected onto the zeroth LLs gives

[TABLE]

II.5 Ferromagnetism at mean-field level

In order to investigate many body instabilities of the zeroth LL, we consider a repulsive, short-range Hubbard interaction in the continuum limit

[TABLE]

Using Eq. S93, we project this term onto the zeroth LLs:

[TABLE]

Note that we have only kept terms which lack rapidly oscillating phase factors such as $e^{\pm 2iK_{F}x}$. We anticipate that in the long-wavelength limit considered here, the terms we kept are the most relevant. We have also ignored any contributions from boundary terms.

Using Eq. S88, we can formally eliminate the degrees of freedom at valley $(1^{\prime})$ in favor of their analogoues at valley $(1)$:

[TABLE]

We take advantage of the fact that the anti-commutators of the projected operators is finite and proportional to half of the density of zeroth LLs per spin (Eq. S89). This way the Hamiltonian reduces to

[TABLE]

We write

[TABLE]

where

[TABLE]

The Hamiltonian can be re-written as

[TABLE]

The Hamiltonian can be decoupled at mean-field level in the ferromagnetic channel. For repulsive interactions, the pairing instability is energetically unfavorable. A possible competing spin-density order is not allowed, since we have eliminated the degrees of freedom at valley $(1^{\prime})$.

Taking the local spin-density and total density as variational parameters, we take the expectation value of $H_{U,P}$. From the expression above, it is clear that a half-filled ground state corresponding to $\rho^{(1)}_{P}=n_{0}/2$ is energetically favored.

Using S87, the spin-part is given by

[TABLE]

The Gaussian term ensures that the dominant contributions are due to momenta $k,q,Q$ which are within all within a scale set by $1/l_{B}$. At half-filling, it is clear that this term is maximally negative at ferromagnetic alignment.

In conclusion, a ferromagnetic ground-state is favored at mean-field level. The resulting Stoner spectrum consists of flat bands at

[TABLE]

III Two dimensional d-wave pairing with strain off the nodal axes

In this section we discuss a 2D d-wave SC with strain applied off the nodal axes of the pairing. In Sec. III.1, we introduce the effective 1D lattice Hamiltonians under uniaxial strain. We also discuss their enhanced symmetry and sublattice structure, which are important differences w.r.t. the models of Sec. II.1. In Sec. III.2, we introduce the low-energy, continuum limit of the BdG eqs. and obtain their solutions. We find two degenerate zeroth LLs in the vicinity of each of the two valleys with BdG coefficients which are subject to the constraint imposed by the combined mirror and p-h symmetry discussed in the main text and below. Sec. III.3 presents a comparison of the analytical solutions with the numerical results. Sec. III.4 is devoted to a discussion of the projection of the lattice operators onto the zeroth LL section, while Sec. III.5 illustrates the ensuing ferromagnetic instability.

III.1 1D Lattice Hamiltonian

We consider the tight-binding Hamiltonian in Eq. S23 without strain and rotate by $\pi/4$. In order to preserve translational symmetry along the rotated $y$ direction at the edges we choose a two-site unit cell. This lattice is illustrated in Fig. S6.

We subsequently introduce a spatially-varying tight-binding parameter $t(x_{i})$, where $x$ is along the diagonal of the un-rotated system. Applying a Fourier transform along the $y$ direction we obtain

[TABLE]

The position index runs over all even values while the momentum is along a folded 1D BZ. Without strain, the effective chain model is periodic by translation via $a\sqrt{2}$ or $i\rightarrow i+1$. Subsequently, we distinguish between two inequivalent sublattices which we label by $A,B$.

We consider the ansatz in Eq. S8 modified to take into account the sublattice structure:

[TABLE]

Using the convention in Eq. S4 where $u_{s\sigma}=\delta_{s\sigma}u_{s}$ and $v_{s\sigma}=(i\sigma_{y})_{s\sigma}v_{s}$, the singlet pairing Hamiltonian becomes independent of the spin indices $\sigma$, which we subsequently drop for simplicity.

The BdG eqs. for sublattice B are given by

[TABLE]

The corresponding eqs. for sublattice $A$ can be obtained by simply exchanging $A,B$ indices.

Finally, we note that the BdG equations for each sublattice are invariant under $k\rightarrow-k$, as discussed in the main text.

III.2 Landau levels in the continuum limit

In the limit of zero strain, the bonding solution of Eqs. S113-S114 exhibits nodes at four momenta $(\pm K_{F},\pm K_{F})$, which we label by valley indices $(1),(2),(1^{\prime}),(2^{\prime})$ in a clockwise direction. The anti-bonding solution is generically gapped and will be ignored in the following.

We extend the ansatz for the pristine case to the strained system by writing the wavefunctions as

[TABLE]

in the vicinity of $k_{y}\approx K_{Fy}+q_{y}$ with $q_{y}$ small. A similar expansion holds for $v_{A/B}$. Note that we are considering bonding solutions where the envelope functions $u^{(\alpha)}_{A}=u^{(\alpha)}_{B}$, while keeping track of the phase difference between the two sublattices. Explicitly, for positive $K_{Fy}$, $\alpha=1$ and $\bar{\alpha}=2^{\prime}$, and $K^{(\alpha)}=K_{Fx}$ while $K^{(\bar{\alpha)}}=-K_{Fx}$. A similar expansion holds for valleys (2), and (1’) with opposite signs for both $K_{Fx}$ and $K_{Fy}$. This form for the ansatz is clearly justified for zero strain, where the pairs of valleys $(1)$ and $(2^{\prime})$ and $(1^{\prime})$ and $(2)$ are decoupled. Any amount of small strain breaks translational symmetry along the $x$ direction and can in principle mix the valleys in pairs. The ansatz continues to hold, as indicated by numerical results which do not assume this decomposition (see main text).

Following the procedure in Sec. II.2, we formally pass to the continuum by expanding in terms of Wannier states. We also drop the $i$, $y$ and the sublattice index. Our ansatz becomes

[TABLE]

We introduced new field operators $\Psi^{(I/II),{\dagger}}_{\sigma}\left(x,q\right)$ in the vicinity of $k_{y}\approx\pm K_{Fy}+q$. For valleys $1^{\prime},2$ we must interchange $I,II$ indices. This decomposition is equivalent to

[TABLE]

With this ansatz, the BdG equations can be separated into valley sectors.

To illustrate the general solution, we consider valley (1) in particular. We expand the envelope functions in real space and the trigonometric functions in momentum space about the nodal points. We furthermore approximate $t(x_{i}+a/\sqrt{2})$ by $t(x_{i})$ since the effect of the strain is assumed to vary slowly on the scale of the lattice. We keep only terms which are first order in the small terms $\delta t(x),\partial_{x},q$. To zeroth order we obtain:

[TABLE]

To first order, we obtain

[TABLE]

We define

[TABLE]

where $K_{F}=|K_{Fx}|=|K_{Fy}|\geq 0$.

Thd BdG eqs. simplify to

[TABLE]

We now focus on the zero-energy solutions. We transform $\tau_{z}\rightarrow\tau_{y}$ and obtain

[TABLE]

These equations are very similar to those for the zeroth LL obtained in Sec. II. By introducing the quantity

[TABLE]

we can write

[TABLE]

The solutions are given by

[TABLE]

We recall that $A(x)\sim x$ by construction. Therefore, only one solution is normalizable since it contains a decaying Gaussian form. This is typical of the zeroth LL wavefunctions for Dirac systems in a Landau gauge.

Without loss of generality, we consider the case where $\tilde{u}$ is the allowed non-trivial solution:

[TABLE]

The equations for valley (2’) can be obtained from Eqs. S126 by changing the sign of $\partial_{x}$. The analogues of Eqs. S127-S128 involve complex conjugation and interchanging $\tilde{u}$ and $\tilde{v}$:

[TABLE]

We introduced a global phase term, which must be chosen to ensure consistency with other symmetries. We can readily check that these are consistent with the mirror symmetry given by Eq. S20. Explicitly, we have

[TABLE]

where we fixed the global phase.

The BdG eqs. at valleys (1’) can be obtained from Eqs. S126 by changing the signs of both $\partial_{x}$ and $q$. The solutions are given by complex-conjugating, interchanging $\tilde{u}$ and $\tilde{v}$ and changing the sign of $q$. These steps are nothing but the p-h transformation in Eq. S12. The explicit solutions read

[TABLE]

and

[TABLE]

Note that we allowed for additional global phases $i$ to ensure time-reversal symmetric solutions at opposite valleys. This does not violate the p-h correspondence since solutions are determined modulo a global phase.

To sum up the discussion thus far, we obtain two independent solutions for $k_{y}\approx K_{Fy}+q_{y}$, where $K_{Fy}$ is positive:

[TABLE]

These satisfy

[TABLE]

These imply that at each valley there are two independent solutions. However, solutions at opposite valleys are not independent, in analogy with the results of Sec. II.3.

III.3 Comparison to numerical solution of the lattice model

In order to compare with the results of the numerical calculation presented in Fig. 4 of the main text, we consider the linear combinations of the analytical solutions found in the previous section

[TABLE]

which imply

[TABLE]

where

[TABLE]

Similar linear combinations can be taken wrt valley (II) operators. These can be obtained via a time-reversal operation

[TABLE]

which gives

[TABLE]

where $u^{(\alpha)}_{s\sigma}=\delta_{s\sigma}u^{(\alpha)}_{s}$ and $v^{(\alpha)}_{s\sigma}=(i\sigma_{y})_{s\sigma}v^{(\alpha)}_{s}$. Hence, the BdG coefficients can be made real.

One can see that these obey the constraints imposed by the mirror symmetry discussed in the main text. Furthermore, these analytical solutions are consistent with the numerical results of Fig. 4 of the main text.

III.4 Projection onto zeroth LLs

The inverse BdG eqs. can be obtained via generalizing Eqs. S76 to the case of two suballtices:

[TABLE]

As discussed in Sec. II.4, we truncate the Hilbert space by considering only $(1)$ and $(2^{\prime})$ solutions. This ensures that the BdG transformation on the zeroth LL sector is well-defined are well-defined. The projected operators are given by

[TABLE]

with similar expressions for $B$. Using Eqs. S134, S136, it is straightforward to verify that

[TABLE]

Similar relations hold for sublattice $B$. These expressions are the analogues of Eq. S82.

Following the procedure set out in Sec. II.4, we consider the real space fields on each sublattice:

[TABLE]

Finally, the real-space fields are determined by

[TABLE]

Note that the overall phases for momenta along the $y$ direction.

III.5 Ferromagnetism at mean-field level

As for the case of strain along a nodal axis, we allow small but finite residual Hubbard interactions. In contrast to the previous case, there are now two independent states per spin at each valley, related to each other by Eqs. S143, S144. As discussed in Sec. III.1 the real-space system is described by a two-site unit cell. The two independent zero-energy LLs per spin at each valley can be used to determine the projections of the two fields $\Psi_{P,A\sigma}$ and $\Psi_{P,B,\sigma}$ defined on each sublattice, as shown in the previous section.

Consider a continuum version of the Hubbard interaction on each sublattice:

[TABLE]

Using Eqs. S169, each sublattice sector maps onto the single-valley case studied in Sec. II.5. Hence, each sublattice is expected to order ferromagnetically at mean-field level.

The large degeneracy due to possible relative orientations of the resulting sublattice moments can be lifted by arbitrarily small inter-sublattice interactions. In view of expected short range, repulsive nature of the residual interactions in a SC, and in the absence of negligible dispersion, we conclude that the most likely instability is still ferromagnetic.

IV Helical Majoranas for 3D equal-spin triplet SCs

In this Section we consider the emergence of helical Majorana modes in 3D equal-spin triplet superconductors under uniaxial strain. Sec. IV.1 introduces a family of effective 1D lattice models. In Sec. IV.2, we obtain helical Majorana states in the continuum limit. In Sec. IV.3 we show that the Fermi fields projected onto the zeroth LLs are emergent real-space Majorana fermions.

IV.1 Lattice Hamiltonian

We consider a 3D lattice model with NN hopping along either $x,y,z$ directions and strain along the $x$ directions. The pairing is chosen to belong to the $\Gamma^{-}_{1}$ representation of the tetragonal $D_{4h}$ point group Sigrist and Ueda (1991). Without strain, this can be written as

[TABLE]

where

[TABLE]

The spectrum is given by

[TABLE]

with Dirac point nodes at $\bm{K}_{1/1^{\prime}}=(0,0,\pm K_{F})$.

The lattice Hamiltonian can be chosen as

[TABLE]

where

[TABLE]

where $\bm{\delta}_{j=x,y,z}$ equals $(a,0,0),(0,a,0),(0,0,a)$, respectively. The pairing part is given by

[TABLE]

Note the difference in sign due to the convention in Eq. S174.

We apply Fourier transforms along the $y$ and $z$ directions. As in the pristine case, the Hamiltonian, and consequently the BdG equations separate into two independent sectors:

[TABLE]

In contrast to the singlet cases of Sec. II and III, the BdG eqs. for each $H_{s}$ retain a dependence on spin index $s$.

IV.2 Landau levels in the continuum limit

We proceed as in Sec. II.2 and consider ansatze

[TABLE]

for the wavefunctions in the vicinity of the two valleys located at $\bm{K}=(0,0\pm K_{F})$ and labeled by $(1)$ and $(1^{\prime})$. A similar expansion is expected to hold for the $v$’s. With this ansatz, the BdG eqs. for each $H_{s}$ further separate into two valley sectors. We proceed by taking the continuum limit. Since the envelope functions in Eq. S184 are assumed to vary slowly on the scale of the lattice, we can expand

[TABLE]

Similarly, we have

[TABLE]

We can likewise expand to first order in $q_{y},q_{z}$

[TABLE]

since $K_{Fy}=0$.

Allowing for $t(x)=t+\delta t(x)$, we write Eq. S182 by retaining terms to first order in either $q_{y},\partial_{x},\delta t(x)$:

[TABLE]

with a similar expression for Eq. S183. We define

[TABLE]

The BdG eqs. can be written as

[TABLE]

or

[TABLE]

Let us consider the solution at valley $(1)$. We assume that $\delta t(x)$ is chosen to generate a positive pseudo-magnetic field. As in Sec. II.3, we apply the transformation in Eq. S49 s.t. $\sigma_{z}\rightarrow\sigma_{y}$. The BdG Hamiltonian becomes

[TABLE]

Following Ref. Balatskii et al., 1986, we can choose the positive-energy solutions as

[TABLE]

where $\psi$’s are the normalized harmonic oscillator wavefunctions and the coefficients are given by

[TABLE]

The overall phase $\text{sgn}(s)$ is added to ensure that the solutions obey time-reversal symmetry. In the un-rotated basis the solutions are

[TABLE]

For the zeroth LL, we have

[TABLE]

and we require $\beta_{0}=1$. The solutions are given by

[TABLE]

Therefore, the chiral spin-up zeroth LL mode is right-moving while it’s spin-down counterpart is left-moving.

We apply a similar procedure for valley $(1^{\prime})$. Namely, we apply the inverse transformation which maps $\sigma_{z}\rightarrow-\sigma_{y}$ to get

[TABLE]

Note the opposite sign for the vector potential wrt valley $(1)$, in accordance with time-reversal symmetry. From Ref. Balatskii et al., 1986, we write the solutions as

[TABLE]

with

[TABLE]

In the un-rotated basis the solutions are

[TABLE]

Formally, the solutions can be converted to the initial sign of the vector potential by sending $q_{z}\rightarrow-q_{z}$ and taking into account the parity of the Hermite polynomials.

For the zeroth LL, $\beta_{0}=0$ and $\alpha_{0}=1$ and the coefficients are given by

[TABLE]

with energies

[TABLE]

The zeroth LLs at valley $(1^{\prime})$ have the same chiralities as their valley $(1)$ correspondents.

The explicit zeroth LL solutions are

[TABLE]

As for the case of singlet pairing, only two out of the above four states are independent since

[TABLE]

IV.3 Projection of fields onto the zeroth LLs

We follow the arguments of Sec. II.4 to obtain the Fermi fields via the inverse BdG transformation:

[TABLE]

In the spin-triplet case $u^{(\alpha)}_{s\sigma}=\delta_{s\sigma}u_{s}$ while $v^{(\alpha)}_{s\sigma}=(-\sigma_{z})_{s\sigma}v_{\bar{s}}$. Only two out of the four degrees-of-freedom in Eqs. S213-S216 are independent. Therefore we must truncate the Hilbert space of the zeroth LL to ensure a well-defined BdG transformation:

[TABLE]

In turn, the projections onto the zeroth LL are not independent:

[TABLE]

As in the case of singlet pairing, these relations are matched by the real-space operators

[TABLE]

We can map the valley $(1^{\prime})$ operators to those at valley $(1)$. From Sec. II.4, these satisfy

[TABLE]

[TABLE]

The projection of the field $\Psi_{\sigma}(\bm{r})$ with contributions from both valleys, onto the zeroth LL is given by

[TABLE]

We deduce that

[TABLE]

This relation clearly signals emergent real-space Majorana fermions. Indeed, these projected fields obey

[TABLE]

The first term represents fast oscillations. We recognize the remaining terms as the contributions from the zeroth LL harmonic oscillator wavefunctions. The last term is a Dirac delta function which is typical of free Majorana fields Elliott and Franz (2015).

The real-space Majorana nature of the projected fields leading to Eq. S231 was illustrated in the continuum. However, this property is not intrinsic to this limit. Indeed, it follows from the general p-h symmetry and topology of the bulk Hamiltonian under strain as illustrated by the general arguments of Sec. I. We can therefore extend the Majorana nature to operators defined on the lattice such that $c^{{\dagger}}_{P,\sigma}(\bm{R})=-i\text{sgn}(\sigma)c_{P,\bar{\sigma}}(\bm{R})$, as mentioned in the main text. These are expected to obey a lattice version of Eq. S230.

V Estimates of strain

V.1 General Formulae

In the continuum limit, the effect of uniaxial strain of the hopping coefficients of our models is given by Nica and Franz (2018); Castro et al. (2009)

[TABLE]

where $t$ is the pristine value of the hopping. $\beta$ is given by

[TABLE]

where $a$ is the nearest-neighbor distance. It is typically associated with the Grüneisen parameter Vozmediano et al. (2010). $u_{xx}$ is a diagonal component of the strain tensor.

In the low-energy theory, $\delta t(x)$ determines an effective vector potential via

[TABLE]

where $e$ is the bare electronic charge. For illustration purposes, we chose the vector potential along the $y$ direction, but generalization to other cases of uniaxial strain are obvious.

In order to produce an approximately uniform pseudo-magnetic field, we require

[TABLE]

From the previous expression, it is clear that we require a strain which, to leading order, varies linearly with distance in the bulk of the sample:

[TABLE]

where $\rho$ is the gradient of the strain. It is convenient to define a dimensionless parameter

[TABLE]

Note that the strain parameter $\epsilon$ discussed in the main text is given by

[TABLE]

All of the relevant scales of the problems can be written in terms of the parameters $t,\Delta,\beta,\tilde{\rho}$, representing hopping coefficient, pairing amplitude, Grüneisen parameter, and rate of strain increase per unit cell, respectively. The pseudo-magnetic length is determined by

[TABLE]

The LL separation in energy is given by

[TABLE]

where we re-introduced factors of $\hbar$ in Eq. S45 with $v_{\Delta}\approx\Delta a/\hbar$.

V.2 Numerical Estimates

Using the formulae of the previous section, we can estimate typical pseudo-magnetic lengths and LL separation in energy using parameters for high-$T_{c}$ materials: $\Delta\approx$30 meV Hashimoto et al. (2014), $a\approx$ 3.9 $\AA$, $t\approx$ 0.38 eV Korshunov et al. (2004), and $\beta\approx 1.2$ Mouallem-Bahout et al. (1994).

Considering a maximum strain of 0.05 achieved over a typical length scale approximately $L=$0.05 mm, we can estimate $\tilde{\rho}$ as

[TABLE]

This gives

[TABLE]