Majorana corner pairs in a two-dimensional $s$-wave cold atomic superfluid

Majorana zero modes (MZMs) have attracted great attention in past decades owing to their non-Abelian exchange statistics and potential applications as topologically protected qubits Kitaev2003 ; Nayak2008 . They also exhibit significant physics in a range of disciplines such as nuclear and particle physics Elliott2015 . Strenuous efforts to search for MZMs are underway in both theories and experiments. In recent years, a variety of schemes to realize Majorana excitations have been proposed Qix2011 ; Hansan2010 ; zhang2008 ; Sato2009 ; Jiang2011 ; Alicea2011 ; Alicea2012 ; Read2000 ; Phong2017 ; Oppen2010 ; Lutchyn2010 ; Jiang2016 ; Sau2010 ; Tewari2012 ; Teemu013 by utilizing $p$-wave superconductors (SCs) or superfluid (SFs) Read2000 ; Phong2017 , or SCs and SFs with effective $p$-wave paring via spin-orbit coupling (SOC) and $s$-wave pairing Hansan2010 ; Qix2011 ; zhang2008 ; Sato2009 ; Jiang2011 . Remarkable experimental progresses have been made in condensed-matter systems Mourik2012 ; Finck2013 ; Nadj2014 ; Xuj2015 ; Wang2016 ; Heq2017 . The experimental realization of SOC in ultracold atomic gases offer another clean platform to explore Majorana physics Lin2011 ; Zhangj2012 ; Wang2012 ; Williams2013 ; Huang2016 ; Meng2016 ; Wuz2016 . In these platforms, the interplay among SOC, Zeeman fields and $s$-wave interactions could produce non-Abelian topological superfluids (TSFs) that host Majorana excitations. There have been several tantalizing proposals for realizing and tuning Majorana excitations, for example, by creating topological defects (such as SF vortices or lattice dislocations) or defect chain Lang2014 ; Deng2015 .

The emergence of Majorana excitations can be intuitively understood by the low-energy theory. A pair of MZMs exist at the kinks where the pairing potential or SOC changes the sign (which corresponds to the sign change of Dirac mass or velocity in Jackiw–Rebbi model). In solid-state materials, the SOC kinks are difficult to tune, while the kinks of pairing potentials can be realized through Josephson junctions in superconducting nanowires, as well as the corners and hinges in recently proposed higher-order topological SCs (TSCs) J. Langbehn2017 ; Benalcazar2017 ; Song2017 ; Fang2017 ; Schindler2018 ; Khalaf2017 ; ZhuX2018 ; Khalaf2018 ; Yan2018 ; Wang2018 ; wangyu2018 ; Liu2018 ; Hsu2018 . In particular, for two-dimensional (2D) second-order TSCs, the bulk topology of the 2D system offers 1D edge modes, which have different topologies for adjacent edges due to the change of the pairing sign, leading to zero-energy Majorana Kramers pairs or Majorana modes at corners Yan2018 ; Wang2018 ; wangyu2018 ; Liu2018 . On the other hand, in cold atomic system, the kinks of pairing potentials may be obtained by soliton excitations Xu2014 ; LiuX2015 . However, these cold atomic systems suffer from dynamical instability in the presence of perturbations. Thus, proposals for realizing robust MZMs in atomic systems, through other manners like SOC kinks, are highly in demand.

In this paper, we propose feasible schemes to realize SOC kinks using trapped ultracold fermionic atoms on a 2D optical lattice, and show that our system supports Majorana pairs in a vortex-free configuration. The main results are listed below:

(i) Effective 1D modes on a rectangular geometry would emerge in the 2D system through engineering on-site potential of the rectangle. In the presence of 1D equal Rashba-Dreeshauls (ERD) SOC, Majorana corner pair emerges at the corner of the rectangle with a proper $s$-wave pairing.

(ii) Each edge of the rectangular defect is characterized by a 1D topological SF in the BDI class, and Majorana pair exists at the interface of two adjacent edges which have different signs of SOC (clockwise or anticlockwise along the defective geometry). Our system is in analog with the higher-order TSCs except that the low-energy 1D model is induced by the defect rectangle rather than the bulk topology, and the MZMs are induced by SOC kinks rather than pairing kinks between two adjacent boundaries Yan2018 ; Wang2018 . Our system is more concise and experimentally friendly since no tricky unconventional pairings like $d$-wave or $s_{\pm}$-wave are required.

