Topological phases of quasi-one-dimensional fermionic atoms with a synthetic gauge field

Majorana fermions are real fermions which are equivalent to their own anti-particles. The pioneering works by Read and Green [1] and Kitaev [2] on the search for their elusive fermions opened an exciting new chapter in condensed matter physics. These original works predicted that the mysterious fermions exist as zero-energy quasi-particles bound at vortices and edges of a spinless p-wave superconductor. Subsequently, tremendous progress has succeeded in extending the platform for realizing Majorana fermions to some categories of the so-called topological superconductors [3–6]. The remarkable consequence of the self-charge conjugate property of the Majorana fermion is non-Abelian braiding statistics, where a pair of Majorana zero modes are created or annihilated by braiding their host vortices [7]. Hence, Majorana fermions possessing non-Abelian braiding statistics can provide a promising platform for fault-tolerant topological quantum computation [8]. Moreover, it has been recently unveiled that zero-energy quasi-particles exhibit multifaceted properties, not only as Majorana fermions but also as odd-frequency Cooper pair correlations [9–16].

An ideal candidate for realizing non-Abelian Majorana zero modes was a chiral p-wave superconductor with half-quantum vortices, where the low-lying quasi-particles are effectively spinless [7]. A half-quantum fluxoid has been observed in mesoscopic annular rings of Sr₂RuO₄ [17], while the half-quantum vortices in superconductors are energetically unstable against integer vortices because of the absence of a screening mechanism of spin current [18]. The other candidate in terms of a chiral p-wave superfluid is the A-phase of superfluid ³He confined to a restricted geometry with sub-micron thickness [19]. The thermodynamic stability of half-quantum vortices in such a superfluid is not trivial, because the Fermi liquid corrections which favor vortices with spin flow rather than mass flow are competitive to the strong coupling correction due to the spin fluctuation which stabilizes integer vortices without spin flow [20].

In contrast, it was demonstrated that a conventional s-wave superconductor can harbor non-Abelian Majorana zero modes [21–25]. The key finding lies in the two-dimensional (2D) Rashba-type spin–orbit interaction in background normal fermions, where the non-Abelian anyons are due to the phase twist of the spin–orbit interaction. The electron bands split by the spin–orbit interaction effectively convert the s-wave pairing to p ± ip pairing. A strong Zeeman field drives the quantum phase transition from a non-topological phase without non-Abelian anyons to a topological phase [23]. This finding provides another approach to the realization of topological quantum computation in condensed matters. Indeed, it has been proposed that topological superconductivity and Majorana fermions can be realized in a one-dimensional (1D) semiconducting wire proximity coupled with an s-wave superconductor [26–31], where semiconductors, such as InSb, have very large g-factor and strong spin–orbit interaction. The signature of Majorana fermions has recently been observed through zero bias conductance peaks in a nanowire topological superconductor [32–34] and unconventional Josephson effect in hybrid superconductor–topological insulator devices [35].

Apart from superconducting materials, cold atoms with a p-wave Feshbach resonance [36, 37] or with a synthetic gauge field [23, 24] offer an alternative playground for Majorana fermions. A spin–orbit coupling with equal Rashba and Dresselhaus strengths can be synthetically induced by applying Raman lasers to atomic gases with hyperfine spin degrees of freedom [38], whose practical scheme was first pointed out by Liu et al [39]. This scheme has recently been implemented using fermionic ⁶Li [40] and ⁴⁰K atoms [41]. In addition, schemes for creating Rashba and Dresselhaus spin–orbit coupling and three-dimensional (3D) analogue to Rashba spin–orbit coupling have theoretically been proposed [42–44]. A 1D geometry with a resonantly interacting Fermi gas can be implemented by using a 2D optical lattice, which was already utilized to search the Fulde–Ferrell–Larkin–Ovchinnikov (FFLO) state in a spin-imbalanced Fermi gas [45]. Hence, cold atoms with Raman laser-induced spin–orbit coupling provide not only a promising platform for realizing Majorana fermions [23, 46–49] but also an ideal system to study competition of various exotic superfluid phases including topological and FFLO phases. Furthermore, it has been proposed that topologically non-trivial phases, such as topological insulating phases, can be realized in a Fermi gas with a non-Abelian and Abelian gauge fields in an optical lattice [50–57].

In this paper, we study topology and quasi-particle spectra of a quasi-1D Fermi gas with a Rashba-type spin–orbit coupling. It has been demonstrated in [46–48] that a pure 1D Fermi gas with spin–orbit coupling is accompanied by exactly zero-energy states bound at the end points of atomic clouds. However, since an actual experiment is performed on a bundle of weakly coupled tubes, intertube tunneling effects are not negligible. We here demonstrate that the intertube coupling plays an important role on determining the topological properties of Fermi gases and the effective Hamiltonian is analogous to nanowire topological superconductors with multibands [58–62]. The existence of non-Abelian Majorana zero modes is ensured by introducing a new topological number associated with a mirror symmetry.

This paper is organized as follows. We begin in section 2 by introducing a tight-binding model for a bundle of 1D Fermi gases with spin–orbit coupling. In section 3, we clarify the topology of such a system, where a 1D winding number is protected by a hidden mirror symmetry. It turns out that this system provides a cold atom analogue to nanowire topological superconductors with multibands. In section 4, based on fully numerical calculations of self-consistent equations, we study the intertube tunneling effect on quasi-particle spectra. The final section is devoted to conclusions and discussions. Throughout this paper, we set ℏ = k_B = 1 and the repeated Greek indices imply the sum over x,y,z.

