Majorana edge state in a number-conserving Fermi gas with tunable p-wave interaction

The remarkable properties and potential applications of Majorana fermions have led to considerable efforts in recent years to realize topological matters that host these excitations. For a number-conserving system, there have been a few proposals, using either coupled-chain models or multi-component system with spin-orbit coupling, to create number fluctuation of fermion pairs in achieving Majorana fermion. In this work, we show that Majorana edge states can occur in a spinless Fermi gas in 1D lattices with tunable $p$-wave interaction. This is facilitated by the conversion between a pair of (open-channel) fermions and a (close-channel) boson, thereby allowing the number fluctuation of fermion pairs in a single chain. This scheme requires neither spin-orbit coupling nor multi-chain setup and can be implemented easily. Using the density-matrix-renormalization-group method, we have identified the Majorana phase in a wide range of parameter regime as well as its associated phase transitions. The topological nature of the Majorana phase manifests itself in a strong edge-edge correlation in an open chain that is robust against disorder, as well as in a non-trivial winding number in the bulk generated by using twisted boundary condition. It is also shown that the Majorana phase in this system can be stable against atom losses due to few-body collisions on the same site, and can be easily identified from the fermion momentum distribution. These results pave the way for probing the intriguing Majorana physics in a simple and stable cold atoms system.


I. INTRODUCTION
Majorana fermions, discovered by Majorana in 1937 [1], has stimulated tremendous research interests over the past decades due to their novel exchange statistics and promise for topological quantum computation [2,3]. A Majorana fermion (or a "Majorana" for short) is an equal magnitude superposition of a fermion operator and its adjoint [λ 1 = f † + f , or λ 2 = i(f † − f )]. It is a mode of excitation rather than a particle in the usual sense. In 2001, Kitaev showed that spinless fermions in a 1D chain coupled to a pairing field will have Majorana fermions at the ends [4]. Efforts to simulate this model in solid state matter have led to the proposal of using 1D semiconducting wires with spin-orbit coupling (SOC) in contact with a superconductor [5,6]. Similar proposals have also been made in the cold atom studies by engineering SOC on an attractive Fermi gas [7][8][9][10][11][12].
Since Majorana fermions also emerge in the number conserving systems such as the Pfaffian quantum Hall state [13,14] and Kitaev's honeycomb spin model [15], there have been questions of whether proximity superconductivity (or lack of number conservation) is necessary for realizing Majaronas in 1D chains. While a single Majorana excitation can not exist in a number conserving system, the correlation of two Majoranas at different locations (i λ 1 i λ 2 j ) is well defined. This provides a natural generalization of the presence of Majorana edge modes in * xlcui@iphy.ac.cn a number conserving system, which is defined as a nonzero correlation of the Majoranas at the opposite end of a finite chain, i λ 1 0 λ 2 N L = 0, in exactly the same way the Majoranas are correlated in the Kitaev model [4]. In Refs. [16][17][18], the authors have studied coupled 1D chains with interchain pair hopping using bosonization methods and have concluded the existence of Majorana edge states in these number conserving systems.
Whether Majorana edge states can exist under number conservation is particularly relevant for their realization with cold atoms, as the latter are number conserving [19][20][21][22][23]. A rigorous proof of their existence was established recently for coupled chain models with inter-chain pair hopping [19][20][21], and for a single chain four-component fermion system with SOC and spin-exchange interaction [23] which has the similar pair hopping physics as in other coupled chain models. All these studies suggest that the number fluctuation of fermion pairs in a single chain is the key to the emergence of Majoranas in number conserving systems.
In this paper, we propose a much simpler scheme for realizing Majorana edge states that makes use of alkali fermions in a single chain without SOC -by simply tuning a single component Fermi gas in a 1D chain to its pwave resonance. The systematic derivation of this model is given in Ref. [24]. In the two-channel description of pwave resonance, two fermions at neighboring lattice sites can convert to a "close channel" boson in one of the two sites, thereby causing the number fluctuation of fermion pairs in the chain. The close channel bosons play the role of proximity superconductor in the Kitaev model, except that they are now quantum mechanical objects. In this model, neither the number of fermion N f nor the number of boson N b is conserved. However, the total number N = 2N b + N f is [See Eq. (1)]. An effective fermionmolecule conversion model was also proposed previously using the laser-assisted pair tunneling [25]. Here through exact numerical calculations, we confirm the Majorana ground state with strong edge-edge correlations in a broad range of paramters : filling factor, boson detuning, and inter-channel coupling (See yellow regions in Fig. 1 and Fig. 12). It is useful to contrast our model with the single-channel model that only consists of spinless fermions with neighboring-site attractions. While both models are number conserving and their mean field theories share the same structure, our exact calculations show that only the resonant two-channel model exhibits strong Majorana correlations. This shows again the essential role played by the number fluctuation of fermion pairs in the chain, which is absent in the single-channel attractive Fermi gas.
Experimentally, a major obstacle for exploring p-wave effect in cold atoms is the severe atom loss, as observed in a 3D Fermi gas [26][27][28][29]. Recent studies have suggested that the p-wave system could be more stable against three-body loss if confined in the quasi-1D geometry [30,31]. In our system, an additional optical lattice is applied along the 1D tube and the space is further discretized. In this case, the atom loss comes from the possibility of finding pairs of boson-fermion or boson-boson at the same site and their collisions at close proximity. Here we show that within the two-channel p-wave model, the probabilities of finding such pairs are very low in a large region of parameter space for the Majorana phase. This provides promising prospects to realize the Majorana phase and to perform studies with a wide range of cold atom techniques.

