Topological quantum phase transitions and edge states in spin-orbital coupled Fermi gases

We study superconducting states in the presence of spin-orbital coupling and Zeeman field. It is found that a phase transition from a Fulde-Ferrell-Larkin-Ovchinnikov state to the topological superconducting state occurs upon increasing the spin-orbital coupling. The nature of this topological phase transition and its critical property are investigated numerically. Physical properties of the topological superconducting phase are also explored. Moreover, the local density of states is calculated, through which the topological feature may be tested experimentally.

T he subject of topological phases in quantum systems has been studied intensively for the past few years due to their nontrivial properties 1 , including significant theoretical research interests on topological superconductors/superfluids (TSC) [2][3][4][5] . A TSC is characterized by a full pairing gap in the bulk and has topologically protected gapless states on edges of the system. Soon after the theoretical prediction, the TSC material was reported experimentally through doping Cu into Be 2 Se 3 (Cu x Be 2 Se 3 ) [6][7][8][9] .
One of the most important features of TSC is the existence of gapless edge states. Especially, the zero energy edge states are usually related to the Majorana Femions (MFs) 1 , who are their own antiparticles. MFs usually obey non-Abelian statistics and have a potential application in quantum computation, so that the realization of MFs is of significant interest.
As is known, a typical TSC is the p 1 ip pairing system 1 . Besides the p-wave pairing, the TSC could also be realized in the system with an s-wave pairing plus the spin-orbital coupling and spin polarization. Actually, it has been shown that the s-wave pairing system with a strong spin-orbital coupling can be equivalent to a p 1 ip system [10][11][12][13][14] . Thus, apart from the intrinsic TSC material such as Cu x Be 2 Se 3 , there may also exist some artificiallymade TSCs. A kind of candidate systems are ultracold atoms, in which both the s-wave superfluidity and spinorbital coupling have been realized experimentally [15][16][17][18][19] . Meanwhile, it was also proposed theoretically that the topological phase and MFs can be realized in a strong spin-orbital coupled material, through coupling to an swave superconductor and an effective Zeeman field 20,21 . Recently, such kind of proposal has been realized experimentally and the signatures of MF excitations have been reported by several groups [22][23][24][25][26] .
Another quite arresting feature in superconducting systems seems to be the possible appearance of Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) states. This kind of states were predicted theoretically in 1960s 27,28 for some superconductors subject to a strong magnetic field, namely, the Cooper pairs have a finite momentum and thus the order parameter varies periodically in real space. However, this long-sought inhomogeneous superconducting state has not been observed experimentally for quite long time, possibly due to the impurity effect 29 or the orbital effect induced by the field 30 .
For low dimensional superconducting systems, the orbital effect is negligibly weak when the magnetic field is parallel to the superconducting layers. Thus they can be promising candidates for realizing the FFLO states. In the past decade, indications for possible FFLO states have been reported in several families of materials [31][32][33][34][35][36][37][38] . Meanwhile, a signature of the FFLO state was also observed in one-dimensional partially spin polarized cold atoms 39 . The properties of FFLO states have attracted renewed interest due to the experimental breakthrough.
As mentioned above, the TSC may be generated in s-wave superconductors with the additional spin orbital coupling and Zeeman field. While the conventional FFLO state is also expected to be generated by applying a Zeeman field on s-wave superconductors, the difference between them is whether the spin-orbital coupling exists or not. Thus it is quite intriguing to study the effect of the spin-orbital coupling on the FFLO state 40 . This issue has attracted considerable interest very recently [41][42][43][44][45][46] . The spin orbital effect would break the FFLO modulation and the topological phase shows up for a moderate-strength spin orbital coupling. Generally, there may exist a phase transition from the FFLO state to the TSC state. The nature of such transition and physical properties near the critical point are of fundamental interest. Moreover, they should be quite important for exploring potential applications of topological phases. In this paper, motivated by the above consideration, we study theoretically the interplay between the spin orbital effect and the exchange field. Based on a lattice model, the nature of the phase transition from the FFLO to TSC is explored numerically. Our results show that there exist two regions of topological states with the chemical potential being m near 24 or 0. Interestingly, near the region of m 5 0, there are two effective energy bands crossing the Fermi energy. This is essentially different from the case of the continuous model 41 . The existence of the two bands is important to realize the FFLO state. Thus we propose that the phase transition from an FFLO state to a topological state may be easier to occur in the lattice model. The properties of topological phase and the MF excitations are also investigated in detail. We also study the local density of states (LDOS) to compare with the experiments.