We here start with the Hamiltonian for spin–orbit coupled two-component fermionic atoms with an s-wave attractive interaction, g,

$$\begin{eqnarray} &&\mathcal{H} = \int \mathrm{d}{\boldsymbol{r}} {\boldsymbol{\Psi}}^{\dagger}({\boldsymbol{r}})\left[ \epsilon ({\boldsymbol{r}}) + \mathcal{S} ({\boldsymbol{r}}) \right]{\boldsymbol{\Psi}}({\boldsymbol{r}}) + g\int \mathrm{d}{\boldsymbol{r}} \psi^{\dagger}_{\uparrow}({\boldsymbol{r}})\psi^{\dagger}_{\downarrow}({\boldsymbol{r}}) \psi _{\downarrow}({\boldsymbol{r}})\psi _{\uparrow}({\boldsymbol{r}}), \end{eqnarray} \tag{ 1 }$$

where Ψ ≡ [ψ_↑,ψ_↓]^T denotes the fermionic field operators with up- and down-spins. The single-particle Hamiltonian density is defined as $\epsilon ({\boldsymbol{r}}) = -\frac {1}{2m}{\boldsymbol{\nabla}}^2 - \mu _{\mathrm { cp}} + V_{\mathrm { pot}}({\boldsymbol{r}}) - h_{\mu }\sigma _{\mu }$ with a confinement potential V_pot and σ_μ being the Pauli matrices in spin space. The Zeeman field h is naturally induced by implementing the spin–orbit coupling through the two-photon Raman process [38]. In equation (1), $\mathcal {S}$ describes the spin–orbit coupling term, which is expressed in general as

$$\begin{eqnarray} &&\mathcal{S}({\boldsymbol{r}}) = \mathrm{i}A_{\mu \nu}\sigma _{\nu}\partial _{\mu}. \end{eqnarray} \tag{ 2 }$$

Within the mean-field approximation, the Hamiltonian in equation (1) can be diagonalized in terms of the quasi-particle states. The quasi-particle states with the energy E_n are obtained by solving the so-called Bogoliubov–de Gennes (BdG) equation [63]

$$\begin{eqnarray} &&\mathcal{H}({\boldsymbol{r}}){\boldsymbol{\varphi}}_n({\boldsymbol{r}}) = E_n {\boldsymbol{\varphi}}_n({\boldsymbol{r}}), \end{eqnarray} \tag{ 3 }$$

where φ = [u_n,↑,u_n,↓,v_n,↑,v_n,↓]^T denotes the wavefunction of quasi-particles, where u_n,σ and v_n,σ describe the wavefunctions of the particle- and hole-components, respectively. The BdG Hamiltonian density is given as

$$\begin{eqnarray} \mathcal{H}({\boldsymbol{r}}) = \left( \begin{array}{@{}cc@{}} \epsilon ({\boldsymbol{r}}) + \mathcal{S} ({\boldsymbol{r}}) & \mathrm{i}\,\sigma _y \Delta ({\boldsymbol{r}}) \\ -\mathrm{i}\,\sigma _y \Delta ({\boldsymbol{r}}) & - \epsilon^{\ast}({\boldsymbol{r}}) - \mathcal{S}^{\ast}({\boldsymbol{r}}) \end{array} \right), \end{eqnarray} \tag{ 4 }$$

where Δ(r) is an s-wave pair potential with a contact interaction g, defined as Δ(r) = g〈ψ_↓(r)ψ_↑(r)〉. Note that the BdG Hamiltonian density in equation (4) holds the particle–hole symmetry, $\tau _{x}\mathcal {H}({\boldsymbol{r}})\tau _x = - \mathcal {H}^{\ast }({\boldsymbol{r}})$, where τ_μ denotes the Pauli matrices in particle–hole space.

In order to isolate Majorana zero modes from the higher energy quasi-particle states, the fermionic atoms are confined by a 2D optical lattice in the y–z plane in addition to a shallow harmonic potential along the x-direction, as shown in figure 1. The system under this confinement potential is regarded as a 2D array of N_y × N_z 1D tubes. The quasi-particle wavefunction is expanded in terms of the functions f_ℓy(y) and f_ℓz(z) localized at (x,y_ℓy,z_ℓz) as ${\boldsymbol{\varphi}} ({\boldsymbol{r}}) = \sum _{{\boldsymbol{\ell}}} {\boldsymbol{\varphi}}_{\boldsymbol{\ell}}(x)f_{\ell _y}(y)f_{\ell _z}(z)$ with ${\boldsymbol{\ell }}=(\ell _{y},\ell _{z})$, where ℓ_y = 1,...,N_y and ℓ_z = 1,...,N_z. Employing the tight-binding approximation in the y–z plane, the Hamiltonian in equation (4) reduces to

$$\begin{eqnarray} \mathcal{H}^{\rm eff}(x) = \left( \begin{array}{@{}cc@{}} \underline{\epsilon}(x) + \underline{\mathcal{S}}(x) & \mathrm{i}\,\sigma _y \underline{\Delta}(x) \\ -\mathrm{i}\,\sigma _y \underline{\Delta}^{\dagger} (x) & - \underline{\epsilon^{\ast}}(x) - \underline{\mathcal{S}}^{\ast}(x) \end{array} \right), \end{eqnarray} \tag{ 5 }$$

