Numerical analysis of spin-orbit coupled one dimensional Fermi gas in the magnetic field

Numerical analysis of spin-orbit coupled one dimensional Fermi gas in the magnetic field

Y.H. Chan Institute of Atomic and Molecular Sciences, Academia Sinica, Taipei 10617, Taiwan

We use the density matrix renormalization group method(DMRG) and the infinite time evolved block decimation method(iTEBD) to investigate the ground states of the spin-orbit coupled Fermi gas in a one dimensional optical lattice with a transverse magnetic field. We discover that the system with attractive interaction can have a polarized insulator(PI), a superconducting phase(SC), a Luther-Emery(LE) phase and a band insulator(BI) phase as we vary the chemical potential and the strength of magnetic field. We find that spin-orbit coupling induces a triplet pairing order at zero momentum with the same critical exponent as that of the singlet pairing one in both the SC and the LE phase. In contrast to the FFLO phase found in the spin imbalanced system without spin-orbit coupling, pairings at finite momentum in these two phases have a larger exponent hence do not dictate the long range behavior. We also find good agreements of the dominant correlations between numerical results and the prediction from the bosonization method. The presence of Majorana fermions is tested. However, unlike results from the mean field study, we do not find positive evidence of Majorana fermions in our system.

I Introduction

Spin-orbit coupling(SOC) in both the condensed matter and the cold atom system has drawn a lot of attention in recent years. It plays a crucial role in the spin Hall effectsMurakami et al. (2004), spintronic devicesŽutić et al. (2004), and novel topological phasesJackeli and Khaliullin (2009); Pesin and Balents (2010). In these systems strong SOC modifies the band structure and the band topology. Synthetic SOC in the cold atom systemDalibard et al. (2011); ZHAI (2012) also gives rise to novel pairing states such as Flude-Ferrell state, and Larkin-Ovchinnikov state with topological properties, where Majorana fermion is found in the soliton excitation in the theoretical investigationWu et al. (2013); Zhang and Yi (2013); Qu et al. (2014); Xu et al. (2014). The interplay between SOC and interaction effects has brought interesting physics in both fields.

The development in the cold atom experiments in the optical lattice has provided new insight in one dimensional quantum many-body systems. In particular there are renewed interest in 1D fermi gasGuan et al. (2013). Novel pairing statesSun et al. (2014); Zhang and Yi (2013); Xu et al. (2014), large spin physicsHo and Yip (2000), dimensional crossoverSun and Bolech (2013) and the dynamical properties has been widely discussed theoretically and experimentally. Thanks to the powerful technique such as Bethe ansatz, bosonization method, and various numerical tools developed for 1D systems a better understanding and a thorough comparison with experiments are possible.

Recently, great efforts has been focused on spin-orbit coupled 1D fermi gas due to the proposed realizationOreg et al. (2010) of Majorana fermion in the nano wire system. Majorana fermion, which has the application as a basic unit in the fault-tolerant topological quantum computingSau et al. (2010); Alicea et al. (2011), is predicted in a 1D -wave superconductorKitaev (2000) yet its direct implementation is not achieved so far. The proposed equivalent model consists of the spin-orbit coupled nano-wire in a transverse magnetic field with the pairing potential induced by the proximity effect of an adjacent -wave superconductor.Oreg et al. (2010) Its experimental implementation was soon established and the evidence of Majorana fermions has also been claimed.Mourik et al. (2012) On the cold atom side implementations of SOC was realized in the bosonic systems.Lin et al. (2009, 2011)The progress in the fermionic systems is also encouraging. Two groups have successfully implemented the SOC with cold fermions and the single particle spectrum can be clearly observed in their experiments.Cheuk et al. (2012); Wang et al. (2012) However, it is not clear if the proximity effect, as a crucial ingredient of Majorana fermions in this proposal can be induced in the cold atom systems.