II. MODEL
Our model is [32] where b † j and f † j create a closed-channel boson and an open-channel fermion at site j, with nearest-neighbor hopping t b and t f respectively; g is the (p-wave) interchannel coupling and ν is the boson detuning. As ν approaches zero, conversion between bosons and fermion pairs for given g will become more frequent due to their energy match. The number of boson ( are not separately conserved, but the sum N = 2N b + N f is. Since bosons are heavier than fermions, t b is smaller than t f . In this paper, we will take t b = 0.2t f and use t f as the energy unit. If b j in Eq. (1) is replaced by a c-number, as in mean field approach, Eq. (1) reduces to the Kitaev model [4],  which has a Majorana phase. The question is whether this phase will survive the quantum fluctuation of the bosons, such that the ground state of Eq. (1) will have the same edge-edge Majorana correlations as in the Kitaev model. This is to be answered in this work.

III. MEAN-FIELD ANALYSIS
We first carry out mean-field analysis to gain a qualitative understanding of the problem. Assuming the bosons fully condense at zero momentum, i.e., b k=0 = √ N b , the Hamiltonian Ω = H − µN in grand canonical ensemble can be written as here d is the lattice spacing; ∆ = 4g √ n b with n b = N b /N L the boson filling and N L the number of lattice sites. Equation (2) can be written as Minimizing the ground state energy of Eq. (2) with respect to ∆ and imposing the number constraint N = 2N b + k f † k f k , we obtain the gap equation and number equations as: and where E k = ξ 2 k + (∆ sin(kd)) 2 is the excitation spectrum, c b = 2n b /n is the boson fraction, and n = N/N L is the total filling. Unlike the Kitaev model where the pwave pairing and the chemical potential are both external inputs [4], here ∆ and µ are determined self-consistently for given g, ν and n.
Since the mean-field Hamiltonian Eq. (2) is in the form of the Kitaev model, the Majorana phase lies in the region |µ| < 2t f as pointed out in Ref. [4]. The topological character of the ground state is specified by the Berry Phase of the ground state χ(k) of H(k), integrated over the Brillouin Zone, i.e. ∆Φ ≡ i we have Here the phase factor e iγ(k) is to keep χ(k) periodic, As seen from Eq. (7), the tip of h(k) traces out a closed curve in the yz-plane as k varies from −π/d to π/d. For |µ| < 2t f , this curve encloses the origin, and we have θ(−π/d) = 0, and θ(π/d) = 2π. In order to ensure the periodicity of χ(k) (Eq. (9)), γ has to satisfy the condition γ(−π/d) − γ(π/d) = π, which gives the Berry phase ∆Φ = π. This is the topologically non-trivial (Majorana) phase. Otherwise for |µ| > 2t f , the curve of h(k) does not enclose the origin, and we have θ(π/d) = θ(−π/d), γ(π/d) = γ(−π/d), and thus ∆Φ = 0. This is the topologically trivial phase. After ∆ and µ are determined self-consistently, we obtain the mean field phase diagram in Fig. 1 over a broad range of ν and n for a fixed coupling g = 1. There are three superfluid phases (all with ∆ = 0): a topological Majorana phase M = SM + W M , consisting a "strong Majorana" region SM and a "weak Majorana" region W M , and two non-topological superfluid phases F HP (fermion hole pairs) and F P P (fermion particle pairs). The reasons of the nomenclature will be explained further later.
The phase boundaries separating the Majorana phase M = SM + W M and the F HP (F P P ) phase is determined by the condition µ = 2t f (µ = −2t f ), as shown by red lines in Fig. 1. Our numerical solutions of Eqs. (5) and (6) show that the upper boundary (µ = 2t f ) takes a sudden upward turn at n = 1. The red dashed line in Fig. 1 marks a crossover from SM to WM regime, when the boson fraction c b continuously decays to a small value 0.01. Later, we show from exact numerical calculations that these two regions can be distinguished from the behavior of edge-edge Majorana correlations.
The F P P phase exists in sufficiently negative ν, where the system is mostly in the Bose condensate, with a dilute gas of fermion particles forming pair superfluid. In contrast, the F HP phases exists for n > 1 and for sufficiently high ν. In this regime, the fermion occupation is more favored than bosons, leading to the nearly full fermion filling n f ∼ 1 with a small fraction of fermion holes. Since the fermion superfluidity essentially relies on the number fluctuations, it can be viewed as the pairing of fermion holes.

IV. DMRG ANALYSIS
We have calculated the ground state properties of the Hamiltonian [Eq. (1)] using density-matrixrenormalization-group (DMRG) method [33,34]. The calculations are done with maximum 800 truncated states and 30 sweeps, and the truncation error is 10 −8 . Since we consider the filling regime with n ≤ 2, we have set the truncated number of bosons at each site as up to four in our simulation.

A. Identification of Majorana phase and Majorana edge state
In this section, we use DMRG to work out the phase diagram of Eq. (1) [(I) and (II) below], and then examine in (III) the presence of Bose condensation (i.e superlfuid nature) in different phases, followed by a study in (IV) of the Majorana correlation in these phases which determines the presence of Majorana edge modes. The section will be ended by a discussion of the momentum distribution which shows the distinct signature of the Majorana phase.
(I) Entanglement entropy. It is known that a thermodynamic phase transition is reflected in a singularity of the entanglement entropy at the transition point [35]. Given the many-body ground state |ψ , the reduced density matrix can be written as ρ L = Tr R |ψ ψ|, with L, R denoting the left and right half of the lattice. Its eigenvalues {λ α } determine the entanglement entropy S = − α λ α ln λ α .
We find that the behavior of S as a function of detuning ν depends on the filling n. For n < 1, we find one cusp in S at ν = ν c1 ( Fig. 2(a)); while for n > 1, we find two cusps in S at ν c1 , ν c2 (ν c1 < ν c2 , see Fig. 2(b)). This indicates one or two transitions by tuning ν. Repeating the calculation for different sample sizes, we have obtained the estimate of the critical values of ν c1 , ν c2 for infinite systems through extrapolation (see the inset of Fig. 2(a)). These phase boundaries are shown by dots connected by solid lines in the phase diagram in Fig. 1.
Here we have chosen g = 1, and we find these phase boundaries are close to those obtained by mean field theory. Accordingly, we adopt the same nomenclature as in mean field theory for the phases obtained from DMRG.
(II) Boson fraction and its variation. We have also verified the phase boundaries by calculating the boson fraction c b and its variation c b ≡ ∂c b /∂ν as a function of ν. The results are shown in Fig. 3. One sees that while c b is continuous in ν, c b has one or two sharp cusps depending on whether the filling n < 1 or n > 1. The locations of these cusps are consistent with those obtained in (I), thus confirming the phase transitions discussed in (I). In addition, the singularity in c b shows that the transitions are of the second order. In Fig. 3, we also compare with the mean-field predictions (thin red lines) and find qualitative agreement.
(III) Condensation of bosons and fermion pairs: The criterion of Bose condensation was generalized to interacting system by Penrose and Onsager [36]. It also applies to finite number systems. The Penrose-Onsager criterion makes use of the property of the single parti- cle density matrix ρ ij ≡ b † i b j evaluated for the state of interest. The state is Bose condensed if ρ has a single maximum eigenvalue ζ 1 of order N b , while the ratios ζ β /ζ 1 for all other eigenvalues ζ β are much less than 1, and tend to zero as N b goes to infinity. The eigenfunction associated with ζ 1 is referred to as the condensate wavefunction. Similarly, one can also define a fermion pair correlation matrix Condensation of fermion pairs [as characterized by C.N. Yang as off-diagonal long-range order, see Ref. [37]] corresponds to a single large eigenvalue of η of order N , as in the bosonic case.
Using DMRG, we have calculated the eigenvalues of the matrices ρ ij and η ij for all the ground states in different phases in Fig. 1. To illustrate the Bose condensation, we have shown our results in Fig. 4(a)-(d) as one increases the detuning from −4 to 2 at n = 0.75. This path takes one from the F P P phase to the SM region and then to the W M region. We see from Fig. 4(a)-(c) that both F P P and SM phases have a distinct large eigenvalue, whereas the large eigenvalue gradually disappears into a continuous (power-law like) distribution of eigenvalues in the W M region, see Fig. 4(d). In contrast, the fermion pair correlation η does not show a large distinct eigenvalue in all cases, and appears to be power-law like. The situation at the F HP phase is similar to that of the F P P phase, i.e. there is a distinct eigenvalue for the boson density matrix, and a power-law like distribution for the fermion pair distribution.
(IV) Edge-edge Majorana correlation. The emergence of Majorana fermion is associated with long range edgeedge correlations. In the Kitaev chain model [4], one defines two sets of Majorana fermion operators , and the Majorana phase can be characterized by the order parameter O ≡ i λ 1 1  directly manifests the correlation between two Majorana modes at different edges [4,[19][20][21]23]. For a number conserving system, |O| is directly reduced to the edgeedge correlation function G (1, N L ), where G is defined as: In Fig. 5, we show the behavior of correlation function G(1, j) as the system passes through different boundaries marked (a) to (d) in Fig. 1. For filling n > 1, as the system enters the SM phase through the phase boundary ν c1 from the F P P phase below, or through the phase boundary ν c2 from the F HP phase above, G(1, j) immediately shows strong revival as j approaches the other edge, which is the hallmark of Majorana edge states. See Fig. 5(c) and Fig. 5(d). In contrast, in the F P P phase (ν < ν c1 ) or F HP phase ( ν > ν c2 ), G(1, j) decreases exponentially fast as j increases, indicating the absence of Majorana edge correlations.
For n < 1, similar revival behavior also shows up as ν across the lower boundary ν c1 (see Fig. 5(a)) . Here, within a small range of ν from −3.5 to −2.5, the edge-edge correlation G(1, N L ) (scaled by G(1, 1)) increases from 10 −3 to as large as 0.6. Continuously increasing ν, G(1, j) crossovers to a different behavior. It decreases slowly without strong revival but reaching a small yet non-zero value as j approaches the other edge, see Fig. 5(b). Specifically, for a large range of ν from 0 to 4, G(1, N L ) decreases from 0.2 to 0.05. Further increasing ν, G(1, N L ) becomes even smaller but still finite (not exponentially small as in the non-Majorana phases F P P and F HP ). In this sense, in Fig. 1 we draw a dashed line for filling n < 1 to show the crossover from (V) Momentum distribution: The behavior of correlation function in spatial space (∼ f † i f j ) directly determines its Fourier transformation, i.e., the momentum distribution of fermions n(k) = f † k f k that can be easily measured in cold atoms experiments. In Fig. 6, we show n(k) from DMRG calculation for three typical values of ν at filling n = 0.75 and 1.25. We can see that in both cases, when ν stays in the SM region(red dashed line in Fig. 6), n(k) features a distinct peak at k = 0 while decays to zero at the Brillouin edge k = π/d. We have checked that such property of n(k) holds true for all SM states in Fig. 1. In contrast, n(k) is roughly a constant in the F HP phase, and has a hole at k = 0 in the F P P phase. n(k) for WM state share similar structure as the SM case, while its peak at k = 0 is not as sharp as the latter.
All these behaviors can be understood from the meanfield theory, where we have the analytical expression n(k) = [1 + (2t f cos(kd) + µ)/E k ]/(2N L ) (see Eq. (6)). It is easily seen that within the Majorana region (|µ| < 2t f ), n(k) is the largest at k = 0 while it is zero at k = ±π/d. In the F HP phase, the detuning ν is large and positive, which leads to a large and positive µ(> 2t f ) in the mean field theory and a finite n(k) at both k = 0 and k = ±π/d. Such detuning makes it costly to create bosons, making the system more free fermion like. In the F P P phase, ν is large and negative, forcing most of the particles into the k = 0 Bose condensate, and leads to a large and negative µ(< −2t f ) in mean field theory. This leads to n(k) = 0 at both k = 0 and k = ±π/d. In particular, at small k, we have n(k) = |∆ sin(kd)/(2µ)| 2 ∼ k 2 .
Given the distinct n(k) for different phases, they can be used as experimental indicators of various phases in current system. Note, however, that n(k) cannot be directly related to the edge mode in Majorana physics, because it contains a large contribution from the bulk. This is evidenced by the fact that the behaviors of n(k) are well accounted for by the results of mean field theory. In this section, we will show that the strong Majorana (SM) phase has robust topological features, in that the edge-edge correlation survives from disorders, and it can host a non-trivial winding number in the bulk system. We will also show that the condensation of bosons and the edge-edge correlation remain strong for increasing system size.
Firstly, we study the robustness of edge-edge correlation in the presence of disorder. Here we impose disorder either in the fermion hopping term (t f ) or in the coupling  (1), and carry out DMRG simulations with OBC. Specifically, t or g are now sitedependent, t f → t f + Dδ j or g → g + Dδ j (j is site index); here δ j ∈ (−1, 1] is a random number, and D is the strength of disorder. In Fig. 7, we show the behavior of correlation function G(1, j) for a typical SM ground state with different disorders in t f [ Fig. 7(a)] and g [ Fig. 7(b)]. We see that small disorder cannot change the strong revival character of G(1, j) even for D reaching 0.25. This shows the edge modes are robust against a fairly large amount of external perturbations. Secondly, in order to demonstrate the topological feature of the bulk system, we impose a twisted phase boundary condition and calculated the resulting winding number. Specifically, we turn on the fermion hopping between the two edges (i = 1 and L) as t * 1L = t L1 = t f e iθ , where θ ∈ (0, 2π] is the twist phase. In the single-particle picture, this corresponds to shifting the momentum as k → k + θ/N , and when θ varies from 0 to 2π, the momentum basis {k} returns to itself and completes a closed loop, so as the Hamiltonian H. For the interacting manybody state, we calculate the winding number of the form

term (g) in the Hamiltonian
with Ψ θ is the ground state with twisted phase θ. In Fig. 8 we show w as a function of detuning ν for several fillings by exactly diagonalizing small size systems. We find that given the filling factor n, by increasing the detuning ν to drive the system from FPP to SM phase, w will have a sudden jump from 0 to a finite value (= nπ) at a critical ν c (in thermodynamic limit ν c is expected to recover the lower boundary as shown in Fig. 1). This signifies a topological transition between the two phases. The finite w continues to the WM phase when further increasing ν.
Note that here w depends on the filling factor, instead of a constant (π) Berry phase in the mean-field analysis (see section III). This difference can be attributed to different ways in introducing a closed path in parameter space. Specifically, in the mean-field analysis, the closed path is completed by moving k through the entire Brillouin Zone, which in the single-particle picture corresponds to fermions occupying a Fermi-sea at full filling n = 1. Here, for interacting many-body system, the closed path is introduced through the twisted boundary and the filling n can be arbitrary. Nevertheless, a remarkable feature of the Majorana phase is that, regardless of the way of introducing closed path, it can always distinguish itself from the trivial phase by producing a non-zero (topological) winding number. Such a non-zero number characterizes the topological nature of the bulk for interacting many-body systems, analogous to the role of π Berry phase in the Kitaev chain under mean-field treatment.
Finally, we study the robustness of the Bose condensation and the edge-edge correlation against increasing the system size. We diagonalize the boson single-particle density matrix ρ ij for a typical SM state with different N L = 24, 40, and 64 in the SM phase. Fig. 9(a) shows that in all three cases, ρ has a distinct largest eigenvalue. We also show the edge-edge correlation function G(1, N L ) for different N L in Fig. 9(b), and the same strong revival at the edge is found in for all lengths studied.

C. Comparison with the single-channel model
It is useful to contrast the Hamiltonian (1) with the single-channel fermion model, (12) where N f,i = f † i f i is the fermion number operator at site i and U < 0 is the attraction between neighboringsite fermions. In this model, the fermion number N f = i N f,i is conserved, unlike in the resonance model (H in Eq. (1)). Yet this model has the same mean field theory as the resonance model, with the mean field gap defined as ∆/2 ≡ U f i f i+1 . This raises the question of whether the single channel model H sc will also have a Majorana ground state in certain parameter regime.
To compare the ground state of single-channel model H sc [Eq. (12)] with that of the resonance model H [Eq. (1))], we shall choose the parameters ({U, n f } in H sc and {ν, n} in H) such that both models have the same fermion density n f and the same mean field gap ∆. With this correspondence, we have calculated the ground state of H sc with OBC using DMRG.
We find that for all the detunings ν in Fig. 1 that cover the SM and F P P phases, the corresponding U in H sc is so negative that the ground state is a droplet with all fermions packed together in a region of the size of N f sites, see Fig. 10(b). Such cluster bound state was also shown previously for few particles [38]. It can be understood by mapping H sc into a quantum spin chain using the Jordan-Wigner transformation, where the occupied (empty) site is mapped to spin-up (spin-down), and the U (< 0) term in Eq. (12) can be mapped to the ferromagnetic Ising interaction. It is then obvious that for large and negative U , the system forms ferromagnetic domains in the ground state, i.e., the occupied and empty sites are spatially well separated as shown in Fig. 10(b). Similar ferromagnetic correlation has also been shown in other 1D systems with p-wave attraction [39,40]. In our calculations, the droplet may appear in different locations, as the energy difference between droplets at different locations is so small that is below our accuracy of our calculation. Clearly, the droplet phase is not the Majorana phase as found in the resonance model H, which exhibits the fermion density distribution as shown in Fig. 10(a).
For weaker attraction U , corresponding to W M or F HP regions in the large detuning limit in Fig. 1, the droplet gives way to a gas phase that covers the entire chain, but still there is no strong edge-edge correlations. To conclude, the single-channel model H sc cannot host strong Majorana character for all couplings U , distinct from the SM phase in Fig. 1 of the resonance model. It is also clear from Fig. 1 that in order to obtain strong Majorana correlations, the detuning ν should stay in a finite region near the two-channel resonance (ν ∼ 0), i.e., when bosons and fermions have comparable proportions and their conversion (or number fluctuation of fermion pairs) is the strongest.
Here we should also remark that to describe p-wave Fermi gas in cold atoms experiments, the two-channel model is more realistic than the single channel model. This is because the p-wave resonance in these systems are generally very narrow, and the closed-channel bosons can take a sizable proportion as measured in a 3D gas near a p-wave resonance [29].

