Flux-Stabilized Majorana Zero Modes in Coupled One-Dimensional Fermi Wires

TL;DR

This paper demonstrates a new method to stabilize Majorana zero modes in coupled 1D fermionic wires by applying a synthetic magnetic flux, overcoming the destructive effects of interchain tunneling.

Contribution

It introduces a flux-based mechanism that preserves Majorana modes in 1D fermionic ladders despite significant interchain tunneling.

Findings

Majorana boundary modes persist with large interchain tunneling

Synthetic magnetic flux induces a robust leg-interchange symmetry

The model is verified using DMRG and exact diagonalization

Abstract

One promising avenue to study one-dimensional ($1$D) topological phases is to realize them in synthetic materials such as cold atomic gases. Intriguingly, it is possible to realize Majorana boundary modes in a $1$D number-conserving system consisting of two fermionic chains coupled only by pair-hopping processes. It is commonly believed that significant interchain single-particle tunneling necessarily destroys these Majorana modes, as it spoils the $\mathbb{Z}_2$ fermion parity symmetry that protects them. In this Letter, we present a new mechanism to overcome this obstacle, by piercing a (synthetic) magnetic $\pi$-flux through each plaquette of the Fermi ladder. Using bosonization, we show that in this case there exists an exact leg-interchange symmetry that is robust to interchain hopping, and acts as fermion parity at long wavelengths. We utilize density matrix renormalization group and exact diagonalization to verify that the resulting model exhibits Majorana boundary modes up to large single-particle tunnelings, comparable to the intrachain hopping strength. Our work highlights the unusual impacts of different topologically trivial band structures on these interaction-driven topological phases, and identifies a distinct route to stabilizing Majorana boundary modes in $1$D fermionic ladders.

Equations

$$H=$$

$$H_K=$$

$$-\sum _{n=0}^{L-1}(t_{\perp }{c}_{n,0}^{†}{c}_{n,1}+\text{H.c.}),$$

$$H_W=$$

$${\theta }_c$$

$${\varphi }_c$$

$$k_{L,\text{II}}+k_{R,\text{II}}-k_{L,\text{I}}-k_{R,\text{I}}=\frac{2\pi }{a}.$$

$$H$$

$$H_c$$

$$H_s$$

$$+\frac{2g_{um}}{(2\pi a{)}^2}{\int }_x\cos \left(2{\varphi }_s(x)\right)-\frac{2g_{bs}}{(2\pi a{)}^2}{\int }_x\cos \left(2{\theta }_s(x)\right)$$

$$-\frac{2g_{mx}}{(2\pi a{)}^2}{\int }_x\cos \left(2{\theta }_s(x)\right)\cdot \cos \left(2{\varphi }_s(x)\right),$$

$$u_s$$

$$K_s$$

$$\frac{d{K}_s(l)}{dl}$$

$$\frac{d{y}_{um}(l)}{dl}$$

$$\frac{d{y}_{bs}(l)}{dl}$$

$$\frac{d{y}_{mx}(l)}{dl}$$

$$H_s\sim H_{\text{s-G}}$$

$$+\frac{g_{um}}{\pi a^2}\cos \left(2{\varphi }_s(x)\right)\}.$$

$$H_K$$

$$\displaystyle=\int^{\infty}_{-\infty}dk\left(\psi^{\dagger}_{k,+}\ \ \ \psi^{\dagger}_{k,-}\right)\left(\begin{array}[]{cc}E_{+}(k,\phi)&0\\ 0&E_{-}(k,\phi)\\ \end{array}\right)\left(\begin{array}[]{c}\psi_{k,+}\\ \psi_{k,-}\\ \end{array}\right),$$

$$\mathbb{M}(k,\phi)=\left(\begin{array}[]{cc}|\chi_{+}\rangle&|\chi_{-}\rangle\\ \end{array}\right)=\left(\begin{array}[]{cc}\frac{-u_{k,\phi}}{\sqrt{1+u_{k,\phi}^{2}}}&\frac{1}{\sqrt{1+u_{k,\phi}^{2}}}\\[1.99997pt] \frac{1}{\sqrt{1+u_{k,\phi}^{2}}}&\frac{u_{k,\phi}}{\sqrt{1+u_{k,\phi}^{2}}}\\ \end{array}\right),\ \ \ \mathbb{M}^{\dagger}\mathbb{H}(k,\phi)\mathbb{M}=\left(\begin{array}[]{cc}E_{+}(k,\phi)&0\\ 0&E_{-}(k,\phi)\\ \end{array}\right),$$

$$u_{k,\varphi }=\frac{1}{t_{\perp }}\left\{2t_{\parallel }\sin (ka)\sin \frac{\varphi }{2}+\sqrt{t_{\perp }^2+{\left(2t_{\parallel }\sin (ka)\sin \frac{\varphi }{2}\right)}^2}\right\}.$$

$${\psi }_{x,\ell }$$

$$+e^{i{k}_{L,\text{II}}x}⟨\ell |{\chi }_{-}(k_{L,\text{II}})⟩{\tilde{\psi }}_{L,\text{II}}(x)+e^{i{k}_{R,\text{II}}x}⟨\ell |{\chi }_{-}(k_{R,\text{II}})⟩{\tilde{\psi }}_{R,\text{II}}(x),$$

$$u_c$$

$$H^{\prime }=$$

$$+\sum _{n=0}^{L-2}\left(W{c}_{n,0}^{†}{c}_{n+1,0}^{†}{c}_{n,1}{c}_{n+1,1}+\text{H.c.}\right),$$