In contrast, the attractive interaction tuned by Feshbach resonance seems to be a more straightforward option for the pairing mechanism in the cold atom system. Several theoretical studies have shown that in a similar two dimensional model the attractive interaction can stabilize the - superconducting phases on the mean field level.Zhang et al. (2008); Sato et al. (2009) For a one dimensional system, it has also been shown to successfully stabilize topological phases in a harmonic trap.Wei and Mueller (2012) However, another study using the bosonization method reaches a different conclusion, in which the interaction merely shifts the value of Luttinger parameters.Fidkowski et al. (2011) They further consider a slightly different scheme where a spin-orbit coupled wire is in contact with a -wave superconductor wire with a power-law decaying order. The main conclusion of their study is that this model with a single spin-orbit coupled wire does not exhibit Majorana degeneracy and the Majorana zero modes can be obtained only with more than two wires.

One recent work also studied a similar system with a longitudinal magnetic field using the Bethe ansatz exact solution and the conformal field theory.Zvyagin and Schlottmann (2013) They obtained the critical exponents of the superfluidity and density waves over a wide range of SOC and found that the superfulid correlation always has a smallest exponent. Their results also show the system has Flude-Ferrell-Larkin-Ovchinnikov(FFLO) type instability with a spin-orbit interaction dependent wave vector which may develop into the FFLO phase when an array of such chains coupled together as in the experiment setup.

Motivated by these results we apply the unbiased numerical methods to solve the full interacting model in one dimension over a wide range of interaction strength. We also study the dominant correlations with the bosonization method. Our numerical results are obtained using the iTEBDVidal (2007) and the DMRGWhite (1992, 1993); Schollwöck (2005) methods. The iTEBD method is convenient for studying the phase diagram in the thermodynamic limit. We keep a virtual dimension up to 1000 and apply the imaginary evolution to reach the ground state. On the other hand the DMRG method is known to be accurate for finite size systems and allow us to determine spin excitation gaps in the case without the external magnetic field. We apply the DMRG method for systems up to 200 sites. The symmetry for particle number conservation is considered. In the DMRG method we keep virtual states and apply seven sweeps. The discarded weight in the eigenvalue of the reduced density matrix with this set of parameters is on the order of in the worst case.

Our main finding is that two different superconducting phases are stablized besides insulating phases found in the strong field or high chemical potential region. Both superconducting phases have a zero momentum triplet pairing order which has the same critical exponent as that of the singlet pairing order. Finite momentum pairing terms as those found in the FFLO phase of the spin imbalanced system without SOC are not dominant in both superconducting phases. We also investigate the phase space which was suggested to be a topological phase with Majorana fermions in the mean field study but no clear evidence of the Majorana fermion is found in our numerical results.

The rest of this work is organized as following: in Sec. II we introduce the model of interest. The phase diagram of the non-interacting system is shown for comparison with that of the interacting one. In Sec. III we follow an earlier workGangadharaiah et al. (2008) to derive the bosonization Hamiltonian and give the expression of the important order parameters. The phase diagram of the interacting system is given in Sec. IV. We also show correlation functions of the different phases and compare their long range behavior with the prediction from the result in Sec. III. Finally, we conclude our study in Sec. V.

Ii model

We considered a spin-orbit coupled fermionic atoms in a transverse magnetic field confined in a one-dimensional optical lattice. The Hamiltonian reads


where is the hopping strength, is the Rashba SOC strengthBychkov and Rashba (1984), is magnetic field in the -direction, is the chemical potential and is the density operator. In this study we take as our energy unit and investigate the (attractive) case.

Figure 1: The phase diagram of the non-interacting system in the - plane. PI and BI stand for a polarized insulator and a band insulator respectively. M4F, M2F, and M4FU are metal phases with four Fermi points, two Fermi points, and all four Fermi points in the upper bands.

In the absence of interaction Eq. 1 can be diagonalized to obtain the energy eigenvalues


and the momentum dependent eigenspinors,




, where is defined as