where the pair potential Δ(r) is transformed to a 2N_yN_z × 2N_yN_z matrix for ℓ and ℓ', $\underline {\Delta }(x)$. The single-particle Hamiltonian density $\underline {\epsilon }(x)$ is a 2N_yN_z × 2N_yN_z matrix given by

$$\begin{eqnarray} &&\fl [\underline{\epsilon} (x)]_{{\boldsymbol{\ell}},{\boldsymbol{\ell}}^{\prime}} = \left(- \frac{1}{2m}\frac{\mathrm{d}^2}{\mathrm{d}x^2} - \mu _{\rm cp} + V_{{\boldsymbol{\ell}}}(x) - h_{\mu}\sigma _{\mu} \right)\delta _{{\boldsymbol{\ell}},{\boldsymbol{\ell}}^{\prime}} - t_y (\delta _{{\boldsymbol{\ell}},{\boldsymbol{\ell}}^{\prime}+ \hat{\boldsymbol{e}}_y} +\delta _{{\boldsymbol{\ell}},{\boldsymbol{\ell}}^{\prime}- \hat{\boldsymbol{e}}_y}) - t_z (\delta _{{\boldsymbol{\ell}},{\boldsymbol{\ell}}^{\prime}+ \hat{\boldsymbol{e}}_z} + \delta _{{\boldsymbol{\ell}},{\boldsymbol{\ell}}^{\prime}- \hat{\boldsymbol{e}}_z}).\nonumber\\ \end{eqnarray} \tag{ 6 }$$

The effective potential V_ℓ(x) is given as $V_{\boldsymbol{\ell}}(x) \equiv \frac {1}{2}m\omega ^2_x x^2 + V_{\boldsymbol{\ell}}$. Here we set ${\boldsymbol{\ell}} = \ell _y \hat {\boldsymbol{e}}_y +\ell _z \hat {\boldsymbol{e}}_z$ with $\hat {\boldsymbol{e}}_y = (1, 0)$ and $\hat {\boldsymbol{e}}_z = (0,1)$. The hopping energies between intertubes are denoted by t_y and t_z. The spin–orbit coupling term also reduces to

$$\begin{eqnarray} &&\fl [\underline{\mathcal{S}} (x)]_{{\boldsymbol{\ell}},{\boldsymbol{\ell}}^{\prime}} = \mathrm{i} \sigma _{\mu}[ A_{x\mu} \partial _x \delta _{{\boldsymbol{\ell}},{\boldsymbol{\ell}}^{\prime}} + \tilde{A}_{y\mu}(\delta _{{\boldsymbol{\ell}},{\boldsymbol{\ell}}^{\prime}+ \hat{\boldsymbol{e}}_y} -\delta _{{\boldsymbol{\ell}},{\boldsymbol{\ell}}^{\prime}- \hat{\boldsymbol{e}}_y}) + \tilde{A}_{z\mu}(\delta _{{\boldsymbol{\ell}},{\boldsymbol{\ell}}^{\prime}+ \hat{\boldsymbol{e}}_z} - \tilde{A}_{z\mu}\delta _{{\boldsymbol{\ell}},{\boldsymbol{\ell}}^{\prime}- \hat{\boldsymbol{e}}_z})], \end{eqnarray} \tag{ 7 }$$

where $\skew3\tilde {A}_{\mu \nu }$ describes an effective non-Abelian gauge field, $\skew3\tilde {A}_{y\mu } \equiv \int \!f^{\ast }_{\ell _z}(z) f^{\ast }_{\ell _y+1}(y)\times A_{y\mu } f_{\ell _y}(y)f_{\ell _z}(z)\,\mathrm {d}y\,\mathrm {d}z$. Note that a quasi-1D Fermi gas has been reported in [45] with a 2D optical lattice potential. The number of tubes is typically about $N_y \times N_z \sim \mathcal {O}(10\times 10)$, so the system should be treated as a finite system.

Zoom In Zoom Out Reset image size

Under the tight-binding approximation, the resulting BdG equation reduces to an effective 1D equation along the x-axis

$$\begin{eqnarray} &&\mathcal{H}^{\rm eff}_{{\boldsymbol{\ell}},{\boldsymbol{\ell}}^{\prime}}(x) {\boldsymbol{\varphi}}_{n,{\boldsymbol{\ell}}^{\prime}}(x) = E_n {\boldsymbol{\varphi}}_{n,{\boldsymbol{\ell}}}(x), \end{eqnarray} \tag{ 8 }$$

which is numerically solved with the finite element method implemented with the discrete variable representation [70]. The BdG equation (8) is self-consistently coupled with the gap equation for the pair potential, $[\underline {\Delta } (x)]_{{\boldsymbol{\ell}},{\boldsymbol{\ell}}^{\prime }} = \Delta _{\boldsymbol{\ell}} (x) \delta _{{\boldsymbol{\ell}},{\boldsymbol{\ell}}^{\prime }}$,

$$\begin{eqnarray} &&\Delta _{\boldsymbol{\ell}}(x) = g \sum _{E_n}\left[ u_{n,\uparrow, {\boldsymbol{\ell}}}(x)v^{\ast}_{n,\downarrow, {\boldsymbol{\ell}}}(x) f(E_n) + u_{n,\downarrow, {\boldsymbol{\ell}}}(x)v^{\ast}_{n,\uparrow, {\boldsymbol{\ell}}}(x) f(-E_n) \right], \end{eqnarray} \tag{ 9 }$$