D. Suppressed atom loss in the strong Majorana phase
Experimental realization and detection of Majorana edge state require low atom loss. For p-wave fermions in the lattice configuration, a previous study showed that the lattice setup will help to reduce inelastic collisional losses compared to free space [41]. The analysis was based on a single-channel model, and the reduced loss can be attributed to the low probability of finding three fermions close to each other outside the lattice sites [41]. For the present p-wave system described by the twochannel lowest-band model, three-fermion collision can be effectively ruled out, while the atom loss is dominantly caused by the fermion-boson or boson-boson collision at the same site. Indeed, previous studies on a continuum gas have shown that the three-body and the four-body loss rates are respectively proportional to the probabili-  ties of finding atom-dimer and dimer-dimer at the same location [42,43], up to a background constant that is determined by the loss rate far from resonance regime.
Here, accordingly we examine the probabilities of finding a pair of boson-fermion and boson-boson at the same site, respectively denoted by P bf and P bb : with Fig. 11, we show P bf and P bb as functions of filling n for ν staying in the lower (ν c1 ) and upper (ν c2 ) boundaries of the SM phase in Fig. 1. We can see that for n 1.25, both probabilities are less than 10%, suggesting the atom loss is well controlled with little atom-dimer and dimer-dimer collisions. The physical reason for these low probabilities is because these configurations do not effectively take advantage of the conversion between boson and fermions (g-term in Eq. 1) to lower the energy. For example, a boson and a fermion on the same site will stop the boson to convert into a fermion pair due to Pauli blocking. Similarly, if two bosons are at the same site, they cannot both convert to fermion pairs. As a result, in general the ground state does not favor the double occupations of boson-fermion or boson-boson at the same site. However, for large fillings, such double occupations are inevitable, as shown by the increasing P bf and P bb with n in Fig. 11. These suggest that the SM phase in Fig. 1 should be stable enough for lower fillings (n 1.25). Now we give an estimation to the loss rate of 40 K and 6 Li fermions in 1D lattices. For 40 K and 6 Li away from Feshbach resonance, the 3D recombination rates are α 3D rec = 10 −25 cm 6 /s [26] and 10 −24 cm 6 /s [27], respectively. For typical transverse confinement length ∼ 50nm and typical 1D density ∼ 10 4 cm −1 , this leads to the decay time about 1s for 40 K [40] and 0.1s for 6 Li when the system is out of the resonance regime (non-interacting limit). In the present case, due to the small probability P bf 10% (for filling less than unity), the actual loss rate will be further reduced by one order of magnitude, i.e, the decay time can extend to 10s and 1s, respectively, for 40 K and 6 Li systems. Considering the typical hopping strength t f about tens to hundreds of Hertz, the time scale for developing the many-body correlation is a few to tens of milliseconds, much shorter than the decay time. We thus expect the Majorana phase can be observed well before severe losses occur in practical cold atoms experiment.