The spin direction of the eigenspinor has a dependence on . Away from , the eigenspinor has a spin component largely aligned in the direction at and in the direction at in the small case while the spin of the eigenspinor aligned in the opposite direction. The nonzero spin overlap between these two bands at the same energy as a consequence of the interplay between SOC and transverse magnetic field introduces triplet spin pairing at zero momentum as we turn on the attractive interaction.

The momentum at Fermi level in the upper band and lower band is denoted by and respectively. We also define




In Fig. 1 we show the phase diagram of the noninteracting system. It is known that the SOC term lifts the spin degeneracy and shifts the band energy minimum to a finite momentum. The transverse magnetic field opens up a gap at the crossing of the two bands at zero momentum. Therefore, depending on the chemical potential (or filling number ) the noninteracting system can go through several phase transitions. We focus on the phase above half-filling since the physics below half-filling is similar due to the particle-hole symmetry. At half filling and we have a phase with four Fermi points, two in the lower band and two in the upper band. As increases such that the Fermi surface sits inside the gap the system enters a phase with two Fermi points in the upper band. Further increasing the filling leads to a phase with four Fermi points all in the upper band. When goes beyond the maximum of the upper band all bands are occupied we thus have a band insulator(BI). In the upper left corner of the phase diagram lies in between the gap opened by a strong magnetic field, which leads to a fully occupied lower band and an empty upper band. In this case we have an insulator at half-filling due to the strong magnetic field and the system has a high magnetization. We will denote it as a polarized insulator(PI).

In the next section we discuss the interaction effect on the system. The case in which there are two Fermi points has been discussed with the bosonization method before.Fidkowski et al. (2011) Here we focus on the filling with four Fermi points. With repulsive interaction a bosonization study suggests a spin density order in the direction of the spin-orbit axis.Gangadharaiah et al. (2008) Following their work we rederive the expression of correlation functions in section III and compare with our numerical results in section IV for the attractive interaction.

Iii Bosonization prediction

We study the interacting system by first linearized the spectrum near the Fermi points, . We focus on the case when both upper and lower band are partially occupied. The case at high filling number where the lower band is fully occupied and all four Fermi points are in the upper band is similar. Due to the difference of spin components in these two bands the dominant momentum component in the correlation function changes but the physical picture remains the same.

The fermionic operator can be written in terms of left-going particles and right-going particles of bands. Gangadharaiah et al. (2008)


We then follow the convention used in the Ref. 33 to bosonize the fermionic operators with


where is the short-distance cutoff and are the Klein factors with . The chiral are expressed in terms of and its dual as


and the chiral densities are written in terms of


The Hamiltonian can then be written in terms of charge and pseudo-spin modes defined as


The quadratic part of the Hamiltonian reads


Following Ref. 33 the only relevant perturbation term is


which comes from the process of a pair of oppositely moving particles in the same subband scattering in to the other subband. The factor describes the spin overlap in the process. Since this term is the only relevant perturbation it determines the phase of the system as we can see in the following discussion.

The Luttinger liquid parameters, and velocities are obtained in the first order of as


where and is assumed.

For the purpose of renormalization group(RG) analysis we define the initial value of interaction parameters as


The RG equations obtained is the standard KT equations35; 36:


Depending on the value of the interaction strength and the RG flow can go to strong coupling or a Luttinger liquid(LL) phase. In Fig. 2 we plot the RG flow of Eq. 18 and initial values of and for several interaction strength at a filling number and . We observe that for these parameters and an attractive interaction the RG flow goes to the strong coupling limit. In fact for our choice of the value is always positive in the M4F region as long as , which means the system will always flow into the strong coupling limit as we can see in Fig. 2. Although for , will have a negative value, we found in these cases is always larger than . Therefore, we conclude that there is no transition to the LL phase in the M4F region when the interaction is turned on.

Figure 2: (Color online) Blue solid lines are RG flows solved from Eq. 18. Red dots are initial values of and at , , and with a step .