(iii) Our system can be realized with currently already established experimental techniques in cold atoms, including 1D ERD SOC by Raman lasers Lin2011 ; Zhangj2012 , single-site addressing in 2D optical lattices Peter2009 ; Bakr2009 ; Bakr2010 ; Sherson2010 ; Weitenberg2011 , and tunable $s$-wave interaction through Feshbach resonance Thorsten2006 ; Chin2010 .

(iv) The Majorana corner pair is robust against shape deformation of the defective rectangle, even to the extent of a defective loop. In addition, the SOC direction dictates on which corners the Majorana pair resides. Therefore the incident direction of Raman lasers can be used to manipulate the Majorana pairs.

The paper is organized as follows. In Sec. II we first introduce the Hamiltonian with 1D ERD SOC on optical lattices, and obtain the phase diagram consisting of metal and $s$-wave SF phases. In Sec. III, we study the case with a defective rectangle, and find Majorana pairs emerge at the corners. We extend the discussion to a ring-shaped geometry in Sec. IV, where the Majorana pairs arise naturally due to soft domain walls of SOC. Finally, we make conclusions and discussions in Sec. V.

We utilize atomic hyperfine states as the pseudospin states $\left|\uparrow\right\rangle$ and $\left|\downarrow\right\rangle$, as illustrated in Fig. 1. The SOC is synthesized by two counter-propagating Raman lasers coupling the two hyperfine states. The single atom motion in 2D real space is described by the Hamiltonian $\hat{H}_{0,a}=\frac{\vec{k}^{2}}{2m_{0}}+\frac{\delta}{2}\sigma_{z}+(\Omega e^{2i\vec{k}_{0}\cdot\vec{r}}\left|\downarrow\right\rangle\left\langle\uparrow\right|+h.c.)$, where the reduced Planck constant $\hbar$ has been set to be $1$, $\Omega\propto\Omega_{1}\Omega_{2}^{\ast}$ is the strength of Raman coupling, and $\vec{k}_{0}=k_{0,x}\vec{e}_{x}+k_{0,y}\vec{e}_{y}$ is the wave vector of the Raman laser. The off-diagonal terms correspond to a spin flip process accompanied by a momentum transfer of $2\vec{k}_{0}$, describing the SOC effects. The detuning term reads $\delta=\omega_{z}-\delta\omega$, where $\omega_{z}>0$ is the energy difference between these two hyperfine states, and $\delta\omega$ denotes the frequency difference between two Raman laser beams. In the following, we assume that other hyperfine levels are far off-resonance under the two-phonon process, for example, by quadratic Zeeman shift. The Hamiltonian is first transformed by a unitary matrix, namely, $\hat{H}_{0,b}=U\hat{H}_{0,a}U^{-1}$, where $U=$diag$(e^{-i\vec{k}_{0}\cdot\vec{r}},e^{i\vec{k}_{0}\cdot\vec{r}})$. We then perform another pseudo-spin rotations $\tilde{U}=e^{-i\frac{\pi}{4}\sigma_{z}}e^{-i\frac{\pi}{4}\sigma_{y}}$ to obtain $\hat{H}_{0}=\tilde{U}\hat{H}_{0,b}\tilde{U}^{\dagger}$, which can be written as

The Raman transition produces a desired ERD SOC. In the following, we assume $\delta=0$ for convenience.

1D system suffers from strong quantum fluctuations, which could eliminate $s$-wave SF order, together with the Majorana modes. This motivates us to investigate the physics in 2D, where quasi-long-range SF order exists below the Berezinskii-Kosterlitz-Thouless (BKT) transition temperature. We concentrate on the lowest (nearly) degenerate bands for constructing a tight-binding model. From Eq. (1), through the operator $\hat{\Psi}\left(\vec{r}\right)=\sum_{i,\sigma}\hat{c}_{i,\sigma}\psi_{\sigma}\left(\vec{r}-\vec{r}_{i}\right)$ with $\psi_{\sigma}\left(\vec{r}-\vec{r}_{i}\right)$ the Wannier function at site $i$, we obtain a second-quantization formula:

where $t_{x}=t_{0}\cos\left(k_{0}\cos\theta\right)$, $t_{y}=t_{0}\cos\left(k_{0}\sin\theta\right)$, $t_{sox}=t_{0}\sin\left(k_{0}\cos\theta\right)$, $t_{soy}=t_{0}\sin\left(k_{0}\sin\theta\right)$, and the bare hopping strength reads $t_{0}=-\int d\vec{r}\psi_{\sigma}^{\ast}\left(\vec{r}-\vec{r}_{i}\right)\left(\frac{k_{x}^{2}+k_{y}^{2}}{2m_{0}}+V_{\mathrm{lat}}\right)\psi_{\sigma}\left(\vec{r}-\vec{r}_{i+1}\right)$. Here, we have chosen the basis $\hat{c}_{i}=\left(\hat{c}_{i,\uparrow},\hat{c}_{i,\downarrow}\right)^{T}$ and denoted $h_{z}=\Omega$. The incident angle $\theta$ is illustrated in Fig. 1. The lattice spacing is set to be $a=1$. Hereafter, we set $t_{0}=1/\cos\left[k_{0}\sin\left(\pi/4\right)\right]=1/\cos(\sqrt{2}\pi/4)$ for convenience. When $\theta=\pi/4$, we have approximately $t_{sox}=t_{soy}\approx 2t_{x}=2t_{y}$.

We consider an attractive SU$(2)$-invariant interaction $\hat{H}_{int}=-\sum_{i}U\hat{n}_{i,\uparrow}\hat{n}_{i,\downarrow}$ and study the superfluid phase under mean-field approach by solving the $s$-wave superfluid order parameter $\Delta_{s}=U\left\langle c_{i,\downarrow}c_{i,\uparrow}\right\rangle$ self-consistently. The Bogoliubov de Genns (BdG) Hamiltonian in the Nambu basis $\Psi_{k}=\left(c_{k\uparrow},c_{k\downarrow},c_{-k\downarrow}^{\dagger},-c_{-k\uparrow}^{\dagger}\right)^{T}$ is described by $\hat{H}_{s}=\sum_{k}\Psi_{k}^{\dagger}H\left(k\right)\Psi_{k}$ with

where $\epsilon_{k}=-2\left(t_{x}\cos k_{x}+t_{y}\cos k_{y}\right)-\mu$, $\gamma_{k}=2\left(t_{sox}\sin k_{x}+t_{soy}\sin k_{y}\right)$, $\sigma$ and $\tau$ are Pauli matrices acting on the spin and particle-hole spaces, respectively. By minimizing free energy with respect to the order parameter $\Delta_{s}$ and chemical potential $\mu$, we may derive the following self-consistent equations

where $n_{f}$ is the particle filling factor, $N_{l}$ is the number of lattice sites, and other parameters are defined as $\beta=1/\left(k_{B}T\right)$ with $k_{B}$ the Boltzmann constant and $T$ the temperature, $\xi_{k,\nu=\pm}=\sqrt{\epsilon_{k}^{2}+\Delta_{s}^{2}+m_{k}^{2}+2\nu g_{k}}$, $g_{k}=\sqrt{\epsilon_{k}^{2}m_{k}^{2}+\Delta_{s}^{2}h_{z}^{2}}$, and $m_{k}=\sqrt{h_{z}^{2}+\left|\gamma_{k}\right|^{2}}$. By numerically solving equations (4) and (5), we obtain phase diagrams at zero temperature for pairing order $\Delta_{s}$ and quasiparticle energy gap $E_{g}$ in Fig. 2 (a) and (b), respectively. Fig. 2 (a) confirms the phase transition from a metal (M) phase to an $s$-wave SF. From panel (b), we find a finite gap for Bogoliubov quasiparticle excitations in proper parameter region in the SF phase. The energy gap also survives on a finite-size sample and could protect Majorana modes from lower extended states.

Given a proper local dip potential, the defect chain enjoys a non-trivial topology, belonging to the BDI class. It can be characterized through a winding number, which is discussed in Appendix A. Similarly, with a local dip $\mu_{d}$, we can get a defect rectangle in the 2D optical lattice as illustrated in Fig. 1, where SOC domain walls (anticlockwise or clockwise) naturally arise at two corners. In the following, we will first focus on the continuum limit to explore the nature of the emerged Majorana pairs, supplemented with self-consistent numerical calculations on a 2D optical lattice.

We assume that with appropriate $\mu_{d}$, the defect rectangle enters the TSF phase while the rest part remains trivial. As a result, we could assume that the topological defect rectangle is isolated from the 2D bulk. The numerics performed on a 2D optical lattice with an imprinted defect-rectangle also supports this assumption later. From Eq. (3), the low-energy Hamiltonian expands around $\vec{k}=\left(0,0\right)$ on edges $\mathrm{m}=$I, II, III, IV (see Fig. 1) and is then given by