$$H^{\prime }=$$

$$+\sum _{n=0}^{L-2}\left(W{c}_{n,0}^{†}{c}_{n+1,0}^{†}{c}_{n,1}{c}_{n+1,1}+\text{H.c.}\right)+\sum _{n=0}^{L-1}\left(h{c}_{n,0}^{†}{c}_{n,0}-h{c}_{n,1}^{†}{c}_{n,1}\right),$$

$$H^{\prime }=$$

$$+\sum _{n=0}^{L-2}\left(W{c}_{n,0}^{†}{c}_{n+1,0}^{†}{c}_{n,1}{c}_{n+1,1}+\text{H.c.}\right),$$

Full text

The classification of topological phases in one dimension Pollmann ; FidkowskiKitaevClass ; ChenGuWen revealed an intriguing array of possible new states of matter, with a variety of types of protected gapless boundary modes. Of these, one class that has generated considerable excitement recently is the one-dimensional ($1$D) topological superconductors first described by Ref. Kitaev . These have protected boundary Majorana zero modes, which harbor non-Abelian statistics AliceaBraiding ; SauPRB ; Hassler and hence are promising candidates for topological quantum computing NayakReview ; SarmaReview .

In practise, to obtain long-range superconducting order in $1$D systems requires inducing superconductivity via coupling to a $3$D bulk superconductor LFu ; Fu ; AliceaPRB ; JDSau ; Lutchyn ; Oreg ; Vazifeh ; Klinovaja ; Braunecker ; Alicea ; Clarke ; experimental progress in this direction has been made in several distinct solid-state systems Mourik ; Das ; Deng ; NPerge ; Albrecht . Interestingly, however, it is also possible to host topological boundary modes in truly $1$D platforms, in spite of the fact that these systems do not support long-range superconducting order, and are in fact gapless Fidkowski ; Sau ; Kraus ; Cheng ; Lang ; Iemini . This opens up the possibility of studying $1$D fermionic topological phases in synthetic materials such as cold atomic gases Zhang ; Sato ; Chen ; Qu ; WZhang ; Liu ; BlochRMP ; Goldman , offering an attractive architecture with increased tunability.

One concrete model in this category was proposed by Ref. Kraus , who showed numerically that even starting from a lattice Hamiltonian with a topologically trivial band structure, a regime bearing the hallmarks of Majorana boundary modes can be accessed in an atomic two-leg ladder by introducing an interleg pair-hopping interaction. These boundary modes are protected by the conserved fermion parity of one of the wires, and are therefore robust provided that the single-particle interleg tunneling ($t_{\perp}$), which breaks this symmetry explicitly, is sufficiently small. A number of proposals Kraus ; Cheng ; Lang ; Iemini for suppressing $t_{\perp}$ such that this regime may be experimentally realized ensued.

In this paper we propose a distinct route to overcome this obstacle: We begin with a different band structure, in which each plaquette of the ladder has a flux of $\pi$. We show that with pair hopping this model also has an interacting topological regime hosting Majorana boundary modes—which in this case are protected by a symmetry preserved even in the presence of finite $t_{\perp}$. We bolster our theory with numerical evidence of Majorana boundary modes over a wide range of $t_{\perp}$. Our findings thus furnish an appealing mechanism to engineer topological boundary modes in a particle-conserving and strongly interacting system without the need to fine-tune single-particle tunneling.

Fermionic flux ladder model.—Motivated by these considerations, we study an interacting two-leg ladder model of spinless fermions in a perpendicular magnetic field described by the following number-conserving Hamiltonian,

[TABLE]

where $c^{(\dagger)}_{n,\ell}$ is the fermionic annihilation $($creation$)$ operator at rung $n$ on the leg $\ell\!=\!0,1$. The intraleg and interleg single-particle tunneling strengths are $t_{\parallel}$ and $t_{\perp}$, respectively, and $W$ (which we take to be negative throughout this work) is the pair-hopping strength. Two essential ingredients of the above model are the synthetic Peierls phase $\phi\!\in\![0,\pi]$ per plaquette, and the interchain pair-hopping interaction $H_{W}$. Previous works have demonstrated the existence of Majorana boundary modes in this Hamiltonian at $\phi\!=\!0$ and $t_{\perp}\!=\!0$ based on a preserved fermion number parity ${\sf P}_{\ell}\!\coloneqq\!(-1)^{N_{\ell}}$, where $N_{\ell}$ is the particle-number operator of a single leg $\ell$ Cheng ; Kraus ; Lang ; Iemini . Here we will show that when $\phi\!=\!\pi$, these boundary modes persist up to $t_{\perp}$ of order $t_{\parallel}$. Note that without pair hopping, the bare band $H_{K}$ in (2) is topologically trivial, so that the model (1) requires interactions to realize the topological regime. This is necessarily the case for isolated $1$D systems, where the total fermion number is conserved—in contrast to models based on e.g., Kitaev or spin-orbit-coupled wires Stoudenmire ; Chan ; FidkowskiKitaev , where the Majorana modes may originate from nontrivial Bogoliubov–de Gennes band structures, as Cooper pairs can be exchanged with a bulk $3$D superconductor.