In order to determine the dominant order of the system we write down here the expression of several correlations in their bosonization forms. The coefficients of different momentum components depend on the value of in Eq. 5. Here we approximate and list only the momentum component with the largest coefficients for each order parameter in the low case where


We note that for general values of , the system can have slow decaying correlations with several different momentum components. This complicates the behaviour of the correlations. For the purpose of comparing with our numerical results in the next section we will only focus on the small field case here. The finite momentum component of spin density() and charge density() terms with the largest coefficients are


The singlet() and triplet() pairing both have a zero momentum term and finite momentum terms


When the RG flow goes to the strong coupling limit in order to minimize the energy, the field in Eq. 15 will take the semiclassical minima or for positive or negative respectively. With attractive interaction we find from the above expression that the system will have non-oscillating singlet and triplet pairing order with the same smallest exponent . The finite momentum terms in the singlet pairing signal a FFLO type instability. However, these terms always have a larger exponent when the system goes to the strong coupling limit. This is quite different from the case without SOC, where FFLO is found in a considerable part of the phase diagram.37 The sub-dominant order in the system is the charge-density wave with an exponent . We will compare these prediction with our numerical results in the next section.

Iv Phase diagram

Figure 3: (Color online) The phase diagram of the interacting system with . SC stands for a superconducting phase. LE is a pseudo-spin gap superconducting phase.(see context)

Our results of the attractively interacting system can be summarized in Fig. 3. We map out the phase diagram in the - plane with attractive interaction . Due to the particle-hole symmetry of the Hamiltonian the system is at exact half-filling when . We will focus on the phase diagram above half-filling. The phase diagram less than half-filling can be obtained by performing a particle-hole transformation.

Two insulating phases similar to those in the non-interacting system are located at two corners. We find a polarized insulator(PI) in the upper-left corner of the phase diagram. Due to the strong transverse field, a gap opened between the upper and the lower band. The ground state on a single site is close to with filling number 1. On the lower-right corner we have a band insulator(BI) with both bands fully occupied.

In the middle region a fully occupied lower band and a partially filled upper band developed into a Luttinger liquid phase when we turn on the interaction. Since this phase has a dominant correlation with superconductor orders we will denote it as a SC phase.

In the lower-left corner we find a phase transition from the SC phase to a another superconducting phase. We identify this phase as the gapped superconducting phase discussed in the earlier section. We will call it as the Luther-Emery phase(LE) due to the similarity of the gapped nature. This phase is developed from the non-interacting M4F region. The M4FU region in the non-interacting system is now connected with the LE phase at small when a finite interaction is turned on. The BI phase at finite extends to the line while another superconducting phase at separates itself from the LE phase. We will show that this phase has a well-defined spin gap and its correlations behave differently from those in the LE phase.

Figure 4: (Color online) Order parameters and as functions of at (blue dot) and (red square). The phase boundary is determined from the discontinuity of slopes.

We determined the phase boundary by calculating the order parameters and as a function of the chemical potential at different field strength with the iTEBD method. In Fig. 4 we show the results at and as an example. At we observe a discontinuity of the slopes of both order parameters at , 1.4, and 1.8, which signals phase transitions at these points. The phase beyond has zero magnetization and a filling number 2. Thus it is a BI. Other three phases show quite different behavior of compressibility and magnetization as changes and do not have simple features to be identified. Similarly we can identify a single transition at in the line. The magnetization saturates to a fixed value when and the filling number is 1, which suggest it is the PI phase as that in the non-interacting system.

iv.1 LE phase

In order to understand those phases other than PI and BI we compute several correlation functions to characterize them. We define the density-density correlation functions,


spin component correlation function,


the pairing correlation function




for on-site singlet pairing, and


for triplet pairing. The corresponding structure factors are Fourier transformation of the above correlation function,


where is the system size.