where $t_{\mathrm{I}}=t_{\mathrm{III}}=t_{y}$, $t_{\mathrm{II}}=t_{\mathrm{IV}}=t_{x}$, $k_{\mathrm{I}}=k_{\mathrm{III}}=k_{y}$, $k_{\mathrm{II}}=k_{\mathrm{IV}}=k_{x}$, $t_{s,\mathrm{I}}=t_{s,\mathrm{III}}=2t_{soy}$ and $t_{s,\mathrm{II}}=t_{s,\mathrm{IV}}=2t_{sox}$. The on-site chemical potential is $\mu_{\mathrm{m}}=\mu+2\left(t_{x}+t_{y}\right)-\mu_{d}$ with $\mu_{d}$ the dip potential, and $\Delta_{s,\mathrm{m}}$ the $s$-wave pairing on each edge. Without loss of generality, we set incident angle of Raman lasers $\theta=\frac{\pi}{4}$ such that $t_{x}=t_{y}=t$, $t_{sox}=t_{soy}=t_{\mathrm{so}}$, $\mu_{\mathrm{m}}=\mu_{\mathrm{edge}}=\mu+4t-\mu_{d}$, and assume the $s$-wave SF order parameter is nearly uniform on the four edges $\Delta_{s,\mathrm{m}}=\Delta_{\mathrm{edge}}$. For later convenience, we take an “edge coordinate” $s$, in which we take the anticlockwise direction as positive. In such a coordinate, the low-energy edge Hamiltonian reads

with $\alpha\left(s\right)=-2t_{so}$, $2t_{so}$, $2t_{so}$, $-2t_{so}$ for edge I-IV respectively. Remarkably, while the terms $\Delta_{\mathrm{edge}}$ and $\mu_{\mathrm{edge}}$ remain the same on the four edges, the effective coupling $\alpha\left(s\right)$ changes sign at two of four corners (the corner between the edges I (III) and II (IV)), forming two SOC domain walls as illustrated in Fig. 1. This will give rise to a Majorana pair if $h_{z}^{2}>\mu_{\mathrm{edge}}^{2}+\Delta_{\mathrm{edge}}^{2}$. Specifically, at the corner between edge I and II (corner $s=0$ in our coordinate), two orthogonal wave functions for MCMs are given by