Symmetry analysis.—Though our system is gapless, its topological boundary modes can be understood using the symmetry classification of $1$D gapped fermionic phases FidkowskiKitaevClass . We therefore begin with a discussion of the underlying symmetries. Below we argue that the Majorana boundary modes are protected by a unitary $\mathbb{Z}_{2}$ leg-interchange symmetry ${\sf L_{s}}$ which takes $c_{n,0}\!\rightarrow\!(-1)^{n+1}c_{n,1},\ c_{n,1}\!\rightarrow\!(-1)^{n+1}c_{n,0}$, where $n$ indexes the site along the chain. This is a symmetry only for $\phi\!=\!\pi$; we will see that its action is equivalent to that of an emergent fermion parity operator that we will describe presently using bosonization. Additionally, for any flux $\phi$, there is an antiunitary time-reversal symmetry that combines complex conjugation with interchanging the two legs of the ladder, obeying ${\sf T}^{2}\!=\!+1$. In principle, this enables eight distinct types of boundary modes FidkowskiKitaevClass , though numerically we observe only one of these, with a single Majorana zero mode at each boundary. Finally, there also exists an overall $U(1)$ fermion number conservation. In our case its chief significance is that, unlike the situation considered by Ref. FidkowskiKitaevClass , only at exactly half-filling is a microscopic (rather than emergent) particle-hole symmetry possible. If $t_{\perp}\!=\!0$, the model has an additional exact $\mathbb{Z}_{2}$ symmetry, corresponding to the fermion parity of a single leg of the ladder, given by ${\sf P}_{\ell}$. It is this symmetry that protects the topological boundary modes observed by Ref. Kraus at zero flux.

Bosonization and renormalization group analysis.—To understand the special role of the $\mathbb{Z}_{2}$ symmetry ${\sf L_{s}}$, we bosonize the model (1). For $\phi\!=\!\pi$, the kinetic Hamiltonian $H_{K}$ has band energies $E_{\textrm{hgr/lwr}}\!=\!\pm\!\sqrt{t^{2}_{\perp}\!+\!4t^{2}_{\parallel}\sin^{2}\!\left(ka\right)}$, with band gap $2t_{\perp}$. Here we consider a system at less than half-filling, with the interaction scale $W$ small relative to the bandwidth, such that we can project out the unoccupied band and focus only on the processes inside the lower band when treating $H_{W}$. Here the Fermi energy intersects the lower band at four separate Fermi points, as shown in Fig. 1.

Linearizing about these four Fermi points results in two right-moving ($R$) and two left-moving ($L$) fermion operators, which we distinguish using a valley index (I or II). Bosonizing these in the usual way, we have $\psi_{\kappa,\nu}\!\sim\!e^{i\varphi_{\kappa,\nu}}$, with $\kappa\!=\!R\ (L)$ and $\nu\!=\!\textrm{I}\ (\textrm{II})$. We define the nonchiral bosonic fields

[TABLE]

where $\theta_{\nu}\!=\!\frac{1}{\sqrt{2}}(\varphi_{R,\nu}\!-\!\varphi_{L,\nu})$, $\phi_{\nu}\!=\!\frac{1}{\sqrt{2}}(\varphi_{R,\nu}\!+\!\varphi_{L,\nu})$, such that $\varphi_{\kappa,\nu}\!=\!\frac{1}{2}[\left(\phi_{c}\!+\!\kappa\theta_{c}\right)\!+\!\nu\!\left(\phi_{s}\!+\!\kappa\theta_{s}\right)]$. After including appropriate Klein factors, the only nontrivial commutators among these nonchiral bosonic fields are $\left[\theta_{c}(x),\phi_{c}(x^{\prime})\right]\!=\!-2i\pi\Theta(x^{\prime}\!-\!x),\ \left[\theta_{s}(x),\phi_{s}(x^{\prime})\right]\!=\!2i\pi\Theta(x\!-\!x^{\prime}),\ \left[\theta_{s}(x),\phi_{c}(x^{\prime})\right]\!=\!-2i\pi$, where $\Theta(x)$ is the Heaviside step function. The corresponding density and current operators are: $\rho_{c}\!=\!(1/\pi)\partial_{x}\theta_{c},\ \rho_{s}\!=\!(1/\pi)\partial_{x}\theta_{s},\ J_{c}\!=\!(1/\pi)\partial_{x}\phi_{c},\ J_{s}\!=\!(1/\pi)\partial_{x}\phi_{s}$.

Bosonizing the interaction term $H_{W}$ produces multiple four-fermion terms, of which only slowly-varying terms contribute in the continuum limit. When $\phi\!=\!\pi$, in addition to the usual momentum-conserving processes, this allows for intervalley umklapp scattering, since the Fermi points obey $k_{L,\textrm{I}}\!=\!-k_{R,\textrm{II}},\ k_{R,\textrm{I}}\!=\!-k_{L,\textrm{II}},\ k_{L,\textrm{II}}\!+\!k_{R,\textrm{II}}\!=\!\pi/a$, such that

[TABLE]

Eq. (6) is valid independent of the chemical potential within the lower band, but it will not hold if $\phi\neq\pi$.

The resulting bosonized form of $H$ decouples into a gapless charge sector and a gapped spin sector:

[TABLE]