Figure 5: (Color online) Structure factors of (a) density , (b) spin density , and (c) spin density at (blue dot), -2(red square), -3(green cross), and .
Figure 6: (Color online) Structure factors of (a) singlet pairing and (b) triplet pairing at (blue dot), -2(red square), -3(green cross), , and .

In Fig. 5 we plot the density wave and the spin density structure factors at , with , -2, and -3 for , which belongs to the phase in the bottom-left corner in Fig. 3. We observe in Fig. 5(a) a cusp develops at momentum of density-density structure factor, which signals a charge-density wave as we increase the attractive interaction to . density structure factor has a kink at and a dip at momentum. The kinks at are smoothed as increases as shown in Fig. 5(b). Meanwhile, we find in Fig. 5(c) the structure factor of density are suppressed in all three interaction strength. From these results we can infer that the dominant order in the density wave part is the charge density one.

In Fig. 6(a) we observe the singlet pairing structure factor develops a cusp at as we increase the strength of attractive interaction, which suggests a strong singlet pairing order at large . Together with the density wave order this corresponds to a strong coupling Luther-Emery(LE) phase discussed in III. We will further study their real space correlations below. A kink at momentum is also observed in the singlet pairing, which also agree with the bosonization result. In Fig. 6(b) we can find a bump at momentum in the triplet pairing structure factors at all interaction strength. However, the supposed dominant component predicted by the bonsonization analysis does not have a prominent peak. This can be understood by a smaller coefficient of the component as we will see when the correlations are plotted in the real space.

The bosonization method suggests with any nonzero attractive interaction the system will flow into strong coupling limit. The exponents of dominant correlation functions can then be determined by a single Luttinger parameters . In order to compare this conclusion with our numerical study we first extract the Luttinger parameter from the density structure factor. It can be shown that in the limit the density structure factor reads38; 39


We shows in Fig. 7(a) the linear fitting result of the density structure factor near at , , and with and . We find that finite size effects do not significantly change the fitting results. The values extract from these two system sizes agree with each other. Thus, we only show the results for in our data below. In Fig. 7(b) we compare the values from numerical results with those calculated in Eq. III and find that they match reasonably well within the interaction strength regime we studied.

Figure 7: (Color online) (a) Structure factor of (blue dot) and (red square) at small . The linear fitting result is shown as a solid line. (b) Luttinger parameter obtained from Eq. III(solid line) and from linear fitting results(blue dot) as a function of .

We now turn to the real space correlations. It is known that the DMRG calculation of real space correlations may suffer from both finite size effects and Friedel oscillation from fixed boundary condition40. One solution widely used is to average over correlations in only the center part of the system, which is taken to be far from the boundary. On the other hand the iTEBD method exploit the translational symmetry of the system and can be used to calculated correlations without worrying about the boundary and finite size effects. In a recent work it has been used to study the correlations of the spinless fermion system, in which good agreement with the bosonization method is obtained41 In the following we will mainly present the correlations from the iTEBD results. The DMRG results will serve as a comparison in some cases.

Figure 8: (Color online) Real space correlation functions (a) (blue dot), (red diamond), (green square) and (b) in log-log scale at , , and . The red dashed line in (b) is a polynomial fitting function with the slope determined from Eq. 32.

In Fig. 8(a) we plot the spin density correlation functions in log-log scale at , , and . We find that both and spin components have a correlation which decays exponentially in the short range but has a polynomially-decaying tail while the correlation decays exponentially. Furthermore, both and modulated with the same periodicity at short distance and has a different modulation similar to that of the density correlation at long distance.

This behavior can be understood from the bosonization results. Since takes the value in the strong coupling limit, the density wave, which has a bosonization form , will have a zero expectation value. This leads to a correlation which decays exponentially. The fluctuation of the conjugate variable also leads to a zero expectation value of , which explains the short-range behavior of them. Besides the short-ranged term, both of them have a term with smaller coefficients described by