E. Effect of inter-channel coupling
The phase diagram shown in Fig. 1 is for coupling g = 1 (in units of hopping t f ). To illustrate the situation for different g, we have worked out the phase diagrams in the (g, ν)-plane for two different fillings n = 0.75 and 1.25 using DMRG. These results together with the mean field predictions are shown in Fig. 12.
From Fig. 12, one sees that the difference between DMRG and mean field results grows with increasing g, and in large g limit the mean-field theory significantly overestimate the SM region (marked by yellow color) compared to DMRG result. This can be attributed to the enhanced quantum fluctuations as increasing g. For filling n < 1 (Fig. 12(a)), the SM phase can always survive in a finite detuning regime, while the lower and upper boundaries both turn upward to higher detunings. In comparison, for filling n > 1 ( Fig. 12(b)), the SM phase finally disappears at a large g c (g c = 5 for n = 1.25). For g > g c the DMRG result suggests that Majorana physics is overwhelmed by certain density waves of bosons and fermions in lattices.
Now we show that it is realistic in experiments to reach the parameter regions (g, ν) of the Majorana phase. To give an example, we shall consider the 1D 40 K Fermi gas in a lattice with depth v ≡ V 0 /E L = 6 (V 0 is the lattice depth and E L = k 2 L /(2m) is the recoil energy). The hopping t f in the lowest band is t f /E L ∼ 0.06. As shown in Ref. [24], (g, ν) can be expressed by g = g ef f Cd −3/2 , ν = −2g 2 /U ef f , where g ef f is related to the effective range r ef f = (mg ef f ) −2 , U ef f is the effective coupling between fermions (see Eqs. (11) and (13) in Ref. [24]), and C is a constant given by the overlap of Wannier functions (C = 0.06 for v = 6, see Fig. 4 in Ref. [24]). Let us consider the regime nearby the first Bloch-wave resonance, where (l o k L ) −1 ≤ 2 (l o is the odd-wave scattering length and k L = π/d is the recoil energy). The range of U ef f and r ef f are shown in Fig. 3 in Ref. [24], from which one can estimate the range of (g, ν) in unit of t f . For instance, in the interaction regime of interest, r ef f can range from 0.75d to 1.5d, so the ratio between the coupling g and the hopping t f can range from 0.15 to 0.25. Similarly, from the information of U ef f one can estimate the range of ν/t f as from −7 to 42. Such a broad range of ν/t f is facilitated by the small value of t f /E L , and it well covers the SM region shown in Fig. 12 for g/t f ∈ [0.15, 0.25]. Therefore the strong Majorana correlation can be achieved in a lattice with v = 6 near a Bloch-wave resonance.