where $\int_{x}\!\equiv\!\int\!dx$ and the associated coupling constants are given by $g_{um}\!=\!-\frac{aW}{\pi^{2}}\!\cos^{2}(k_{R,\textrm{II}}a)(\sin^{4}\!\frac{\xi}{2}\!+\!\cos^{4}\!\frac{\xi}{2}),\ g_{bs}\!=\!-\frac{aW}{2\pi^{2}}\!\sin^{2}(k_{R,\textrm{II}}a)\sin^{2}\!\xi,\ \mbox{and}\ g_{mx}\!=\!-\frac{aW}{2\pi^{2}}(\sin^{4}\!\frac{\xi}{2}\!+\!\cos^{4}\!\frac{\xi}{2})$. Here the wavefunction in the lower band at $k_{R,\textrm{II}}$, has the form $(\cos\!\frac{\xi}{2},\ \sin\!\frac{\xi}{2})$ in the leg (i.e., $0,1$) basis SupplementalMaterial . Note that the $g_{um}$-term favors Cooper pairing while the competing $g_{bs}$-term favors a valley density-wave-type order. Moreover, owing to the density-density interactions, the velocities and Luttinger parameters are renormalized. Particularly, for the spin channel

[TABLE]

where u\!=\!\frac{1}{2}v_{F}\!\equiv\!\frac{1}{2}\frac{dE_{\textrm{lwr}}}{dk}\big{|}_{k=k_{R,\textrm{II}}}\!>\!0 and $g\!=\!\frac{aW}{2(2\pi)^{3}}\!\sin^{2}\!\xi$.

As the charge sector is gapless, we will concentrate on the spin sector, which is responsible for the topological Majorana boundary modes. To determine which of the sine-Gordon terms in $H_{s}$ dominates, we use the one-loop renormalization group (RG) flow equations:

[TABLE]

where $y_{um}\!=\!\frac{g_{um}}{2u\pi},\ y_{bs}\!=\!\frac{g_{bs}}{2u\pi},$ and $y_{mx}\!=\!\frac{g_{mx}}{2u\pi}$ are the dimensionless coupling constants. The $g_{mx}$-channel, involving both $\theta_{s}$ and $\phi_{s}$, is power-counting irrelevant for any value of $K_{s}$ and can be neglected. Accordingly, the gap-opening competition is between $g_{um}$ and $g_{bs}$. If we take $W$ to be negative, then for $3\pi/(4a)\!>\!k_{R,\textrm{II}}\!>\!\pi/(2a)$ we have $K_{s}\!>\!1$ and $g_{um}$ is the only relevant coupling. In this regime the long-wavelength effective Hamiltonian is therefore expected to be

[TABLE]

Note that for $W\!<\!0,\ g_{um}\!>\!0$ and the coefficient of the $\cos\!\left(2\phi_{s}\right)$ term is positive, contrary to the usual sine-Gordon model. As explained below, this sign flip, which arises from the particular umklapp processes that occur in the band structure for $\phi\!=\!\pi$, is crucial to ensuring that the staggered leg-interchange symmetry ${\sf L_{s}}$ acts like a fermion parity operator. Indeed, for $W\!>\!0$ and $3\pi/(4a)\!>\!k_{R,\textrm{II}}\!>\!\pi/(2a)$ it is the negative coupling $g_{bs}$ that is most relevant; this phase shows no numerical signatures of topological boundary modes, as discussed in the Supplemental Material SupplementalMaterial .

$\mathbb{Z}_{2}$* symmetry and Majorana operators.*—Following Ref. Clarke we can further derive the bosonized Majorana-boundary-mode operators by specifying the vacuum Hamiltonian density ${\sf H_{vac}}(x)\!=\!-\frac{2M_{\infty}}{\pi a}\!\sin\!\left(\theta_{s}(x)\right)\!\cos\!\left(\theta_{c}(x)\right)$, which amounts to sending the fermion mass outside the system to $+\infty$ on both legs of the ladder Keselman . This fixes $\theta_{s}(x\!\in\!(-\infty,-\frac{L}{2}))\!=\!\pm\frac{\pi}{2}\!+\!2\hat{n}^{(1)}_{\theta_{s}}\pi$ and $\theta_{s}(x\!\in\!(\frac{L}{2},\infty))\!=\!\pm\frac{\pi}{2}\!+\!2\hat{n}^{(2)}_{\theta_{s}}\pi$. In the gapped bulk of the ladder, $\phi_{s}(x\!\in\!(-\frac{L}{2},\frac{L}{2}))\!=\!-\frac{\pi}{2}\!+\!\hat{n}_{\phi_{s}}\pi$, where $\hat{n}^{(1,2)}_{\theta_{s}},\ \hat{n}_{\phi_{s}}$ are integer operators. The two domain-wall Majorana operators are then: $\gamma_{L}\!\simeq\!e^{i\pi(\hat{n}^{(1)}_{\theta_{s}}+\hat{n}_{\phi_{s}})},\ \gamma_{R}\!\simeq\!e^{i\pi(\hat{n}^{(2)}_{\theta_{s}}+\hat{n}_{\phi_{s}})}$. These satisfy the Majorana relations $\gamma^{\dagger}_{L/R}\!=\!\gamma_{L/R},\ \gamma^{2}_{L/R}\!=\!1,\ \{\gamma_{L},\gamma_{R}\}\!=\!0$, and define an emergent fermion parity ${\sf P_{v}}\!=\!i\gamma_{L}\gamma_{R}$. This acts on the spin fields within the gapped system via ${\sf P_{v}}\phi_{s}{\sf P^{-1}_{v}}\!=\!\phi_{s}\!-\!\pi,\ {\sf P_{v}}\theta_{s}{\sf P^{-1}_{v}}\!=\!\theta_{s}$. In particular, when $g_{um}$ dominates the RG flow, ${\sf P_{v}}$ maps between the system’s two classical ground states; when $g_{bs}$ dominates the flow, it acts trivially on the two classical ground states.