In the strong coupling limit takes a finite value and both terms decay polynomially. However, since , these terms show up only after a certain distance where the exponentially-decaying terms vanish. The bosonization result also explains different modulation periods at the short range and the long range. We also check this behavior at stronger interaction (not shown), where we find a similar result. Due to the shorter correlation length at , the terms start to dominate at a shorter distance.

The charge density correlations and singlet pairing correlation are now given by


since the in Eq. 23 takes it semiclassical value in the LE phase. In Fig. 8(b) and 9 we show the correlation functions and their asymptotic forms. The used in the asymptotic forms is obtained from Eq. 32. The constant is determined by fitting the maximum values of the correlations in the oscillations. We can clearly see that the power law behavior is well-described by Eq. 34. Furthermore, the same periodicity of the charge density correlation and the long range term of the and correlation seen in Eq. 34 and Eq. 33 also agrees with our numerical results. Similarly, the non-oscillating component of the singlet paring correlation agrees nicely with Eq. 35 as we show in Fig. 9(a).

Figure 9: (Color online) Real space correlation functions (a) and (b) determined by the DMRG (red square) or the iTEBD (blue dot) method in log-log scale with the same parameter as in Fig. 8. The green dashed lines in (a) and (b) show the asymptotic behavior determined by the bosonization method.

In Fig.9(b) we plot the modulus of the real part of the triplet correlation. We have shown earlier the non-oscillating part of the triplet correlation possesses exactly the same exponent as that of the singlet correlation from the bosonization results. Besides, it also has a momentum term with the largest coefficient at small field, which can readily be seen in Fig. 6(b). Due to the fluctuation of in the LE phase this term decays exponentially. This explains what we observe in Fig.9(b). The triplet correlation has a oscillating part which decays rather faster and vanishes after 30 sites. In the long distance the component dominates. We find the long range behavior can be well described by the power law function with the exponent as we show by the green dashed line in which the constant is obtained by fitting only part of the correlation.

We also demonstrate similar results computed with the DMRG method. For the DMRG calculation we keep states and work with a system size . In order to reduce the boundary effect only the correlations in the central 100 sites are averaged. In Fig. 9(a) and (b) we plot the singlet and the triplet correlation calculated from both the iTEBD and the DMRG method. As we can see two results agree well in the short distance while the DMRG data decay a bit faster in the long range part of the singlet correlation. We suspect this may due to the limited virtual dimension we used here and the finite size effect. Since for the iTEBD method with the same computational resource a larger number of can be used it gives better results in the long distance. Typically, we keep (300) states in the iTEBD calculation of the singlet(triplet) correlation.

Although we have a fairly good agreement between the bosonization results and the numerical one at a strong interaction strength here, we also find that at a smaller the power low exponents for the singlet pairing correlation are always more negative than , which suggests a finite value is required to fit the correlation functions. Since no signature of phase transition is found in our study as we change over a wide range we suspect our numerical methods could not capture the real exponents of the correlation functions because of the finite size effect from the limited size and the limited virtual dimension . We note that the correlation length of the system at small interaction could be quite large, which means in order to get the correct exponents both the system size and the inherent correlation length determined by the virtual dimension should be larger than the correlation length of order parameters. This can only be achieved by keeping more states in the DMRG or the iTEBD.

Our numerical results demonstrate that the LE phase is a superconducting phase with a dominant singlet pairing at momentum. Although the pairing at finite momentum does occur it has a larger exponent than that of the zero momentum term, which suggests the FFLO phase found in the spin imbalanced system without the SOC is not the main instability here in the spin-orbit coupled system. We also showed that the triplet pairing correlation has the same component but the amplitude of the correlations is several orders smaller than that of the singlet one. Hence, the short-ranged behavior of the triplet pairing is dominated by a finite momentum term. We further check this conclusion with a stronger interaction . We find that similar behaviors are seen at the same filling number and external field but the exponentially decaying terms vanishes within tens of sites. This can be understood as a shorter correlation length and a larger gap is developed with a stronger interaction. No further transition is found in the LE phase as we increase the strength up to .