where f(E) = 1/(e^E/T  + 1) is the Fermi distribution function at a temperature T. In addition, the chemical potential μ is determined so as to preserve the total particle number

$$\begin{eqnarray} &&N = \sum _{E_n, \sigma}\int \mathrm{d}{\boldsymbol{r}}\left[ \left| u_{n,\sigma, {\boldsymbol{\ell}}}(x)\right|^2f(E_n) + \left| v_{n,\sigma, {\boldsymbol{\ell}}}(x)\right|^2f(-E_n) \right]. \end{eqnarray} \tag{ 10 }$$

The sum in equations (9) and (10) is taken over E_n∈[0,E_c], where E_c denotes the energy cutoff. The effective 1D coupling constant, g, in equation (9) is expressed in terms of an effective 1D scattering length a_1D as [64]

$$\begin{eqnarray} &&g = - \frac{2}{ma_{\rm 1D}} = - \frac{2}{m}\frac{a_{\rm 3D}/d_{\perp}}{Aa_{\rm 3D}-d_{\perp}}, \end{eqnarray} \tag{ 11 }$$

where a_3D denotes a 3D scattering length and d_⊥ is the characteristic harmonic oscillator length in y- and z-axis. The constant A is given as A ≈ 1.0326. From equation (11), the 1D scattering length a_1D is expressed in terms of a_3D as a_1D = d_⊥(A − d_⊥/a_3D). In the definition of a_1D, the sign is opposite to that of a_3D. Hence, the Cooper pairing state in a 1D Fermi gas can be stabilized by a positive a_1D which corresponds to a_3D < 0 realized in ⁶Li atoms. It is also seen from equation (11) that the pairing interaction, g, can be controlled by changing the characteristic length scale of the confinement potential in the y–z plane, d_⊥.

Before closing this section, we mention the validity of the mean-field theory in a quasi-1D Fermi gas. For a pure 1D system, in general, the quantum fluctuation plays a critical role, which violates the long-range ordering. Hence, the mean-field approximation employed in the current work might not work very well in a pure 1D system. Based on the direct comparison of the mean-field theory with the Bethe ansatz solutions, however, Liu et al [64] demonstrated that the mean-field theory provides a useful description in weakly or moderately interacting regimes of a pure 1D Fermi gas. Furthermore, the quantum fluctuation is suppressed by introducing the intertube tunneling t_y and t_z, which involves the crossover of the single-particle dispersion from 1D to 3D.

We here consider a 2D Rashba-type spin–orbit interaction

$$\begin{eqnarray} &&A_{\mu\nu}\sigma _{\nu}\partial _{\mu} = \kappa _{x} \sigma _y \partial _x - \kappa _y \sigma _x \partial _y. \end{eqnarray} \tag{ 12 }$$

To clarify the topological property of the effective Hamiltonian, we here ignore the shallow trap potential along the x-direction, i.e. ω_x → 0. Then, the effective Hamiltonian in equation (5) is rewritten with −i∂_x → k as