Strictly speaking the operator ${\sf P_{v}}$ does not correspond to any microscopic symmetry for $|t_{\perp}|\!>\!0$. (For $t_{\perp}\!=\!0$, it can be interpreted as the fermion parity of one of the ladder’s two legs.) However, at $\phi\!=\!\pi$ its action is equivalent to that of the microscopic $\mathbb{Z}_{2}$ symmetry ${\sf L_{s}}$, which acts on the bosonized field via ${\sf L_{s}}\phi_{s}(x){\sf L^{-1}_{s}}\!=\!-\phi_{s}(x)\ (\mbox{mod}\ 2\pi)$. Specifically, since $g_{um}$ is positive, both ${\sf P_{v}}$ and ${\sf L_{s}}$ map between the two classical minima $\left|\pm\pi/2\right\rangle$ of the sine-Gordon potential, where $\phi_{s}(x)\left|\pm\pi/2\right\rangle\!=\!\left(\pm\pi/2\ \mbox{mod}\ 2\pi\right)\left|\pm\pi/2\right\rangle$. Thus within the ground-state manifold, the action of ${\sf P_{v}}$ is equivalent to that of the exact symmetry of ${\sf L_{s}}$. Indeed, in a finite-size system, due to the proliferation of instanton processes ChenBurnell ; Fidkowski , the ground states of $H_{\textrm{s-G}}$ are split into $\left|\pm\right\rangle\!=\!\frac{1}{\sqrt{2}}\!\left(\left|\pi/2\right\rangle\!\pm\!\left|-\pi/2\right\rangle\right)$ whose energy difference is exponentially small in the system’s size. The lowest two eigenstates therefore obey ${\sf P_{v}}\left|\pm\right\rangle\!=\!\pm\left|\pm\right\rangle,\ {\sf L_{s}}\left|\pm\right\rangle\!=\!\pm\left|\pm\right\rangle$, and can be labeled by their eigenvalues under the microscopic symmetry ${\sf L_{s}}$.

We note that for flux $\phi\!\neq\!\pi$, the umklapp scattering leading to Eq. (16) is no longer an effective zero (lattice)-momentum transfer process, and $\cos\!\left(2\phi_{s}(x)\right)$ is replaced by $\cos\!\left(2\phi_{s}(x)\!-\!\delta x\right)$, where $\delta$ parameterizes the deviation of the flux from $\pi$. In this case we expect a commensurate-incommensurate transition which destroys the topological boundary modes. As described in the Supplemental Material SupplementalMaterial , numerical evidence suggests that this transition leads to a phase in which $\theta_{s}$ is locked.

Comparing $\phi\!=\!\pi$ with $\phi\!=\!0$.—As noted above, the topological state found in Ref. Kraus at $\phi\!=\!0,\ t_{\perp}\!=\!0$ is protected by the symmetry ${\sf P_{\ell}}$ corresponding to the conserved fermion parity of one of the ladder’s legs. For $t_{\perp}\!=\!0$, the operator ${\sf P_{v}}$ constructed from the Majorana boundary modes is exactly equal to ${\sf P_{\ell}}$, and therefore corresponds to an exact symmetry that protects the resulting topological boundary modes. In comparison, at $\phi\!=\!\pi$ we have found that ${\sf P_{v}}$’s action on the classical ground states is identical to that of the leg-exchange symmetry ${\sf L_{s}}$, which is not violated by finite $t_{\perp}$. Hence it is the particular action of ${\sf L_{s}}$ at this point which gives the topological boundary modes their enhanced stability. This result highlights the fact that different choices of topologically trivial band structures can have profound implications for interaction-driven topological phases.

In addition to these symmetry considerations, introducing $t_{\perp}$ at $\phi=0$ has a significantly different impact on the band structure and umklapp processes than doing so at $\phi\!=\!\pi$, making the latter state more stable. For $\phi\!=\!t_{\perp}\!=\!0$, the two chains’ bands are identical. This overlap favors the umklapp-scattering processes that lead to the topological phase. However, increasing $t_{\perp}$ at $\phi\!=\!0$ separates the two bands and creates a Fermi-velocity mismatch, which disfavors these processes, ultimately causing the spin gap to close at a small but finite value of $t_{\perp}$ Kraus . By contrast, when $\phi\!=\!\pi$, the two valleys where the chains’ bands cross the Fermi surface are maximally separated by a wave vector $\pi/a$, and remain symmetric with a finite $t_{\perp}$. In this case the Fermi-velocity mismatch is absent, and the spin gap persists up to large values of $t_{\perp}$. Thus the different nontopological band structures play a key role in making the topological Majorana modes robust (fragile) against $t_{\perp}$ at $\phi\!=\!\pi$ ($\phi\!=\!0$), consistent with the fact that $t_{\perp}$ breaks ${\sf P_{\ell}}$ but preserves ${\sf L_{s}}$.