iv.2 Zero field SC phase

The ground state at the line is a spin gapped superconducting phase. At zero field spin components of non-interacting bands do not mix with each other and the spin up and the spin down atoms just have the opposite momentum. Once we turns on the interaction the formation of singlet pairs dominates and brings in a spin excitation gap.

Figure 10: (Color online) Spin gaps at , and (blue dot), (green square), and (red diamond). The extrapolations have finite intercept in all cases.

To verify the existence of a spin gap at zero field we can measure it with the DMRG method at finite sizes. This can be done by assigning the as an additional good quantum number of the system. The spin gap is given by the energy difference of ground states with different , namely

where is the ground state energy of given quantum numbers and . We measure the at and and extrapolate them to the thermodynamic limit with a linear function. Figure 10 shows spin gaps as a function of at and , and . We can see that spin gaps decrease rather fast as we lower the interaction strength but remain finite in the thermodynamic limit. At the reads . This conclusion remains at the smallest interaction strength we studied where we have at .

Figure 11: (Color online) (a) Real space (blue dot), (red diamond), (green square) and (brown cross) density correlation functions. (b) (blue dot) and (red cross) correlation functions at , and .

We could roughly estimates the correlation length from spin gaps data. The estimation is given by33

The at is roughly at the order of 1, which in turns gives a correlation length at . This may justify our conclusion of the failure to reproduce the correct exponents in the correlation functions with the DMRG or iTEBD method at small interaction strength.

The spin gap can also be inferred from the long-ranged behavior of the correlation functions. In Fig. 11, we plot the spin, charge density, and pairing correlations in the real space with the iTEBD method. We found in Fig. 11(a) that all three spin density correlations decay exponentially. Since symmetry is restored without the external field, and have exactly the same correlations. The triplet pairing correlation in Fig. 11(b) also decays exponentially. These results agree with the existence of a spin gap from the DMRG calculation.

When comparing correlations with those at shown in Fig. 8 and 9, we found the ground state at zero field also has a dominant singlet pairing at and a sub-dominant charge density correlation at . Significant difference can be seen at the , density and the triplet pairing correlations. At a finite field all correlations decay polynomially at the long distance while those at zero field simply have exponentially-decaying correlations and no slow decaying tails are observed up to sites. These results remain true at a stronger interaction , in which slow decaying correlations could reveal themselves at a shorter distance due to faster decaying exponential terms if they exist. We also find that this spin gapped phase transits into the LE phase as long as a small field turns on. This is quite different from the case without the SOC, in which a spin gapped superconducting phase survives a small finite field when the system has a strong attractive interaction.37

iv.3 SC phase

Figure 12: (Color online) Luttinger parameter obtained from Eq. 36(solid line) and from linear fitting results(red dot) as functions of at and .

At a higher magnetic filed and a larger filling number the system enters into another superconducting phase. The free system has a fully occupied lower band and two Fermi points in the upper band. As discussed in Ref. 27, interaction effect merely shifts the free Luttinger parameter . We find that the dependence of on reads33


where is the fermi velocity in the upper band. In Fig. 12 we compare the dependence of on from Eq. 36 and the numerical results obtained from the structure factors of the charge density wave at and . We find a systematic overestimate from numerical results at . The relative error is about . With even stronger attractive interaction our numerical data shows a much higher than those determined by Eq. 36 we suspect the bosonization method fails in this strong coupling region.

The bosonized spin density and charge density operators have a component proportional to and respectively. The momentum component of the charge density operator is


Their real space correlations then read


Since for the attractive interaction the term dominates and both correlations have a decay exponent -2.

Figure 13: (Color online) Real space spin density waves (a) and charge density waves (b) correlations in log-log scale at , , and . The red solid lines in (a) and (b) show the asymptotic behavior determined by the bosonization method.
Figure 14: (Color online) Real space singlet pairing (a) and triplet pairing (b) correlations in log-log scale at the same parameters as those in Fig. 13. The red dashed lines in (a) and (b) are given by the bosonization method.