V. SUMMARY AND DISCUSSION
In summary, we have shown that the spinless Fermi gas in a 1D optical lattice near a p-wave resonance can have Majorana ground state over a sizable range of parameter space that are experimentally accessible. Our scheme makes use of the intrinsic property of cold atoms with double channels and requires neither spin-orbit coupling nor multi-chain setup. Our work, together with other multi-chain studies, show the number fluctuation of fermion pairs are crucial for the formation of Majorana phase. In comparison, we demonstrate that the singlechannel fermions with neighboring-site attraction have no strong Majorana features.
In identifying the phase boundaries between the Majorana phase and other trivial superfluid phases, we have examined a number of different physical quantities, including the entanglement entropy, the boson fraction and edge-edge correlation, which give rise to consistent results as shown in Fig.1. In the practical detection of Majorana phase, the low probability of dimer-fermion and dimerdimer pair at the same site will help to reduce atom loss. In addition, it is proposed to identify various phases in the present system from the momentum distribution of fermions, which can be easily measured experimentally.
Our results can be directly tested in the 1D cold atomic gases of 40 K or 6 Li fermions.
Finally, we further summarize our characterization of Majorana edge state in interacting many-body systems. In this work, we show that the phases that we labeled to be Majorana (SM phase) exhibit the following properties identical to the (number non-conserving) Kiteav chain: (i) The ground state in an open chain exhibits a strong edge-edge correlation that is robust against various kind of disorder. (iii) The corresponding ground state in the bulk has a non-zero winding number, distinguished from the nearby phases which has zero winding number. (iii) The phase diagram of our number-conserving Majorana phase is remarkably closed to that of the mean-field theory, which is the Kitaev model (see Fig. 1). (iv) The properties of our Majorana state can also be interpreted from the Bosonalization method, which has been carried out in Ref. [44] for a similar boson-fermion model. A second-order topological transition was found, consistent with our findings as shown in Fig. 1.
All above evidences (i)-(iv) show that much of the essential physics of Majorana state exhibited in the number non-conserving Kitaev chain also appear in our number conserving model. In other studies of Majorana physics in number-conserving models [19][20][21]23], a ground state degeneracy between different number parity sectors has been established. Our system corresponds to one of the fixed number parity states (i.e. either odd or even fermion number) and we have focused on the Majorana correlation function. We shall explore that the physical effects related to the long range coherence of the Majorana correlation in future studies.