Numerical verification.—Simulations based on density matrix renormalization group (DMRG) White and exact diagonalization (ED) Sandvik have been performed to solve the lattice model (1) at $\phi\!=\!\pi$ and $t_{\perp}\!=\!0.5t_{\parallel}$. The numerical outcomes provide strong evidence supporting our theoretical predictions. Fig. 2(a) demonstrates that in the low-population region ($N/L\!=\!1/3$), when pair hopping is strong ($W\!=\!-1.7t_{\parallel}$), the energy gap between the ground state and the $1$st excited state closes exponentially as the ladder’s size $L$ increases. This is in contrast with a power-law gap closing resulting from the gapless charge sector. As anticipated, these two nearly degenerate eigenstates are distinguished by their eigenvalues of ${\sf L_{s}}$: The ground ($1$st excited) state has eigenvalue $+1$ ($-1$). The resulting ground-state manifold is further separated from the rest of the spectrum by a gap which only decreases inversely with $L$ (Fig. 2(b)). Moreover, the topological Majorana boundary modes can be characterized via the nonlocal correlations Kraus ; Lang ; Iemini in the single-particle Green functions $G_{mn}\!\coloneqq\!\langle c^{\dagger}_{m,0}c_{n,0}\rangle\!=\!\langle c^{\dagger}_{m,1}c_{n,1}\rangle$. These are apparent for a range of $t_{\perp}$ values in the topological regime (Figs. 2(c)–(d)). The presence of the edge states also gives rise to a two-fold degeneracy in the entanglement spectrum (ES) Li ; Pollmann on the central bond (Fig. 2(f)). These DMRG results have been confirmed by ED for small system sizes.

As discussed in the Supplemental Material SupplementalMaterial , adding a small ${\sf L_{s}}$-symmetry-breaking perturbation to (1) causes the degeneracy in the ES to lift, and the energy gap of $E_{1}\!-\!E_{0}$ to deviate from the exponentially small splitting observed in Fig. 2(a). This supports our claim that the topological boundary modes are protected by ${\sf L_{s}}$ at $\phi\!=\!\pi$ and $t_{\perp}\!\neq\!0$.

To fully characterize the topological regime, we also show evidence of a phase transition out of the topological regime at $t_{\perp,\textrm{c}}\!\approx\!2.5t_{\parallel}$. Fig. 2(e) shows a value $t_{\perp}\!>\!t_{\perp,\textrm{c}}$, for which the nonlocal correlations in $G_{mn}$ are absent. Further, as shown by Fig. 2(f), the degeneracy in the ES apparent for $t_{\perp}\!<\!t_{\perp,\textrm{c}}$ is lifted and the lowest level is clearly nondegenerate at $t_{\perp}\!=\!3.0t_{\parallel}$. Our numerics suggest that the transition is toward a state with density-wave order at large $t_{\perp}$: For $t_{\perp}\!>\!t_{\perp,\textrm{c}}$, the fermion-density profile $\rho_{n}\!\coloneqq\!\langle c^{\dagger}_{n,0}c_{n,0}\rangle\!=\!\langle c^{\dagger}_{n,1}c_{n,1}\rangle$ in Fig. 2(g) evolves from a uniform distribution toward an oscillatory pattern. Additionally, for $t_{\perp}\!\geq\!3.0t_{\parallel}$, the two lowest energy eigenstates split and both have ${\sf L_{s}}$ eigenvalues of $+1$ SupplementalMaterial . All the above observations suggest that once $t_{\perp}\!>\!t_{\perp,\textrm{c}}$, the Majorana boundary modes disappear.

Experimental feasibility.—In cold atom laboratories, currently a great deal of effort has been devoted to simulating strong synthetic magnetic fields in optical lattices by Raman-assisted tunnelings and lattice-shaking techniques, which mimic standard Peierls substitutions through attaching an Aharanov–Bohm-like complex phase to the nearest-neighbor single-particle hopping Aidelsburger ; AidelsburgerFlux ; MiyakeFlux ; Atala ; Mancini ; Kennedy ; Struck ; Livi . The remarkable tunability of these artificial gauge fluxes has paved an experimental avenue to exploring the Bose–Einstein condensation in Harper–Hofstadter model of $2$D bosonic gases Kennedy . Chiral edge states due to particle’s cyclotron motion stirred by uniform fluxes have also been detected in quasi-$1$D fermionic ladders and Hall ribbons, where a sizable flux $\phi$ can reach up to $1.31\pi$ per plaquette by utilizing the technology of optical atomic clocks Livi ; Mancini . In Ref. Kraus an atomic scheme has also been outlined for creating the more intricate pair-hopping interaction. These developments make the prospect of realizing flux-stabilized Majorana zero modes a possibility in the near future.

To summarize, we have shown that by threading $\pi$-flux through each plaquette, interaction-driven Majorana bound states in fermionic ladders may be stabilized in the presence of single-particle interleg tunneling. En route we have established a connection between a microscopic $\mathbb{Z}_{2}$ leg-exchange symmetry present only at this flux value and the action of an emergent fermion parity operator in the long-wavelength bosonized theory. We have also highlighted the advantages of the $\pi$-flux state in fostering umklapp processes which generate the spin gap enabling this topological regime. Our theory has been substantiated by extensive DMRG and ED calculations.

C. C. thanks M. D. Schulz for discussions. Part of the simulation is developed from ALPS package Dolfi . C. C. and F. J. B. are supported by NSF-DMR $1352271$ and by the Sloan Foundation FG-$2015$-$65927$. W. Y. and Y. C. are supported by the National Natural Science Foundation of China (Grant Nos. $11625416$, $11474064$, and $11274069$). C. S. T. is supported by the Robert A. Welch Foundation under Grant No. E-$1146$.

C. C. and W. Y. contributed equally to this work.

Appendix A Band structure

Here we provide more details on the band structure of $H_{K}$, and show how it enters the bosonized form of the full Hamiltonian. The noninteracting kinetic Hamiltonian,

[TABLE]

can be diagonalized by introducing the unitary matrix $\mathbb{M}(k,\phi)$,

[TABLE]

where the component in the matrix elements assumes

[TABLE]