Results
We start from a model Hamiltonian that includes the spin-orbital coupling, the Zeeman field coupling, and the s-wave pairing term, which is given by with y i 5 (y i" , y i# ) T , where s n are the identity (n 5 0) and Pauli matrix (n 5 1, 2, 3), respectively. l is the spin-orbital coupling strength and h represents an effective Zeeman field. Here we study different phases and the phase transition at zero temperature. If not specifically indicated, the parameters are chosen as h 5 0.6 and m 5 0, D i is obtained self-consistently with the pairing potential V 5 2.4 (See methods). We first illustrate numerically the topological feature of the present model. In the topological state, if we consider the edge effect to occur along the x-direction and the order parameter to be uniform along the y-direction, the two-dimensional Hamiltonian can be reduced to a quasi-one dimensional model [See methods, Eqs. The energy bands are plotted in Figs. 1(a) and 1(b). As is shown, for both spin-orbital coupling strengths we considered, there exists an energy gap and one gapless state. We have checked numerically that the spatially distribution of the gapless states are at the system edge. The existence of the edge states indicates the topological feature of our present model. The gap magnitude depends on the spinorbital coupling strength; it decreases significantly as the spin orbital coupling decreases.
The order parameters with different spin-orbital coupling strengths are plotted in Figs. 1(c) and 1(d). As is seen, they are uniform in the bulk and exhibit a boundary effect near the edge. For the stronger spin orbital coupling strength l 5 0.5, an obvious boundary effect occurs only within the ten lattice constants away from the edge [ Fig. 1(c)]. While when l decreases, the size of more lattice constants is affected by the boundary. For l 5 0.3, as is seen in Fig. 1(d), the order parameter oscillates significantly, with the magnitude of the oscillation decreasing when away from the boundary. At more than one hundred lattice constants from the boundary, the order parameter recovers to the bulk one. Previously, it was also reported that the wave function may oscillate near the vortex core and system edge 47 . As proposed in Ref. 47, it might make it difficult to realize a single qubit gate for universal topological quantum computation using the tunneling between two vortices. On the other hand, such oscillation is of fundamental interest. As the spin-orbital coupling decreases further, we expect that the oscillated lattice should increase further and expand to the whole system at certain critical spin orbital coupling strength. Then the system will transit to the FFLO state for weaker spin orbital coupling. The above expecta- www.nature.com/scientificreports SCIENTIFIC REPORTS | 4 : 5218 | DOI: 10.1038/srep05218 tion will be verified by the numerical results in the quasi-one-dimensional lattice, which will be presented in the following.
To conduct a more systematic and accurate numerical study on the nature of the phase transition, we perform numerical calculations on the three-leg ladder system (with the size 400 3 3) in the following. The low energy eigenvalues of the Hamiltonian are plotted in Fig. 2 with l 5 0.5. Here we use a fully self-consistent calculation [Eqs. (3) and (4) in methods section], with the open boundary condition along the x-direction. As is seen in Fig. 2, there exist two degenerate zero-energy eigenvalues protected by a minigap about 0.16. In the BdG framework, the eigenvalues 6E usually come from one physical particle. In the case of the zero energy mode, one physical particle corresponds to one pair Majorana Fermions. The two MFs are spatially separated and protected by the minigap, thus they cannot be removed by local weak perturbations. The existence of the protected zero-energy state also confirms the topological feature of the system with the present parameter.
We elucidate numerically the nature of the phase transition driven by the spin-orbital coupling. The order parameters with different spin-orbital strengths are plotted in Fig. 3. As is seen from Fig. 3(a), the order parameter oscillates periodically in the whole space for the small spin-orbital strength (l 5 0.15). Such oscillation is different from the above one shown in the topological state [ Fig. 1(d)]. Firstly, the order parameter changes the sign within one periodical lattice, thus the phase of the order parameter changes and also varies periodically. Secondly, the magnitude of the oscillation does not change in the bulk. So the system is in the FFLO state for this weak spin-orbital coupling. As l increases to l 5 0.17, the result of order parameter changes qualitatively. As is seen, there exists an obvious boundary effect, namely, the order parameter oscillates strongly near the boundary while it oscillates very weakly in the bulk. The phase of the order parameter does not change when away from the boundary. On the other hand, the order parameter still oscillates in the whole space. Thus this coupling strength l 5 0.17 should be the critical coupling strength between the FFLO state and the TSC phase. As l increases further [ Fig. 3(c)], the oscillation region decreases and the order parameter is uniform in the bulk. For a larger spin-orbital strength, the oscillation behavior disappears and the order parameter recovers to the bulk value within several lattice-constants away from the boundary. The above results of the TSC state are qualitatively the same as those of the two-dimensional lattice shown in Fig. 1. Now let us study the properties of the TSC phase. As shown in Fig. 1, these exists a gapless edge state in this phase for the twodimensional systems. For the present model, there should be two MFs associated with the zero energy state, expressed as c 1 and c 2 .
Here the zero-energy fermion C { can be obtained from the BdG equation, thus the MF states can be studied numerically. The numerical results of the spatial distribution of the two MF states are presented in Fig. 4. At the critical point, as is seen in Figs. 4(a) and 4(b), the distributions of the two MFs are the same so that they shall annihilate to an ordinary fermion. The distribution curve oscillates in the whole lattices. And the oscillation is stronger near the boundary. This is similar as the case of the order parameter. As the spin orbital strength increases to l 5 0.2, the two MFs are completely separated and locate near the two boundaries, respectively. The distribution curve still oscillates with the depth about 100 lattice size. As the spin-orbital   strength increases further, the oscillation disappears and the two MFs are located near the boundaries. These results may be verified by experiments and may be useful when exploring the application of MFs in the topological quantum computation. Finally, let us look into the LDOS to disclose the existence and distribution of the zero mode. The LDOS spectra as a function of the energy at the boundary for different spin orbital coupling are plotted in Fig. 5(a). As is shown, at the critical point (l 5 0.17), one can see a weak zero energy peak. As the spin-orbital strength increases, the intensity of the zero energy peak increases. For the strong strength (l 5 0.5), the spectra weight of the zero energy peak is quite strong. Actually, the zero bias peak has also been observed by very recent experiment and it was believed to be signature of the MFs 26 . The LDOS spectra as a function of the real space position are presented in Figs. 5(b-d). Here the LDOS spectra may be qualitatively consistent with the distributions of the two MF states. This is understandable because the contributions from non-zero energy spectra are very small due the existence of the minigap. The LDOS spectra also oscillates near the boundary for l 5 0.17 and l 5 0.2, which might be verified by later experiments. The numerical results for LDOS may establish a useful link for theoretical calculations and experimental observations.

Discussion
We can explain qualitatively the origin of the phase transition and the oscillation of the order parameter in the topological phase. In the momentum space, the two dimensional normal state Hamiltonian  In the presence of the spin orbital coupling, the renormalized normal state energy bands are expressed as The quasiparticle operators can also be obtained through diagalizing the Hamiltonian [Eq. (2)]. Both are expressed as the superposition of the spin up and spin down electrons. As a result, the intra-band pairing becomes possible, with the weight being determined by l.
The inter-band pairing would generate the FFLO-type modulation.
The intra-band one would generate zero-momentum Cooper pair. And here the Cooper pairs are made up from two spinless quasiparticles, so that the case of intra-band pairing is equivalent to two-band p-wave superconductor with opposite chirality.
Let us first brief an idea about the topological phase based on the continuous-type model [48][49][50] , with the band dispersion expressed as ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi h 2 +a k j j 2 q . A basic point to obtain a topological phase is that the Zeeman field splits the two spin-orbital bands with a gap 2 h. One may choose the chemical potential to place the Fermi level inside the gap. As a result, the electrons only occupy the lower band, while the upper unoccupied band almost plays no role. Based on the above point, the parameter region for the appearance of the topological phase has been obtained. One can conclude from the band dispersion that the Fermi level only crosses the lower band when the chemical potential is less than the exchange field (j m j , h). Thus the system is equivalent to a one-band p 1 ip pairing system. With a pairing gap D and an effective exchange field h, the parameter region for the topological phase is hw ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi D 2 zm 2 q [48][49][50] .  For the lattice-type model, the energy bands are expressed by Eq. (3). The topological behavior should be determined by the four points (0, 0), (p, 0), (0, p), and (p, p). For the case of quasi-one dimensional system with finite odd lattice size along y direction, the topological behavior is determined by the C 5 (0, 0) and X 5 (p, 0) points. We plot in Fig. 6 the band dispersions for E 6 along the k y 5 0 and k y 5 2p/3 cuts to discuss the nature of the topological phase in our system. Here the C point is the bottom of the energy bands with the energies 24 6 h. If one sets the Fermi level between them, then only the lower band crosses the level. As a result, this system is equivalent to a one-band p 1 ip system, similar to the continuous model as discussed above. The topological state near m 5 24 has been studied previously 11 and is not concerned with in the present work. On the other hand, the energy reaches the maximum at X point along k y 5 0 cut. If we set the chemical potential to satisfy jmj, h, then the Fermi level is also inside the gap between the lower band and upper band. With a particle-hole transformation, the upper band is occupied with holes, while the lower band is unoccupied. This leads to the topological feature of the system. However, here the system is indeed an effective two-band system. Actually, the lower band also crosses the Fermi level along the other direction of the Brillouin zone, e.g., along the k y 5 2p/3 direction, as seen in Fig. 6(b). This would not affect the topological feature of the system. The existence of two sheets of Fermi surface for the case m 5 0 is essential for the appearance of the FFLO state when the spin-orbital coupling strength decreases or disappears.
Then one can understand the nature of the phase transition. For a two-band system with the bands expressed by Eq.(3), the FFLO state comes from the interband pairing. The topological state come from the intra-band pairing. The origin of the phase transition is due to the competing of these two kinds of pairing. As the spin-orbital coupling increases, the intra-band pairing dominates over the inter-band one, so that the long range FFLO order is broken and the phase transition occurs. However, even in the topological phase, the local FFLO modulation may still survive as the spin orbital coupling is not strong enough. This will make the order parameter oscillate near the edge, as shown in Fig. 3.
The appearance of the topological state for the chemical potential near 24 and 0 can be checked numerically. Here we set the pairing gap D and the exchange field h to be D 5 0.3 and h 5 0.6. We consider the open (periodic) boundary conditions along the x (y) direction. The energy bands for the 100 3 100 two dimensional lattice and 400 3 3 quasi-one-dimensional lattice are presented in Fig. 7. As is seen in Figs. 7(a-c), for the chemical potential m 5 24, the energy bands are gapped in the bulk while the gapless edge state exists at the k y 5 0 point, indicating that the system is in the topological state. The gapless state disappears for m 5 23, thus the system is in the topological trivial state. For the case of m 5 0, the gapless states appear at k y 5 0 and k y 5 p. These two gapless mode are contributed by the (p, 0) and (0, p) points. So there are two topological regions. We have also checked numerically that the system is in the tropological state if the parameters satisfy hw ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi mz4 ð Þ 2 zD 2 q , or hw ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi m 2 zD 2 q . For the quasi-one dimensional lattice, as seen in Figs. 7(d-f), the zero-mode exists for m 5 24 and m 5 0. Both are protected by a minigap. We have also checked numerically that the parameter regimes for the topological state are the same as the case of the two-dimensional lattice. And here we focus our study on the region hw ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi m 2 zD 2 q . Our main results shown in Figs. (1-5) are robust for this parameter region, which has been checked numerically (not presented here).
Recently, the interplay between the FFLO state and topological state in the presence of the c-axis Zeeman field has also been studied 41,42 . We now compare our results with theirs. A primary difference has been revealed in our above discussion, i.e., there exist two bands crossing the Fermi level in the topological state. As a result, as the spin-orbital coupling strength decreases, the topological state can smoothly connect to the FFLO state and a second order phase transition occurs. For the model used in Refs. 41, 42, the realization of the topological state requires that the Fermi level crosses only lower band. While the FFLO state requires two separated Fermi surfaces. In order to obtain a direct transition from the FFLO state to the topological state, a very strong pairing potential is usually required to overcome the gap between the Fermi level and the upper band. Thus the phase diagram should be quite different from ours. And the large pairing gap would also prevent the topological state to appear. Generally, there may exist a phase transition from the FFLO state to the topological trivial superconducting state. The direct transition from the FFLO to topological state may only occur for very large pairing potential and very large exchange field. Besides this important difference, there are some other technique differences between our work and others. Namely, it was argued in Ref. 41 that the study of one-dimensional system is problematic according to the Mermin-Wagner theorem [51][52][53] . Thus the coupling chains are considered. The FFLO modulation occurs only along the chain direction and the order parameter is assumed to be uniform along its perpendicular direction. The BdG equations are solved in the momentum space so that certain local inhomogeneous feature is neglected. This is also different from our model and actually the gap fluctuation is also an important part of the present work. Some interesting signatures for the phase transition are provided based on the present work which may be detected by later experiments. Ref. 42 mainly discussed the numerical results for the pure one-dimensional system and weakly coupled chains. Actually, the topological phase near m 5 0 region could not appear for the one-dimensional model. It appears near m 5 22, which is the bottom of the normal state band. Thus the system is also a one-band system, and a direct transition from the FFLO state to the topological state can hardly occur.
It was proposed that the topological order may form before the disappearance of the FFLO state. Thus a special region, named as ''topological FFLO state'' might appear in the phase diagram 41,44,45 . Usually the topological order can be determined strictly through studying the Pfaffian topological invariant. In Refs. 41, 44, 45, the model includes both in-plane Zeeman field and c-axis Zeeman field. Such a model would favor an FF state with D(i) 5 D 0 exp(iQ ? R i ), in which the order parameter magnitude is spatially uniform. The phase factor in the order parameter can be gauged out and the BdG hamiltonian can be transformed to the momentum space 44 . This technique can be used to investigate the Pfaffian topological invariant, and a topological FF state has indeed been verified numerically 41,44 . In the present work, only c-axis Zeeman field is considered. The quasiparticle bands are degenerate with the momentum k and 2k. For this case, usually LO state with D(i) 5 D 0 cos(Q ? R i ) is favored. The translational symmetry is broken and thus it is difficult to study numerically the possible topological order in the LO state. While if the topological phase could persist into the LO state, the possible ''topological-LO'' state is quite interesting, which may be worthwhile studying in future.
It is also meaningful to compare our results with previous numerical and analytical results on the topological states for the continuous model. A main character for the topological state is the existence of the zero energy states protected by a minigap. The magnitude of the minigap is an important parameter that can determine the robustness of the MF state against various fluctuations and disorder/impurity effect. Generally the minigap should depend strongly on the magnitude of the bulk pairing gap. With a uniform pairing gap, the minigap can approximately be obtained in the momentum space. Namely, the Hamiltonian [Eq. (4) in the method section] can be rewritten in the momentum space considering periodic boundaries. The quasiparticle bands for the Hamiltonian can be obtained numerically, denoted as E i (k). The quasiparticle bands are fully gapped and the minigap should be the minimum value of j E i (k) j. In the present work, we solve Eq.(4) in a fully self-consistent method and obtain the minigaps numerically in the real space. For a very large lattice-size, the same result is obtained for the two methods. In real space, the gap is size dependent. We define the bulk gap D 0 to be the gap at the site (200,2). Numerically as the other parameters (h, l, m) are fixed, the bulk gap can be controlled through varying the pairing potential V. We show in Figs. 8(a-c) the minigap as a function of the bulk gap. For all of the parameters we considered, the nonzero minigap exists only when the relation hw ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi D 2 0 zm 2 q is satisfied. This numerical conclusion is consistent with previous analytical results 48,49 and numerical ones for a continuous model 47 . Here the minigap increases linearly with the bulk gap as the bulk gap is small and then decreases linearly for large bulk gap. This result is also qualitatively consistent with that of the continuous model 14,47,48,49 . The scale coefficient for small bulk gap almost does not depend on the chemical potential [ Fig. 8(a)]. It decreases slightly as the exchange field increases [ Fig. 8(b)]. While it depends strongly on the spin-orbital strength l, as seen in Fig. 8(c). It increases as l increases. The relation between the minigap/the bulk gap) and the spin-orbital strength is presented in Fig. 8(d). The ratio of the minigap to the bulk gap with a fixed pairing potential V 5 2.4 is  plotted in the inset of Fig. 8(d). As is seen, the ratio increases monotonously with l and saturates to about 1.0 at a strong limit of spinorbital coupling. While the magnitudes of the bulk gap and minigap are non-monotonous as the spin-orbital strength increases. The minigap reaches the maximum for about l 5 0.4.
Finally, we would like to discuss how the present work is relevant to a real system. Generally, our model might be simulated by cold atoms. Recently, the spin-orbital coupling has been realized in one-dimensional systems with equal Rashba and Dresselhaus strengths [17][18][19] . While the pure Rashba-type spin-orbital coupling is used in the present work. Although the experimental realization of the Rashba-type spin-orbital coupling is still awaited, there have been several theoretical proposals to generate it in cold atoms very recently [54][55][56] , which seem to be quite promising. The s-wave pairing superfluid state has also been realized in cold atoms 15 . We also wish to point out that there always exist the external trap in cold atom systems, which was addressed in Ref. 57 but is not considered in the present work. The effect of the trap is also interesting and may deserve a further study. While, in principle, our main results are robust again small fluctuations, and we expect that the existence of external trap does not change our main conclusions. At last, we emphasize that we here mainly focus on the phase transition in the lattice system. The numerical results are rather different from those obtained with a continuous model. The existence of the lattice seems to play a rather important role to reach our main conclusion.
In summary, we have studied the physical properties and competing superconducting phase based on a lattice model that includes the s-wave pairing, spin-orbital coupling, and the Zeeman field term. The phase transition from the FFLO state to the TSC state induced by the spin orbital coupling has been revealed and discussed. The difference between our lattice model and a previous continuous model are discussed in detail. It has been proposed that the transition from the FFLO state to topological state is easier to occur for the lattice model. We have also explored the properties of the critical point, the minigap, and the excitations of MFs in the TSC state, providing a helpful insight in profound understanding of topological superconductivity and potential applications. The LDOS spectra have been calculated, serving as a link between our theoretical analysis and experimental observations.

Methods
The Hamiltonian can be diagonalized by solving the Bogoliubov-de Gennes (BdG) equations as where Y n (r) denotes the order column vector with Y n r ð Þ~u n r: ,u n r; ,v n r; ,v n r: T , and the order parameter D(r) is determined self-consistently D r ð Þ~V 2 X n v n r: v nÃ r; tanh with V the pairing strength.
Considering the order parameter to be uniform along the y-direction, the twodimensional Hamiltonian can be reduced to a quasi-one dimensional model, expressed as,

H~X
x where e(k y ) 5 22 cos k y 2 m and T 5 s 0 2 ils 1 . The BdG equation and self-consistently determined order parameter are written as, The local density of states (LDOS) can be calculated as The delta function d(E) is taken as d 5 x/[p(E 2 1 x 2 )], with the quasiparticle damping x 5 0.01.