Not surprisingly, the most divergent susceptibility is that of the pairing sectors. The singlet and triplet pairing operators have a zero momentum term , which leads to a decay exponent in the real space correlations. With for the attractive interaction we have dominant singlet and triplet pairing orders in this phase.

We show the real space correlations of and in Fig. 13(a) and (b) at and in the log-log scale. The power law functions with exponents -2 are also shown for comparison with constants and determined by fitting the data at . It can be clearly seen that these correlations are well-described by the power law function in Eq. 38 and Eq. 39. The singlet and triplet correlations in real space are aslo demonstrated in Fig. 14. Both correlation are described by power law functions with the same exponents . The is obtained by fitting the part of the charge density wave structure factor.

The interaction term in the Hamiltonian Eq. 1 under the mean field treatment with the singlet pairing decoupling looks exactly the same as the pairing term induced by the proximity effect. It can be shown that this region is the topological SC phase in the model with the proximity effect. To check the topological nature of this phase in the current model, we test the existence of Majorana fermions following the criteria used in Ref. 42. First, Majorana fermions are energy zero modes of the system. Therefore, the energy difference in even and odd particle number sectors will vanish if there are Majorana fermions. We scan the ground state energy as a function of particle numbers with fixed and and do not observe any degeneracy in . Second, the entanglement spectrum is found to have a two-fold degeneracy in the topological phase.43, 42 We expect it would also be the case if there were Majorana fermions in our system. However, we do not find any evidence of such degeneracy. Our results suggests that Majorana fermions might not exist in the current system.

V Conclusion

To sum up, we have numerically studied the phase diagram of a spin-orbit coupled interacting Fermi gas in the 1D optical lattice with a transverse magnetic field. We find that the intrinsic attractive interaction in the system can not induce the Majorana fermions, unlike the more well-known example in which a pairing is induced by the proximity effect. Instead, besides a half-filled insulating phase in the strong field limit and a band insulator we identify a LE phase in the weak field and another superconducting phase in the strong field.

The LE phase has dominant singlet and triplet pairing orders. Although they decay with the same exponent, the singlet pairing term has a larger coefficient. The sub-dominant correlation is that of the charge density wave. The Luttinger parameters extracted from our numerical results matches reasonably well to those determined by the bosonization method. The correlation functions obtained from numerical methods are in good agreement with the asymptotic forms predicted in the bosonization method at strong interaction. We also find pairing terms with finite momentum. However, in contrast to the FFLO phase found in the case without the spin-orbit coupling these terms always decay faster than the zero-momentum pairing term. We compute spin gaps at limit by the DMRG method. Our results show relatively small spin gaps at weak interaction, which suggest a large correlation length in this region and may cause our failure to capture the correct exponents of the correlation functions in this limit.

The mean field phase at high field region was suggested to be a topological superconducting phase hosting Majorana fermion. In our numerical study we again find a superconducting phase with equally dominant singlet and triplet pairing order. This phase can be well-described by one Luttinger parameter in all the attractive interaction strength we studied, which agrees with an earlier bosonization study27. To search for evidence of Majorana fermions we follow the criteria used in an earlier DMRG study.42 We compute the ground state energy difference in the even and odd parity sectors and the entanglement spectrum at the cut of half the system. However, no supportive evidence is found in our study. It would be an interesting future work to find numerical evidence of the Majorana fermion in a ladder system as suggested by Ref. 27 or systems with long range interaction44.


We are grateful for the discussion with Ching-Kai Chiu, Professor Kai Sun, Professor Luming Duan, and Professor Jesko Sirker. We especially appreciate Kuei Sun’s suggestion on the manuscript. We acknowledge support by a Thematic Project at Academia Sinica. Part of the work has been developed by using the DMRG code released within the ”Powder with Power” project (


Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description