From now on, we focus on the case of $\phi=\pi$, with the chemical potential in the lower band. Linearizing the band structure about the four Fermi points in the lower band gives rise to four independent but mutually interacting chiral-fermion branches, which we label $(L,\textrm{I}),\ (R,\textrm{I}),\ (L,\textrm{II}),$ and $(R,\textrm{II})$ (see Fig. 1 in the main text). The fermion operator in real space is related to these four chiral fermions via:

[TABLE]

where $x$ is the position along the ladder, and $\ell=0,1$ indexes the leg. Here we have defined $\langle\ell=0|\equiv\left(1\ \ 0\right)$, $\langle\ell=1|\equiv\left(0\ \ 1\right)$, and |\chi_{-}(k)\rangle\equiv|\chi_{-}(k,\phi=\pi)\rangle=\left(\begin{array}[]{c}\frac{1}{\sqrt{1+u^{2}_{k,\pi}}}\\[1.99997pt] \frac{u_{k,\pi}}{\sqrt{1+u^{2}_{k,\pi}}}\\ \end{array}\right)=\left(\begin{array}[]{c}\cos\frac{\xi(k)}{2}\\[3.00003pt] \sin\frac{\xi(k)}{2}\\ \end{array}\right).

The Hamiltonian in the main text is obtained by substituting these expressions into the pair-hopping term $H_{W}$, and performing a standard bosonization analysis of the result. Following these calculations through, one finds that the matrix elements $\langle\ell|\chi_{-}\rangle$ enter the resulting Luttinger parameters, as well as the sine-Gordon couplings $g_{um},\ g_{bs},\ g_{mx}$. The relevant expressions in the spin sector are given in the main text; in the charge sector the associated Luttinger parameters are defined by

[TABLE]

where the matrix element enters $u_{c}$ and $K_{c}$ via $g\!=\!\frac{aW}{2(2\pi)^{3}}\sin^{2}\xi$.

It follows that the microscopic parameters in our bosonized model are sensitive to the precise initial band structure, as well as to the initial value of the chemical potential within the lower band. In particular, the wrong choice of chemical potential can lead to $K_{s}<1$, and an RG flow that tends to take the system away from the topological regime.

Appendix B Topological quantum phase transition

Though at $\phi=\pi$ our model is robust to the single-particle interchain tunneling up to $t_{\perp}\approx t_{\parallel}$, Figs. 2(c)–(e) in the main text suggest that by $t_{\perp}=3t_{\parallel}$, the system has undergone a transition to a nontopological regime. In particular, for this value we see that the nonlocal correlation in the single-particle Green function, an indicator of the presence of Majorana boundary modes, disappears. Further, the double degeneracy in the entanglement spectrum, another indicator of the topological boundary modes, disappears by $t_{\perp}=2.5t_{\parallel}$. Finally, the spatial profile of the local fermion density (shown in Fig. 2(g) in the main text) shows clear evidence of a fermion density-wave-type order for $t_{\perp}=3t_{\parallel}$. In comparison, the bulk fermion density is uniform for $t_{\perp}\leq t_{\parallel}$, where the nonlocal correlation indicates the presence of Majorana boundary modes.

As further evidence of this transition, Fig. 3 shows the low-energy spectrum of the $\pi$-flux ladder as a function of $t_{\perp}/t_{\parallel}$ at fixed pair-hopping strength $W$. At approximately $t_{\perp}=2.8t_{\parallel}$ a level crossing occurs between two states with opposite $\mathbb{Z}_{2}$ quantum numbers under the staggered leg-interchange symmetry ${\sf L_{s}}$. For smaller $t_{\perp}$ the two lowest energy states have opposite $\mathbb{Z}_{2}$ quantum numbers; whereas for larger $t_{\perp}$ both have the same ($+1$) $\mathbb{Z}_{2}$ eigenvalue. As explained in the main text this $\mathbb{Z}_{2}$ quantum number is associated with the fermion parity in the topological region, suggesting that this level crossing indicates a transition out of the topological regime. This level crossing is consistent with what we expect for a transition from a topological regime to a fermion-density-wave order: In the latter case we expect the lowest two states to be related by translation, and hence to have the same parity under the leg-exchange symmetry ${\sf L_{s}}$.

Appendix C The pi-flux ladder for positive $W$

The RG analysis in the main text suggests that for finite $t_{\perp}$ when $W>0,\ \phi=\pi$, the most relevant gapping channel of $H_{s}$ for low fermion density $\left(3\pi/(4a)>k_{R,\textrm{II}}>\pi/(2a)\right)$ is the backscattering term $g_{bs}$. In this section we present the corresponding DMRG results to illustrate that the resulting phase exhibits no signatures of topological protection at finite $t_{\perp}$.

In particular, we evaluate the entanglement spectra (ES) for $W=+1.7t_{\parallel}$ at several different values of $t_{\perp}$. The results in Fig. 4(a) demonstrate the significant splitting of the ES two-fold degeneracy at the lowest level even for the smallest nonzero values of $t_{\perp}\approx 0.1t_{\parallel}$. This is in sharp contrast to the situation of $W=-1.7t_{\parallel}$, where the even degeneracy in ES persists up to $t_{\perp}\approx 2.0t_{\parallel}$ (see Fig. 2(f) in the main text). Further, for parameters where this degeneracy is lifted, the local fermion-density profile $\rho_{n,\ell}$ exhibits weak modulations along the chain, as can be seen from Fig. 4(b). Here we define $\rho_{n,0}=\langle c^{\dagger}_{n,0}c_{n,0}\rangle,\ \rho_{n,1}=\langle c^{\dagger}_{n,1}c_{n,1}\rangle$ for the upper and lower chains, respectively.