Here, $C_{\pm}$ are normalization constants and $e^{i\phi_{\pm}}=\frac{h_{z}\left[i\left(\alpha\eta_{\pm}\mp\Delta_{\mathrm{edge}}\right)-\left(t\eta_{\pm}^{2}+\mu_{\mathrm{edge}}\right)\right]}{\left(t\eta_{\pm}^{2}+\mu_{\mathrm{edge}}\right)^{2}+\left(\alpha\eta_{\pm}\mp\Delta_{\mathrm{edge}}\right)^{2}}$, where $\eta_{\pm}=\frac{1}{2}\sqrt{-\frac{2\varkappa}{3}+\delta}\mp\frac{1}{2}\sqrt{-\frac{4\varkappa}{3a}-\delta+\zeta_{\pm}}>0$, $\zeta_{\pm}=\mp 2d/(a\sqrt{-2c/(3a)+\delta}$, $\delta=\frac{\sqrt[3]{2}\delta_{1}}{3a\sqrt[3]{\delta_{2}+\sqrt{-4\delta_{1}^{3}+\delta_{2}^{2}}}}+\frac{\sqrt[3]{\delta_{2}+\sqrt{-4\delta_{1}^{3}+\delta_{2}^{2}}}}{3\sqrt[3]{2a}}$, $\delta_{1}=\varkappa^{2}+12ae$, $\delta_{2}=2\varkappa^{3}+27ad^{2}-72a\varkappa e$, $a=t^{2}$, $\varkappa=\alpha^{2}+2t\mu_{\mathrm{edge}}$, $d=-2\alpha\Delta_{\mathrm{edge}}$ and $e=\Delta_{\mathrm{edge}}^{2}-h_{z}^{2}+\mu_{\mathrm{edge}}^{2}$. The vectors $\left|y_{\pm}\right\rangle_{\sigma}$ and $\left|y_{\pm}\right\rangle_{\tau}$ are eigenstates of operators $\sigma_{y}$ and $\tau_{y}$, respectively. Following similar approach, we could also find two Majorana modes at the corner between edges III and IV (see Appendix B for details). We emphasize that as long as the four edges are in the TSF phase, the very existence of Majorana pairs is robust against the fluctuations of chemical potential and SF order parameter.

With the above understanding of continuum systems, we now proceed to study the discrete cases on an optical lattice shown in Fig. 1. The total Hamiltonian now becomes

where $i\in\square$ enumerate each site with the dip potential (a rectangular geometry in this case). The local $s$-wave superfluid order parameter in real space is determined in a self-consistent manner Jiang2016 , as well as the quasiparticle energy spectra and wave functions. On the defect rectangle, the system is topological once $h_{z}>\sqrt{\tilde{\mu}^{2}+\Delta_{s}^{2}}$ and $\tilde{\mu}=\mu+2t_{x}+2t_{y}-\mu_{d}$. In our self-consistent numerical calculations, we take the lattice sizes $n_{x}=n_{y}=30$, and the defect rectangle is given by $n_{x}^{d}=n_{y}^{d}=22$. The SF order parameter $\Delta_{s,i}$ is shown in Fig. 3 (a), which has a constant phase across the entire system. Fig. 3 (b) shows the quasiparticle energy spectrum, where four Majorana bound states (two Majorana corner pairs) exist in the energy gap. A small energy splitting is observed as a result of finite-size effect. Fig. 3 (c) shows the density distribution of the bounded Majorana corner states, which clearly demonstrates its localization at the corners of the defect rectangle. The Majorana corner pairs are robust against the perturbations of chemical potential and SF order parameter that preserve chiral symmetry. We have confirmed this point by numerical calculations.

The incident angle of Raman lasers can change the SOC and the nearest-neighbor hopping, and thus alter the Majorana bound states. Figs. 4 (a) and (b) illustrate the corresponding phase diagram with respect to $\mu$-$\theta$ and $\mu_{d}$-$\theta$, where Majorana corner pairs exist in the topological region (T). It is found that the Majorana pairs is also robust to certain variation of the incident angle $\theta$. We remark that if the sign of $\theta$ is reversed, the Majorana pairs appear at another two corners (the interfaces of II-III and I-IV) as shown in Fig. 3(d), which can be compared with Fig. 3(c). Thus, our proposed setup provides better tunability for manipulating Majorana bound states.

In this section, we study the case with a ring-shaped defect line. Here, the optical lattice is removed and we focus on the low-energy effective 1D model for simplicity. The effective model is illustrated in Fig. 5 (a) and we find that soft domain walls of SOC naturally arise on the ring, which leads to the emergence of Majorana corner pairs.

Without loss of generality, we assume the momentum kick by Raman lasers is along the $x$ direction. Under a spin-rotation $\sigma_{y}\rightarrow$ $\sigma_{x}$ with $k_{0,y}=0$ in Eq. (1), the SOC has the form $\alpha_{0}k_{x}\sigma_{x}$, where the coupling constant is given by the ratio of laser wavevector and atomic mass, i.e., $\alpha_{0}=k_{0}/m_{0}$. Hence, in the continuum limit, the effective Hamiltonian reads

where $\Delta_{s}$ is an $s$-wave SF order. For simplicity, we set $\Delta_{s}$ to be real. The relation $h_{z}>\sqrt{\mu^{2}+\Delta_{s}^{2}}$ holds in the topological regions.

In a polar coordinate $(\rho,\phi)$, the above Hamiltonian becomes a function of polar angle $\phi$ on a ring with given radii $\rho$, i.e.,

where $\eta=\frac{1}{2m_{0}\rho^{2}}$, $\tilde{\alpha}_{0}=\frac{\alpha_{0}}{\rho}$ and $\mu^{\prime}=\mu-\frac{k_{0}^{2}}{2m_{0}}$ (more details are discussed in Appendix 7). The Hamiltonian $\mathcal{H}\left(\phi\right)$ has particle-hole symmetry $\mathcal{PH}\left(\phi\right)\mathcal{P}^{-1}=-\mathcal{H}\left(-\phi\right)$ where $\mathcal{P}=\sigma_{y}\tau_{y}\mathcal{K}$ and $\mathcal{K}$ denotes the complex conjugation. It also preserves a generalized time-reversal symmetry $\mathcal{TH}\left(\phi\right)\mathcal{T}^{-1}=\mathcal{H}\left(-\phi\right)$ with $\mathcal{T}=\sigma_{z}\mathcal{K}$. The combination of $\mathcal{P}$ and $\mathcal{T}$ leads to the chiral symmetry $\mathcal{C}$: $\mathcal{CH}\left(\phi\right)\mathcal{C}^{-1}=-\mathcal{H}\left(\phi\right)$, with $\mathcal{C}=\mathcal{PT}=i\sigma_{x}\tau_{y}$. Therefore, the Hamiltonian belongs to $\mathrm{BDI}$ class and can be characterized by a $\mathbb{Z}$ topological invariant Altland1997 ; Schnyder2008 .

In Eq. (11), the SOC $\alpha\left(\phi\right)=\tilde{\alpha}_{0}\sin\phi$ changes sign at  $\phi=0,\pi$, as shown in Fig. 5 (b). Specifically, we have $\alpha\left(\phi\right)>$ $0$ if $\phi\in\left(0,\pi\right)$ and $\alpha\left(\phi\right)<0$ if $\phi\in\left(\pi,2\pi\right)$. Hence, the system can be divided into two segments. Both belong to the $\mathrm{BDI}$ class but possess opposite topological invariant. The interfaces are determined by $\phi=0$ and $\pi$, corresponding to two “soft” domain walls in the sense that the SOC term changes smoothly across these two points. From Eq. (11), the Hamiltonian $\mathcal{H}\left(\phi\right)$ is invariant under a $2n\pi$ rotation if $n$ is an integer. Therefore, to solve the eigenvalues of $\mathcal{H}\left(\phi\right)$, we assume the following trial solution,

where $u_{\nu=a,b,c,d}\left(\phi\right)=\sum_{m}\nu_{m}e^{im\phi}$ and $m$ is an integer. By solving the Schrödinger equation $\mathcal{H}\left(\phi\right)\Phi\left(\phi\right)=E\Phi\left(\phi\right)$, the eigenvalues are obtained as shown in Fig. 5 (c). See Appendix 7 for more details. It is clear that four Majorana modes emerge (with an numerical error about $E\approx 10^{-4}$). One Majorana corner pair consisting of two Majorana modes localizes at $\phi=0$, and the other pair localizes at $\phi=\pi$, as illustrated in Fig. 5 (a). This is also demonstrated by the particle density distribution $\rho_{{}_{MZ}}$ of Majorana modes, as shown in Fig. 5 (d). We remark that a toroidal Bose-Einstein condensate has been created in an all-optical trap Ramanathan2011 . We expect our scheme could be reached with similar techniques and additional Raman lasers.

From the effective low-energy theory of TSFs, it is well-known that Majorana modes would emerge if the sign of the Dirac mass changes and most previous proposals are based on this principle. In this paper, we propose an alternative approach to implement Majorana modes (Majorana corner pairs) through tuning the effective SOC. By loading Fermi gases on 2D optical lattices subjected to a 1D ERD SOC, we can find a SF phase under appropriate $s$-wave interaction and Zeeman field. Using single-site addressing techniques, we could engineer defective geometries, which are topologically non-trivial, on the 2D optical lattice. From the viewpoint of low-energy theories, a defect rectangle consists of two TSFs characterized by distinct topological invariants whose sign is determined by the sign of SOC in edge coordinate. Obviously, the sign of SOC changes at two corners on the defect rectangle. At the interface of two distinct TSF, a topologically protected Majorana pair naturally arises according to the index theorem. For TSF with 1D ERD SOC on a ring, two soft SOC domain walls exist, and two Majorana pairs also appear near the domain walls. In principle, as long as two effective 1D SFs are topological with different topological invariants $w=\pm 1$, the Majorana pair will emerge at the interface. It is robust as long as the perturbations preserve three underlying symmetries ($\mathcal{P}$, $\mathcal{T}$, $\mathcal{C}$) of the system.

We emphasize that the Majorana corner pair in the context differs from those in second order TSCs in two dimensions. First, for second order TSCs with time reversal symmetry Yan2018 ; Wang2018 , 1D edge modes evolve from the higher-dimensional bulk of the topological insulators. However, in our scheme, the 1D modes originate from the defect geometry. Second, in a higher-order TSC, a momentum-dependent SC pairing ($s_{\pm}$ or $d$-wave) leads to Dirac mass kink at the corner of the sample, and then induces the Majorana Kramers pair. In contrary, our proposal utilizes the sign reverse of effective SOC on the edges and lacks Kramers degeneracy.

In summary, we propose a distinct scheme to implement Majorana pairs in an atomic platform. The coordinate of Majorana pair depends on the position of SOC domain wall which can be tuned by the directions of the Raman laser beams. Moreover, our system is free of dynamical instability such that the MZM has a longer lifetime. Our work opens the possibility of implementing robust Majorana pairs and the associated non-Abelian braiding in cold atoms.