$$\begin{eqnarray} &&\fl \mathcal{H}^{\rm eff}_{{\boldsymbol{\ell}}, {\boldsymbol{\ell}}'}(k) = \left[ \epsilon^{(0)}_{{\boldsymbol{\ell}},{\boldsymbol{\ell}}^{\prime}}(k) - \left\{ h_x\sigma _x +h_z \sigma _z +\kappa _x k \sigma _y \right\} \delta _{{\boldsymbol{\ell}},{\boldsymbol{\ell}}^{\prime}}\right]\tau _z - \Delta _{{\boldsymbol{\ell}},{\boldsymbol{\ell}}^{\prime}}\sigma _y \tau _y\nonumber \\ &&+ \left[ h_y\sigma_y \delta _{{\boldsymbol{\ell}},{\boldsymbol{\ell}}^{\prime}} + i \tilde{\kappa}_y\sigma _x \left\{\delta _{{\boldsymbol{\ell}},{\boldsymbol{\ell}}^{\prime}+ \hat{\boldsymbol{e}}_y} -\delta _{{\boldsymbol{\ell}},{\boldsymbol{\ell}}^{\prime}- \hat{\boldsymbol{e}}_y}\right\} \right] \tau_0, \end{eqnarray} \tag{ 13 }$$

where ⁽⁰⁾_ℓ,ℓ' describes the single-particle Hamiltonian density without the Zeeman term. Here, without the loss of generality, Δ_ℓ is assumed to be real. We have also assumed h_y = 0. It is seen from equation (13) that the array of 1D tubes with a spin–orbit interaction is analogous to a semiconductor–superconductor nanowire with Nth electron bands [62], where N ≡ N_y × N_z.

We find that our system supports two different topological numbers. The first one is the 1D ${\mathbb Z}_2$ topological number for class D topological phase. Because of superfluidity, the BdG Hamiltonian (13) has the particle–hole symmetry, which allows us to define the 1D ${\mathbb Z }_2$ topological number. The 1D ${\mathbb Z}_2$ number is defined as

$$\begin{eqnarray} &&\nu=\frac{1}{\pi}\int_{-\pi}^{\pi}\mathrm{d}k A(k)+\mbox{mod 2}, \end{eqnarray} \tag{ 14 }$$

where A(k) is the geometrical phase

$$\begin{eqnarray} &&A(k)=\mathrm{i}\sum_{E_n(k)<0}\sum _{\boldsymbol{\ell}} \langle {\boldsymbol{\varphi}}_{n,{\boldsymbol{\ell}}}(k)|\partial_k {\boldsymbol{\varphi}}_{n,{\boldsymbol{\ell}}}(k)\rangle \end{eqnarray} \tag{ 15 }$$

with |φ_n,ℓ(k)〉 the Bloch wavefunction of a negative energy state of ${\cal H}^{\mathrm { eff}}(k)$. When ν is odd (even), the system is topologically non-trivial (trivial).

The second topological number comes from a remnant of a mirror reflection symmetry of the system. If one temporarily neglects the Zeeman fields h, our system is invariant under the mirror reflection to the zx-plane, as well as the time reversal. Once the Zeeman fields are applied, the mirror symmetry is lost, but a combination of the mirror reflection and the time reversal is still preserved if h_y = 0. Consequently, the Hamiltonian $\mathcal {H}_{\mathrm { eff}}(k)$ with h_y = 0 holds the following Z₂ symmetry:

$$\begin{eqnarray} &&\mathcal{T}\mathcal{M}_{zx}\mathcal{H}_{\rm eff}(k) \mathcal{M}_{zx}^{\dagger}\mathcal{T}^{-1} = \mathcal{H}^{\ast}_{\rm eff}(-k), \end{eqnarray} \tag{ 16 }$$

where $\mathcal {T} = \mathrm {i}\,\sigma _y K$ is the time-reversal operator with the complex conjugate operator K and $\mathcal {M}_{zx}=\mathrm {i}\,\sigma _yU$ the mirror reflection operator. Here U is the operator flipping the y-component of ℓ = (ℓ_y,ℓ_z),

$$\begin{eqnarray} &&U_{{\boldsymbol{\ell}}, {\boldsymbol{\ell}}'}=\delta_{\ell_y, N_y+1-\ell_y'} \,\delta_{\ell_z,\ell_z'}. \end{eqnarray} \tag{ 17 }$$

Combining the Z₂ symmetry with the particle–hole symmetry, $\mathcal {C}\mathcal {H}_{\mathrm { eff}}(k)\mathcal {C}^{-1} = - \mathcal {H}^{\ast }_{\mathrm { eff}}(-k)$, we define the chiral symmetry operator, $\Gamma = \mathcal {T}\mathcal {M}_{zx}\mathcal {C} = \tau _xU$, where $\mathcal {C} = \tau _x K$ denotes the particle–hole transformation operator with complex conjugation operator K. Then, it turns out that Γ is anti-commutable with the effective Hamiltonian

$$\begin{eqnarray} &&\left\{ \Gamma, \mathcal{H}_{\rm eff} (k)\right\} = 0. \end{eqnarray} \tag{ 18 }$$

This implies that the BdG Hamiltonian $\mathcal {H}_{\mathrm { eff}}(k)$ holds the chiral symmetry. Then, the 1D winding number is defined as [65–68]

$$\begin{eqnarray} &&w = -\frac{1}{4\pi \mathrm{i}} \int^{\infty}_{-\infty} \mathrm{d}k\, {\rm tr}\left[\Gamma{\cal H}^{-1}_{\rm eff}(k) \partial_k {\cal H}_{\rm eff}(k)\right] \end{eqnarray} \tag{ 19 }$$

which takes an integer. A similar 1D winding number was considered for 2D and pure 1D Rashba superconductors, where Γ in equation (19) is replaced by τ_x [65, 69]. The above expression (19) is a generalization of these cases into multi-tube systems.

As a consequence of the bulk-edge correspondence, these two 1D topological numbers ensure the existence of zero-energy states appearing in the end points of 1D tubes. Here we note that the parities of these two topological numbers coincide with each other

$$\begin{eqnarray} &&(-1)^{\nu}=(-1)^w, \end{eqnarray} \tag{ 20 }$$

which implies that w can be non-zero even when ν is trivial, but the opposite is not true. Therefore, the actual number of the zero-energy states is determined by w unless the Z₂ symmetry (16) is broken macroscopically. In addition, the particle–hole symmetry of superfluidity results in the Majorana property of the zero-energy state, where the creation operator γ^†_E=0 is equivalent to its own annihilation, γ_E=0 = γ^†_E=0 [66, 67]. In summary, the winding number w ensures the existence of Majorana zero modes, whereas once the Z₂ symmetry is broken by turning on h_y for example, w becomes ill-defined and the 1D ${\mathbb Z}_2$ number ν in equation (14) determines the topological stability of the Majorana zero modes.

In the case of a 2D disc geometry, where the 2D optical lattice potential is absent, the topological property is characterized by the first Chern number [24]. The non-trivial value of the Chern number ensures the existence of the gapless state localized at the circumference of the atomic cloud. In addition, the edge states carry the macroscopic mass current.

Let us start with a pure 1D system with $t_y = t_z = \tilde {\kappa }_y = 0$. Here, the BdG equation (8) coupled with the gap equation (9) is numerically solved with the set of parameters: T = 0, E_c = 4E_F and κ_x = 1. The Zeeman field is applied along the $\hat {\boldsymbol{z}}$-axis, which does not break the Z₂ symmetry in equation (16): h = (0,0,h). In realistic situations [38, 49], the strength of the spin–orbit coupling, κ_x, depends on the wavelength of applied lasers. Throughout this work, we fix the pairing interaction, γ = 1.4, where $\gamma = \frac {1}{\pi \sqrt {N}}(\frac {d}{a_{\mathrm { 1D}}})$ is the dimensionless coupling constant [64], with $d = \sqrt {1/(m\omega _x)}$ being the harmonic oscillator length. The total particle number is fixed to be N = 200 in each tube, where the Fermi energy per one tube is E_F = 100ω_x.

Before going to numerical results, we first clarify the effect of the mean-field potential associated with the on-site interaction g, which is dropped in the BdG Hamiltonian density in equation (4). The on-site interaction term changes the single-particle Hamiltonian density in equation (4) to (r) → (r) + gρ(r) + gm(r)σ_z, where ρ(r) and m(r) denote the local particle and spin densities. Within the local density approximation, the maximum value of the local density ρ(x) in an ideal Fermi gas confined to a 1D harmonic oscillator is estimated as $\max \rho (x) \equiv \rho _0 = 2\sqrt {N}/\pi d$. The dimensionless parameter γ gives a rough estimation about the ratio of the local potential and the Fermi energy, γ ∼ gρ₀/E_F, where in our calculated system, μ ∼ E_F. Since the potential term changes the local chemical potential and the local Zeeman field, it quantitatively alters quasi-particles with finite energies and the critical Zeeman field above which Majorana zero modes appear. However, the on-site interaction term does not affect the topological properties associated with Majorana zero modes, where the mirror symmetry is preserved. Therefore, we here ignore the effect of the on-site interaction term in the BdG equation.

For the pure 1D system, topology of the BdG Hamiltonian for each tube is characterized by the winding number w in equation (19) with N_y = N_z = 1. With spatially constant Δ and μ, w is given by

$$\begin{eqnarray} &&w=\frac{1}{2\pi\, \mathrm{i}}\int_{-\infty}^{\infty}\mathrm{d}k \partial_k \ln\left[\det {\cal Q}(k)\right], \end{eqnarray} \tag{ 21 }$$

where $\mathcal {Q}(k) \equiv \epsilon ^{(0)}(k) - ( h_x\sigma _x + h_z\sigma _z + \kappa k \sigma _y)+ \mathrm {i}\,\Delta \sigma _y$. Then, it is found that w = 1 when $|h| > h_{\mathrm { c}} \equiv \sqrt {\mu ^2 + \Delta ^2}$. Therefore, the lower magnetic field regime is the topologically trivial phase and the critical field at h_c involves the topological phase transition. For our system, however, a more careful consideration is needed. Since fermionic atoms are confined by a trap potential, a spatially inhomogeneous superfluid is realized naturally. In addition, the inhomogeneous pair potential Δ(x), which is self-consistently determined by the gap equation and the BdG equation, depends on the Zeeman fields significantly. In contrast to semiconductor–superconductor junction systems, these two characteristics cannot be neglected. This means that within the local density approximation, the critical field h_c and 1D winding number w should be replaced by $h_{\mathrm { c}}(x) \equiv \sqrt {\mu ^2(x) + \Delta ^2(x)}$ and w(x), where μ(x) is the local chemical potential including the confinement potential along the x-axis. The inhomogeneity and self-consistency of Δ(x) and μ(x) play a critical role on the topological property.

Figure 2(a) shows the spatial profile of Δ(x) at h = 0.20E_F and 0.36E_F, where the harmonic trap potential along the x-axis, $\frac {1}{2}m\omega ^2_x x^2$, is taken into account. In the case of h = 0.36E_F, the intermediate region between 8.8d ≲ |x| ≲ 14.6d, where the pair potential has a dip, becomes topologically non-trivial, i.e. it satisfies h > h_c(x) and w(x) ≠ 0, while the inner region within |x| ≲ 8.8d is not. It is also found that the outermost region is not topological again, because the local chemical potential, $\mu (x) \equiv \mu -\frac {1}{2}m\omega ^2_xx^2$, changes its sign at the Thomas–Fermi edge $x = x_{\mathrm { TF}} \equiv d\sqrt {2\mu /\omega _x} \sim 14 d$ and its magnitude becomes large so as h < h_c(x).

Zoom In Zoom Out Reset image size

In the absence of the hopping between tubes, four zero modes appear in each tube. At h = 0.36E_F, |w| = 1 can be realized in the region within 8.8d ≲ |x| ≲ 14.6d, otherwise w = 0. As shown in figure 2(b), the lowest and second lowest energy states are tightly bound at the phase boundaries at x ∼ ± 8d and ±15d. We summarize the field dependence of the energy spectrum E_n in figure 3(a). In our calculated system, the level spacing due to the harmonic potential along the x-axis is given as 0.01E_F. Thus, in figure 3(a), the quasi-particles with E > 0.01E_F can be regarded as the 'continuum states' and those with E < 0.01E_F are referred to as the bound states. The quasi-particle states having E ≪ 0.01E_F can be referred to as the 'zero-energy' states. It is seen from figure 3(a) that the low-lying eigenenergies go to zero as h increases, because of the interference between two zero modes localized at x/d ∼ ± 8 and ±15 [70, 71]. In the Zeeman field regime higher than 0.35E_F, the zero-energy states are split into two branches. It turns out that the upper (lower) branches correspond to the quasi-particle states bound at the inner (outer) edges at x ∼ ± 8d (± 15d), as displayed in figure 2(b). As h further increases, the amplitude of the pair potential Δ(x) decreases. This implies that the wavefunction of the zero-energy states spreads, giving rise to the hybridization of zero-energy states.

Zoom In Zoom Out Reset image size

Now, let us clarify how intertube tunneling affects the low-lying quasi-particle spectra. The bundle of 1D tubes is coupled with each other through the hoppings t_y and t_z and the spin–orbit interactions $\tilde {\kappa }_y$. In this situation, the topological winding number w is not defined for each tube, but is defined only for a whole system of tubes. As a result, some of the Majorana zero modes become non-zero modes, as shown later. Note that if the hopping t_z is small enough, it does not change the winding number and the topological property because of the two dimensionality of the Rashba-type spin–orbit coupling. This implies that zero-energy states are dispersionless and insensitive to N_z. Figure 3(b) shows the field dependence of low-lying quasi-particle energies for N_y = 2 and 3, where we fix the parameters, t_y = 0.01E_F and $\tilde {\kappa }_y/\kappa _x = 0.5$. In the case of N_y = 2, the zero-energy states are lifted to finite energies by their hybridization through the intertube tunneling t_y ≠ 0. In contrast, as shown in figure 3(b), the zero-energy state in the case of N_y = 3 survives even in the presence of the finite intertube tunneling. Figure 3(c) summarizes the low-lying spectra as a function of N_y, where h is fixed to be h/E_F = 0.36 which corresponds to the topological phase in the case of N_y = 1. It is clear that the zero-energy states exist when N_y is odd, while in the case of even N_y's, the intertube tunneling lifts the low-lying states from the zero energy.

To understand the low-lying energy states in the case of odd N_y's, which remain zero modes, we show in figure 4 the local density of states (LDOS) at the site ℓ defined as

$$\begin{eqnarray} \mathcal{N}_{\sigma, {\boldsymbol{\ell}}}(x,E) &=& \sum _{E_n> 0} \left[ \left| u_{n,\sigma,{\boldsymbol{\ell}}}(x)\right|^2\delta(E-E_n) + \left| v_{n,\sigma,{\boldsymbol{\ell}}}(x)\right|^2\delta(E+E_n) \right]. \end{eqnarray} \tag{ 22 }$$

The LDOS for up-spins (down-spins) in the case of N_y = 3 is displayed in figures 4(a) and (c) ((b) and (d)). It is clearly seen from figures 4(a) and (b) that the LDOS for up-spins at ℓ_y = 1 and 3 is accompanied by the zero-energy states which are bound at the end points of the tube (x/d ≈ ± 15), while $\mathcal {N}_{\downarrow }$ has sharp peaks at x/d ≈ ± 8 corresponding to the inner phase boundaries between w(x) = 0 (non-topological region) and 1 (topological region) in the case of N_y = 1. In contrast, the intertube tunneling splits the zero-energy states into the positive and the negative energy states in the LDOS at ℓ_y = 2, which have a mini-gap with ±0.02E_F. In figure 5, we also show the amplitude of the zero modes for N_y = 3, 5 and 7, respectively. In all cases, the wavefunctions have large amplitudes at the tubes located at odd ℓ_y's. At even ℓ_y's, their amplitudes are almost negligible.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

These intertube tunneling effects can be understood as follows. As was shown above, when one neglects the intertube couplings t_y, t_z and $\tilde {\kappa }_y$, each tube supports four Majorana zero modes localized at x ∼ ± 8d and ±15d. Now let us denote one of them (say, the zero mode localized at x ∼ ± 15d) as γ_ℓy (ℓ_y = 1,...,N_y), and consider how the intertube couplings affect them. When t_y and $\tilde {\kappa }_y$ are turned on, the zero modes on neighboring tubes are coupled by the intertube tunneling

$$\begin{eqnarray} &&{\cal H}=\mathrm{i}\,t\gamma_1\gamma_2+\mathrm{i}\,t\gamma_2\gamma_3+\cdots+\mathrm{i}\,t\gamma_{N_y-1}\gamma_{N_y}, \end{eqnarray} \tag{ 23 }$$

where t denotes the induced tunneling coupling. Note that t is real since γ_ℓy is a Majorana zero mode satisfying γ_ℓy = γ^†_ℓy. Equation (23) is rewritten as ${\cal H }=\Gamma ^{\dagger }\hat {H}\Gamma /2$ with

$$\begin{eqnarray} \hat{H}= \left( \begin{array}{@{}ccccc@{}} 0 & \mathrm{i}\,t & 0 & \cdots & 0\\ -\mathrm{i}\,t & 0 & \mathrm{i}\,t & \ddots & \vdots \\ 0 & -\mathrm{i}\,t & 0 & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & \mathrm{i}\,t \\ 0 & \cdots & 0 & -\mathrm{i}\,t & 0 \end{array} \right) \end{eqnarray} \tag{ 24 }$$

and Γ = (γ₁,γ₂,...,γ_Ny − 1,γ_Ny)^t. Diagonalizing the N_y × N_y matrix $\hat {H}$, one can examine the effects of the intertube tunneling.

It can be easily shown that $\hat {H}$ has a single zero eigenvalue for odd N_y's, while it does not for even N_y's. This result naturally explains why Majorana zero modes survive only for odd N_y's. One also finds that the zero eigenstate of $\hat {H}$ has the following form:

$$\begin{eqnarray} &&(1,0, -1)^{\mathrm{t}} \quad \mbox{for $N_y=3$}, \nonumber\\ &&(1,0, -1,0, 1)^{\mathrm{t}} \quad \mbox{for $N_y=5$},\\ &&(1,0, -1,0, 1, 0, -1)^{\mathrm{t}} \quad \mbox{for $N_y=7$},\nonumber \end{eqnarray} \tag{ 25 }$$

which explains qualitatively why the remaining Majorana zero modes illustrated in figure 5 have large amplitudes on tubes at odd ℓ_y's.

The robustness of the zero-energy states against the intertube tunneling is also understood by the topological number w. As we mentioned above, even in the presence of intertube tunneling, the winding number w is well defined for a whole system of tubes. Since one can turn off the intertube tunneling without the bulk gap closing, the value of w can be evaluated by setting $t_y=t_z=\tilde {\kappa }_y=0$ in equation (19). Then one obtains

$$\begin{eqnarray} &&w={\rm tr}\,U\frac{1}{2\pi \mathrm{i}}\int_{-\infty}^{\infty}\mathrm{d}k \partial_k \ln \left[\det{\cal Q}(k)\right]. \end{eqnarray} \tag{ 26 }$$

Noting

$$\begin{eqnarray} {\rm tr}\,U=\left\{ \begin{array}{@{}cc@{}} 0 & \mbox{for even $N_y$'s},\\ N_z & \mbox{for odd $N_y$'s}, \end{array} \right. \end{eqnarray} \tag{ 27 }$$

one can valuate w as

$$\begin{eqnarray} |w|=\left\{\begin{array}{@{}ll@{}} N_z & \mbox{for odd $N_y$'s},\\ 0 & \mbox{for even $N_y$'s} \end{array} \right. \end{eqnarray} \tag{ 28 }$$

when $\sqrt {h_x^2+h_z^2}>h_{\mathrm { c}}$. This implies that Majorana zero modes survive for odd N_y's.

In this paper, we have investigated the effect of intertube tunneling in a quasi-1D Fermi gas with a Rashba-type spin–orbit coupling. From the argument based on the symmetry of the effective Hamiltonian, the topological property has been studied. In the absence of the intertube tunneling, the 1D winding number in equation (19) ensures the existence of Majorana zero modes bound at the end points of each tube. It also ensures the existence of the Majorana zero modes in the presence of the intertube tunneling if the number of tubes in the y-direction is odd. Using full numerical calculations of self-consistent equations, we have confirmed that this topological property is clearly reflected in low-lying quasi-particle states. These behaviors of low-lying quasi-particles associated with the Z₂-symmetry protected topology might be detectable through momentum-resolved radio-frequency spectroscopy [72, 73].

Here, we have considered the 2D Rashba spin–orbit coupling which has not been realized in atomic gases yet. It should be mentioned that the results obtained in this work are not straightforwardly applicable to Fermi gases under the realistic spin–orbit coupling with equal Rashba and Dresselhaus strengths [38, 40, 41]. Because of an additional symmetry specific to the equal Rashba and Dresselhaus spin–orbit coupling, the later situation is accompanied by zero-energy states regardless of even–odd parity of N_y. It is also important to mention that topological superfluidity protected by mirror symmetry can be affected by the orientation of the applied Zeeman field, because the mirror symmetry is explicitly broken. The topological property in a system with the breaking of the mirror symmetry can be associated with the 1D $\mathbb {Z}_2$ number defined in equation (14). The details will be reported elsewhere [74].

Finally, we would like to mention a generalization of the present consideration to semiconductor–superconductor nanowire systems. The 1D winding number (19) introduced in this paper is also applicable to semiconductor–superconductor nanowire with multichannels. If we consider the nanowire extending in the x-direction on top of an s-wave superconductor in the xy-plane, the system is naturally supposed to be invariant under the mirror reflection, y → − y, to the xz-plane. This mirror symmetry could be broken under Zeeman fields, but the Z₂ symmetry (16) remains if the Zeeman fields are applied in the x- or z-direction. Then, the topological number w in equation (19) is defined in the same manner. From arguments similar to the above, one finds that w is non-zero if the Zeeman field h satisfies |h| > h_c and the number of channels in the y-direction of the nanowire is odd. Indeed, under this condition, |w| is equal to the number of channels in the z-direction of the nanowire. Note that, in contrast to the 1D ${\mathbb Z}_2$ number in equation (14), w can be non-zero even when the total number of channels in the nanowire is even, since it is given by the sum of the channels in the z- and y-direction. In addition to the fermionic gas system studied in this paper, the local density operator of the Majorana zero modes vanishes [67]. This implies that the coupling between the Majorana zero modes and non-magnetic local disorder potential also vanishes, and thus the Majorana zero modes are stable against weak non-magnetic disorders.

The authors are grateful to Takuto Kawakami for fruitful discussions and comments. This work was supported by a Grant-in-Aid for Scientific Research from MEXT of Japan, 'Topological Quantum Phenomena' no. 22103005 and JPSJ no. 22540383.