Appendix D Leg-interchange-symmetry-breaking perturbations

Here, we demonstrate numerically that the $\mathbb{Z}_{2}$ leg-exchange symmetry and the flux $\phi=\pi$ are both central to protecting the topological features of our model. Specifically, we investigate three different ways of breaking the leg-interchange symmetry, and show that small (but finite) values of all symmetry-breaking perturbations lift the degeneracies in the ES associated with the topological boundary modes. In addition, we perform an analysis of the effect of perturbing the flux away from $\phi=\pi$, which lifts the degeneracies in the entanglement spectrum and leads to a fermion-density-wave order.

D.1 Making the intrachain hoppings different

We first consider the result of explicitly breaking the leg-interchange symmetry by choosing $t_{\parallel}$ to be different on the two legs of the ladder. We thus modify the Hamiltonian in the following way:

[TABLE]

where we have set $\phi=\pi$, and the added $\delta_{t}$-terms explicitly break the ${\sf L_{s}}$ symmetry.

The arguments given in the main text assert that breaking ${\sf L_{s}}$ should lift the degeneracy associated with the topological boundary modes, since there is no longer any reason for the fermion parity even and fermion parity odd states to be proximate in energy. This is apparent in the evolution of the ground-state entanglement spectra (shown in Fig. 5(a) for a system of $L=36,\ N=12$) as a function of $\delta_{t}/t_{\parallel}$. We observe that for fixed $W=-1.7t_{\parallel},\ t_{\perp}=t_{\parallel}$, the two-fold degeneracy in the lowest level of the ES has completely disappeared by $\delta_{t}=0.15t_{\parallel}$. This suggests that the symmetry ${\sf L_{s}}$ is, as claimed, integral to protecting the topological boundary modes. Fig. 5(b) further supports this, by showing that for sufficiently large $\delta_{t}$ ($\delta_{t}=0.2t_{\parallel}$), the fermion density on the lower chain tends to develop weak spatial modulations, similar to those observed for positive $W$ in the previous section. For smaller values ($\delta_{t}=0.05t_{\parallel}$), the density profiles appear to remain uniform in the bulk, though the relative populations of the two chains are generically different when the leg-interchange symmetry is broken.

It is somewhat surprising that for the smallest values of $\delta_{t}/t_{\parallel}$ the degeneracy in the ES and the uniform fermion density persist, suggesting that the topological boundary modes are only destroyed at finite values of the symmetry-breaking parameter. A similar effect was observed in Ref. Kraus , where the authors saw signatures of topological features in the $\phi=0$ chain for small but finite interchain single-particle hopping—the relevant symmetry-breaking perturbation in that case. This may be because the spin sector is gapped and completely decoupled from the gapless charge sector; hence its behavior is robust to sufficiently small perturbations.

D.2 Making imbalanced local potentials

As a second test, we consider adding a chemical potential that is imbalanced between the two chains of the Fermi ladder. We thus take the Hamiltonian to be

[TABLE]

where the last $h$-terms also break the ${\sf L_{s}}$ symmetry explicitly.

Via DMRG, we find that the effect of $h$-terms (see Fig. 6) is similar to that of the $\delta_{t}$-terms, both of which tend to destroy the topological Majorana boundary modes as indicated by the splitting (or nondegeneracy) in the lowest level of ES. This perturbation also has a similar effect on the fermion-density profile, first changing the relative occupancies of the two chains, and at higher values leading to density modulations in the higher density chain. The two perturbations differ quantitatively, however: In this case, the two-fold ES degeneracy survives up to $h/t_{\parallel}$ between about $0.2$ and $0.25$. The other parameters are kept the same as that in Fig. 5: $W=-1.7t_{\parallel},\ t_{\perp}=t_{\parallel},\ \phi=\pi,\ L=36,\ N=12$.

D.3 Making the flux deviate from pi

Finally, we check numerically the consequences of the deviations of $\phi$ from the specific value of $\pi$. The corresponding Hamiltonian takes the following form,

[TABLE]

where the extra $\delta_{\phi}$-phases break the ${\sf L_{s}}$ symmetry explicitly.

Qualitatively, our results are similar to the two cases described above: A small perturbation $\delta_{\phi}$ is sufficient to lift the degeneracy of the ES, signaling a loss of the topological boundary modes. However, as demonstrated in Fig. 7(a), for $t_{\perp}/t_{\parallel}=1$ there is a significant quantitative difference: The topological ES degeneracy is significantly more sensitive to the $\delta_{\phi}$-perturbations than to the $\delta_{t}$- and $h$-terms. The lowest-level two-fold degeneracy in the ES becomes split at a small value of $\delta_{\phi}/\pi\lesssim 0.03$, indicating the instability of the Majorana boundary modes. We note, however, that this quantitative effect is necessarily sensitive to the magnitude of $t_{\perp}$, since for $t_{\perp}=0$ the flux has less impact on the band structure. Thus for small interchain single-particle tunneling strengths, the ES degeneracy survives to larger values of $\delta_{\phi}/\pi$.

As discussed in the main text, as the flux deviates from $\pi$ at finite $t_{\perp}$, a commensurate-incommensurate transition is expected. In this case, our numerics suggest that this transition is into a regime with a fermion density-wave-type order. This is shown in Fig. 7(b) for $W=-1.7t_{\parallel},\ t_{\perp}=t_{\parallel},\ \phi=\pi,\ L=36,\ N=12$.