Strongly interacting Bose-Fermi mixtures in one dimension

We study one-dimensional strongly interacting Bose-Fermi mixtures by both the exact Bethe-ansatz method and variational perturbation theory within the degenerate ground state subspace of the system in the infinitely repulsive limit. Based on the exact solution of the one-dimensional Bose-Fermi gas with equal boson-boson and boson-fermion interaction strengths, we demonstrate that the ground state energy is degenerate for different Bose-Fermi configurations and the degeneracy is lifted when the interaction deviates the infinitely interacting limit. We then show that the ground properties in the strongly interacting regime can be well characterized by using the variational perturbation method within the degenerate ground state subspace, which can be applied to deal with more general cases with anisotropic interactions and in external traps. Our results indicate that the total ground-state density profile in the strongly repulsive regime behaves like the polarized noninteracting fermions, whereas the density distributions of bosons and fermions display different properties for different Bose-Fermi configurations and are sensitive to the anisotropy of interactions.


I. INTRODUCTION
One-dimensional (1D) quantum gases have attracted renewed attention during the past decades due to the experimental progress in trapping and manipulating cold atomic systems [1,2]. By loading bosonic or (and) fermioic cold atoms in 1D waveguides, one may realize bosonic or fermionic gases (the Bose-Fermi mixtures). In comparison with the bosonic and fermionc gases, the Bose-Fermi mixtures are particularly interesting as they rarely occur in nature but are accessible in current cold atomic experiments [3][4][5]. Theoretical investigations have unveiled the Bose-Fermi mixtures displaying rich phase diagrams and interesting excitation properties [6][7][8][9][10][11][12][13][14]. The Bose-Fermi mixtures also provide a platform to realize and study physical properties of the Bose-Fermi supersymmetry [15,16].
The tunability of interaction strengths between ultracold atoms has provided unprecedented opportunities for investigating intriguing many-body physics in 1D quantum gases in the entire parameter regime [17][18][19][20][21], with the symbolic experimental progresses in the realization of Tonks-Girardeau (TG) gas [20,21] and super-Tonks-Girardeau (STG) gas [22]. Furthermore, recent experiments on few-particle atomic systems with the controllability of precise atom numbers open access to studying few-body physics and the size-dependent evolution from few-body to manybody systems [23,24]. Besides the bosonic TG gases, the fermionization of the interacting fermion system has also been observed in the few-particle fermion system [25].
It is well known that the TG gas corresponds to a 1D Bose gas with infinitely repulsive interactions between bosonic particles, which was well understood by using the Bose-Fermi mapping proposed by Girardeau in his seminal work more than fifty years ago [26]. Motivated by the cold atomic experimental progress, the Bose-Fermi mapping was also generalized to study 1D multi-component quantum gases in the strongly repulsive limit [27][28][29][30]. Different from the single-component bosonic TG gas, the ground state of multi-component quantum gases in the infinitely repulsive limit is highly degenerate with the degree of degeneracy given by the number of distinct species (spin) configurations [27][28][29][30][31][32][33]. Slightly away from the infinitely repulsive limit, a perturbation theory within the degenerate subspace can be constructed using with the inverse of the interaction strength as a small parameter [34][35][36][37][38][39]. Particularly, the effective spin-exchange Hamiltonian, which describes the spin dynamics of the spin-1/2 boson and fermion systems in the strongly repulsive regime, has been derived [36][37][38][39][40]. The trapped multicomponent systems in the full interaction regime have also been studied by exact solutions [41,42], highly accurate numerical diagonalization methods [43][44][45][46][47] and variational methods [48][49][50][51][52].
While most of the previous works concentrated on the strongly interacting bosonic and fermioic systems [27][28][29][31][32][33][34][35][36][37][38][39][40][41][42][43][44][45][46], the Bose-Fermi mixtures in the strongly interacting regime are not well studied except for the limiting case with infinite repulsion [27,30]. It is still not clear how the degenerate ground states split when the interaction deviates the infinite repulsion limit. Another interesting question is how the anisotropy of boson-boson interaction and bosonfermion interaction affects the physical properties of the strongly interacting mixture system? Aiming to answering the above questions, in this work we shall focus our study on the ground state properties of the 1D Bose-Fermi mixtures in the strongly repulsive regime. In order to get some exact results and provide a touchstone for the following results based on the perturbation theory within the degenerate subspace, we first consider the isotropic case with equal boson-boson and boson-fermion interactions, which is exactly solvable by the Bethe-ansatz method in the absence of external potentials. Next we derive the universal energy relation for the Bose-Fermi mixtures and then present our perturbation theory within the degenerate subspace. By comparing the variational result with the exact result, we find that they agree very well in the strongly repulsive regime. The variational perturbation theory is then applied to deal with the anisotropic mixtures with different boson-boson and boson-fermion interaction strengths trapped in a harmonic trap. Our results indicate that the anisotropic parameter has a significant effect on the ground state density distribution of the Bose and Fermi components, although it almost does not affect the total density distribution of the mixtures in the strongly repulsive regime.

II. MODEL AND EXACT SOLUTIONS FOR ISOTROPIC MIXTURES
We consider the 1D interacting Bose-Fermi mixtures described by the Hamiltonian Here, Ψ b and Ψ f denote the bosonic and fermionic annihilation operators, respectively. The mixtures are confined in external traps V b (x) and V f (x). The boson and boson or boson and fermion are interacting through the contact interaction with different strengths g bb and g bf (define anisotropy η = g bf /g bb ), while the interaction between fermions is forbidden by the Pauli exclusion principle. The model (1) in the absence of external potential, i.e., with 10,53]. While the first condition is approximately fulfilled for isotopes of atoms (for example 171 Yb and 172 Yb, and 86 Rb and 87 Rb), the second condition can be realized by tuning the interaction strength via Feshbach resonances. In this work, we only consider the case with m b = m f ≡ m, but relax the restriction g bb = g bf when we discuss the anisotropic case with η = 1 by using the variational perturbation theory. Few-body systems with different inter-component and intra-component interaction strengths have been studied in some recent works [48][49][50][51][52].
Next we shall consider the exactly solvable model with equal masses and equal repulsive boson-boson and bosonfermion interaction strengths. To keep consistent with the traditional references for the integrable Bose-Fermi model, we set c = mg/ 2 . Consider the Hilbert space spanned by N particles, then the eigenvalue problem reduces to solve the Schrödinger equation with first quantized Hamiltonian Among the N particles, there are N b bosons and the rest of them are fermions. Under periodic boundary condition, the Bethe-ansatz equations (BAEs) are given by where M = N b . We have set = 2m = 1 and θ(x) = tan −1 x c/2 . The quantum numbers I j and J α are integers or half integers, depending on the parity of N b and N . In general, k j s are called as quasi-momenta, and Λ α s are called as rapidities. For c > 0, both k j s and Λ α s take real solutions. The eigenenergies are given by E = j rewritten in order of k j /c by using the Taylor expansion where we have used tan −1 (−x) = − tan −1 (x) and d(tan −1 (x))/dx = 1/(1 + x 2 ). Keeping to the first order of k j /c, the quasi-momenta can be calculated: where Here we have used relations j k j = 0 and α θ(Λ α ) = 0, which are true if I j and J α are symmetrical and are always fulfilled in the case of ground state. By using 2πJ α = 2N θ(Λ α ), which is just the second BAE under the first order Taylor expansion, we have ζ ≈ M α=1 4 tan 2 (πJα/N )+1 . In thermodynamical limit, the sum can be calculated via integration and we have ζ = 2M + 2N π sin( Mπ N ). The energy of the mixture gas in the strongly repulsive limit is thus given by In the limit c → ∞, the ground energy is identical to that of a polarized N-fermion system and thus is degenerate for different Bose-Fermi configurations. When c deviates the infinitely repulsive limit, the degeneracy is lifted as the value of ζ is dependent on the number of bosons. We note that ζ( for systems with fixed total particle numbers. For attractive interactions with c < 0, the BAEs can still have real solutions, which describe the scattering states of the attractive systems. The lowest states with real solutions in the strongly attractive regime correspond to the STG state, similar to STG states in the attractive Bose systems [54][55][56][57] and Fermi systems [31]. We have In The degeneracy is also lifted when c deviates the infinitely attractive limit, and we have Fig.1(a), we demonstrate the ground state energies versus the inversion of interaction strength 1/c for Bose-Fermi mixtures with total particle number N = 4. Totally there are five cases with different Bose-Fermi configurations: 4B, 3B1F , 2B2F , 1B3F and 4F , where the configuration nBmF represents the system composed of n bosons and m fermions. These different states are degenerate at the exact Girardeau's mapping point The first two cases make no difference since no Fermi statistic should be considered. We can clearly see that the degeneracy of the ground state energy is lifted when the system deviates the exact mapping point. At the repulsive side, i.e., in the TG regime, the ground state energy decreases with the decrease of c. Systems with more bosons have relative lower energies at this regime. While at the attractive side, i.e., in the STG regime, the energy increases with the decrease of |c|, and systems with more bosons have relative higher energies. To the linear order, they are perfectly described by the above expansion formula from the BAEs. To be more intuitive, Fig.1(b) shows the distribution of quasi-momenta with respect to the interaction strength for the 2B2F system in the whole repulsive regime. At the non-interacting limit, c = 0, the quasi-momenta should be 0 for bosons and ± π L for fermions. As the interaction increases, the quasi-momenta distribution becomes wider and finally at c = ∞, we have k j = 2πI j /L. The whole system behaves like a polarized fermion system with quasi-momenta ± 3π L , ± π L , which is consistent with the result by the generalized Bose-Fermi mapping [27].

III. UNIVERSAL ENERGY RELATION FOR THE 1D BOSE-FERMI MIXTURE
Next we shall derive the universal energy relation for the 1D Bose-Fermi mixtures, which can be viewed as a generalization of Tan's universal energy relation for the spin-1/2 Fermi gas [34,[58][59][60]. For convenience, we only consider the case with m a = m b = m and V a (x) = V b (x) = V (x) and the Hamiltonian can be rewritten as the following form where σ i = b, f represent the bosonic and fermionic components, respectively, δ b σi,σj = 1 only if σ i =σ j =b and δ σi,−σj = 1 when ( Here the summation to the spin index is assumed. Let Ψ(x 1 , σ 1 ; x 2 , σ 2 ; ...; x N , σ N ) be the normalized eigenstate of the system, which fulfills the Schrodinger equation: Denote x 1 , x 2 as the coordinates of two interacting particle with interaction g bb or g bf , depending on whether they are two bosons or one boson and one fermion. Suppose that x 3 , .., x N do not coincide with x 1 and x 2 . In terms of center-of-mass and relative coordinates, The equation can be rewritten as This is just the boundary condition required for two interacting particles. At small r, we have the Taylor expansion [34]: where X = (x 3 , σ 3 ; ...x N , σ N ), a bb = − 2 2 mg bb and a bf = − 2 2 mg bf are the effective 1D scattering lengths, and Similarly, for a different coupling g ′ with energy E ′ , we have From Integrating this equation over the coordinates and summing over all the interacting pairs, we have where dX = dx 3 ...dx N , N bb P = N b (N b − 1)/2 for boson-boson interactions and N bf P = N b N f for boson-fermion interactions. Taking the limit of g ′ → g and Ψ ′ → Ψ, we get where I bb and I bf are the contacts defined by From Eq. (18), we have The above relations are the universal energy relations of the Bose-Fermi mixtures.
For the Bose-Fermi mixtures with fixed anisotropy parameter η, we can get the expression for the energy near the infinitely interacting limit where E ∞ represents the energy of the system at the TG limit. Some results directly follow from this formula. First, since I bb ∝ N b (N b − 1)/2 and I bf ∝ N b N f , for the isotropic interacting case with η = 1 and fixed particle numbers N , we can conclude that systems with more bosons have lower energies. When the boson-boson interaction dominates, η ≪ 1, the second term in the contact matrix requires the system to have the lowest energy at N b = N f . When the boson-fermion interaction dominates, η ≫ 1, the first term requires the state with more bosons having lower energy. Generally, for η < 1, the states with the relative ratio N b /N ∝ 1 2−η have the lowest energy. Our universal relations derived in this section coincide with the exact results given by the BAEs.

IV. VARIATIONAL PERTURBATION METHOD
The Bethe-ansatz method is powerful but only limited to the integrable case with g bb = g bf and V (x) = 0. While for a trapped system or the system with anisotropic interactions g bb = g bf , we need develop a more general method to address the problem. We notice that there exists one exact soluble point 1/g = 0, i.e., the infinitely repulsive limit g bb = g bf = ∞, in which the system can be mapped to the polarized fermionic system and the many-body wave function can be constructed from the single-particle fermionic wave function, taking into account the statistics of exchange between bosons or fermions [27,30]. This is the central idea of the generalization of Girardeau's Bose-Fermi mapping to the multicomponent systems [27][28][29][30]. Since no symmetry is required for the exchange between bosonic and fermionc particles, there exists degeneracy for the eigenenergy. For the system composed of N b bosons and N f spinless fermions, the degeneracy is D = which corresponds to different configurations of N b bosons in N single-particle states.
Once the system deviates from the infinite repulsion limit, the degeneracy of ground state manifold is lifted. As long as the system is still in the strongly interacting regime with the inverse of interaction strengths much smaller than the single particle level spacing, we can utilize the degenerative perturbation to determine energy splitting. For our system, at infinite interaction strength, the ground state is D-fold degenerate and the energy can be expressed as E = N l=1 ǫ j l where j l s are N different integers and ǫ j l is the j l -th single particle energy level in the trapping potential V (x) with corresponding wave function φ j l (x). In this paper, we focus on the ground state where j l = l. We should stress that our method can also be applied to excited states which depend on the j l we choose. By the generalized Bose-Fermi mapping [27][28][29][30], the coordinate part of the many-body wave function can be constructed from the anti-symmetric Slater determinant where P is a permutation of the integers (1, 2, ...N ) and sgn(P ) = ±1 for even and odd permutations. For the infinite interaction strengths the particles can not penetrate each other and the real space can be divided into N ! distinct parts. Introducing where α is a sequence of [1, 2, ...N ]. Totally we can construct N ! orthogonal basis ψ A θ α . Once we consider constrains of the Bose-Fermi statistics, that is, exchange between fermions should contribute a minus sign while no sign will emerge for the exchange between bosons, the number of allowed eigenfunctions is reduced to D = which gives the degeneracy of the considered state. So these constructed states can be used as the basis of the degenerate space. Define permutation operators for bosons and fermions as P b and P f , respectively, we can get the D normalized and orthogonal basis of the degenerate subspace as follows: where (−1) P b = ±1 for even and odd permutations between bosons. For the case with a finite large interaction, the eigenfunction should approach the infinite one smoothly when g goes to infinity and it therefore should go into the degenerate subspace. Introducing the projection operator: P deg = α |ψ α ψ α |, then it is reasonable to expand the eigenfunction at finite interaction as with α |a 2 α | = 1. The Bose-Bose contact I bb and Bose-Fermi contact I bf can be calculated using the reduced contact matrices J bb and J bf as: where − → a = (a 1 , a 2 ...a D ) T and the reduced matrices for Bose-Bose and Bose-Fermi interaction are defined as Then the energy can be read as: Finally, we define total contact as J = J bb + η −1 J bf . This is the key quantity which determines the behavior of our system. We can clearly see that the quantum state is largely dependent on the ratio η of the two interaction strengths. The contact can be determined via the variational principle [34,35]. Let L = E − λ( α a * α a α − 1), the variational principle requires: From the above equation we can see λ and α are eigenvalue and eigen-vector of the total D × D contact matrix J [34]. The diagonalization procedure gives the splitting energy of original degenerate manifold: To show the validity and power of the method, we compare the results with expansions given by BAEs for g bb = g bf in uniform space with periodic boundary condition. As 1/g = 0, k j = 2π/I j , the totally anti-symmetric wave- ]. Take N = 3 and N = 4 cases as examples. The eigenvalues for the ground states given by the BAEs to the first order Taylor expansions are 48 4 π 2 m 2 L 3 and 160 4 π 2 m 2 L 3 , respectively. For the variational perturbation method, the largest eigenvalue by the diagonalization of the contact matrix produces the same result for the same system.

V. APPLICATION TO FEW-PARTICLE SYSTEMS IN A HARMONIC TRAP
In the above section, we have generally discussed the variational perturbation method and introduced the reduced contact matrix to simplify the strongly interacting problem to the matrix diagonalization. In this section, we apply this method to study few-particle Bose-Fermi mixtures with N = 3 and N = 4 in a 1D harmonic trap, i.e., V (x) = 1 2 mω 2 x 2 . As the trap potential preserves the inversion symmetry, the eigenstates are characterized by distinct parities. Our discussion is valid for finite but strong interaction strengths. The interaction g is larger than any other scale in the problem, i.e., g/ ωa ω ≫ 1. The single-particle state can be expressed as φ i (x) = is the Hermite polynomial and ξ = x/a ω ≡ x/ /mω. For the N-particle system, the Slater determinant composed of the N lowest eigen-functions is ∆ = [61]. We first take the 2B1F case as an example. At the infinite repulsion limit, the ground state is 3-fold degenerate. Denote x 1 , x 2 and x 3 as coordinates of bosons and fermions, respectively. We can define the following three distinct subspace basis respecting exchange statistics: An explicit calculation gives the contact matrices for the boson-boson interaction and boson-fermion interaction as By solving the eigen-equation (31) for the contact matrix, we can directly get the variational wavefunctions and energies. We note that the wavefunctions obtained in our scheme automatically fulfill the symmetry of parity. For example, for the isotropic case with η = 1, we get three eigenvalues of total contact matrix: . The corresponding normalized eigenvectors (a 1 , a 2 , a 3 ) T are: 1 It is clear that the ground state and the second excited state have even parity while the first excited state has odd parity as the even parity require a 1 = a 3 and the odd parity state requires a 2 = 0 and a 1 = −a 3 . Since the interaction terms do not change the parity of the eigenstate, the parity of the many-body state does not change even we tune the interaction strength continuously to the non-interacting limit, and thus the three many-body states can be adiabatically connected to their noninteracting limits [35], labeled by (3, 0, 0),(2, 1, 0) and (1, 2, 0), respectively. Here (n 1 , n 2 , n 3 ) means that the occupation numbers on the three lowest single-particle states are n 1 ,n 2 and n 3 , respectively.
The density profiles can be calculated from the reduced sing-particle density matrices as: where X ′ = (x 2 , ..., x N ) and X ′′ = (x 1 , ...x N −1 ). The diagonal elements are nothing but the single-particle density profiles, i.e., ρ b (x) = n b (x, x) and ρ f (x) = n f (x, x). In Fig.2 (a)-(c), we show the density distribution for the 2B1F system with different interaction anisotropy η. For η ≪ 1, the boson-boson interaction dominates. As expected, the bosons will be repelled to the two wings while the single fermion mainly locates in the middle regime. For η ≫ 1, the boson-fermion interaction dominates, the overlap between bosons and fermion must be the smallest to avoid the strong inter-species interaction. The bosons can locate on the middle regime and the fermion must be repelled from the center to form two peaks. Particularly, for the case of η = 1, the two interactions equally compete and we have n b = 2 3 n G and n f = 1 3 n G , where n G = N i |φ i (x)| 2 denotes the density distribution in the TG limit. On the other hand, for the 1B2F case, the fermions are repelled to the wings with the bosons located at the trap center as shown in Fig.2 (d). In all cases, the totally density distribution is nearly as the same as n G .
Next, we consider the equal-mixing mixtures composed of 2 bosons and 2 fermions in the harmonic trap, i.e., the 2B2F case. At the infinitely repulsive limit, the ground state is 6-fold degenerate. Denote x 1 , x 2 and x 3 , x 4 as coordinates of bosons and fermions, respectively. We can define the six distinct subspace basis as follows according to exchange statistics: Here C 1 = 0.5938 and C 2 = 0.7796 are integral constants. The parity of the state can be directly read out from the eigenvectors of the contact matrix. For the even parity state, a 1 = −a 6 , a 2 = −a 5 and a 3 = a 4 = 0 while for the odd parity state the coefficient should satisfy: a 1 = a 6 , a 2 = a 5 . The diagonalization of the total contact matrix demonstrates that there always exist four odd and two even parity states in the six degenerate subspaces. Similar to 2B1F case, we can label these states by their adiabatically connections with the non-interacting states. From the ground state to the fifth excited state, the corresponding occupations in the single-particle orbital are (3,1,0,0), (2,2,0,0), (1,3,0,0), (2,1,1,0), (1,2,1,0) and (1,1,2,0), respectively. The density profiles are shown in Fig.3 for different interaction ratio η, exhibiting a bit of difference from the 2B1F case. First, when η ≪ 1, the boson-boson interaction dominates and the two bosons behave like hard-core bosons. As shown in Fig.3(a), for the 2B2F system with η = 0.1, the density distribution of bosons has almost the same distribution as that of fermions and we have n b ≈ 1 2 n G and n f ≈ 1 2 n G . As η increases, the boson-fermion interaction gradually dominates and bosons and fermions will repel each other. The fermions are repelled from the harmonic trap center while the bosons will eventually localize at the trap center for η = 10. We notice that while each component in the three different cases has quite different density distribution, the total density distribution is nearly as the same as n G .
It is interesting to indicate that at the limit g bb → ∞ but with finite g bf , corresponding to the case of η = 0, the system becomes a mixture of hard-core bosons and fermions [62]. In this case, the system can be mapped to a spin-1/2 Fermi gas by a generalized Bose-Fermi transformation: where n σ (x) = Ψ † f,σ (x) Ψ f,σ (x), and Ψ † f,σ (x) ( Ψ f,σ (x)) is the creation (annihilation) operator at location x for σ-component fermions (σ =↑, ↓). The second mapping in the above equations is introduced to enforce the fermion operators Ψ f,↑ and Ψ f,↓ fulfilling the anti-commutation relation {Ψ f,↑ , Ψ f,↓ } = 0. By this mapping, it is known that the density distributions of the Bose-Fermi mixture in the limit of g bb → ∞ are identical to the distributions of the corresponding spin-1/2 Fermi gas, and thus we have n b = n f = n tot /2 (n tot = n b + n f ) when N b = N f from the symmetry requirement of the exchange invariance for ↑ and ↓ fermions. From the above analysis, it is not hard to understand why we have n b ≈ n f for the equal-mixing system with η ≪ 1 as shown in Fig.3 (a). Also the density distributions shown in Fig.2 (a) for the 2B1F system with η ≪ 1 are consistent with the distributions of spin-1/2 Fermi gas in Ref. [29], according to the above mapping in the limit of η = 0.

VI. SUMMARY AND OUTLOOKS
In summary, we have studied the properties of 1D Bose-Fermi mixtures at the strongly repulsively limit. For the exactly solvable model with equal boson-boson and boson-fermion interactions, we give the analytical expression for the ground state energy in the strongly interacting regime, which clearly indicates that the ground state energy is dependent on the Bose-Fermi configuration and the degeneracy in the infinitely repulsive limit is lifted when the interaction deviates this limit. For the general case with different boson-boson and boson-fermion interactions, we derive the universal energy relation of the mixture system and then study the few-particle systems in harmonic traps by the variational perturbation theory within the degenerate ground state subspace of the system in the infinitely repulsive limit. Our results show that the total ground-state density profile in the strongly repulsive regime is not sensitive to the anisotropy of interactions and Bose-Fermi configurations, which however have significant effects on the density distributions of bosons and fermions. The species-dependent density distributions may be experimentally detected in a similar way as in the recent experiment for a 1D two-component fermionc system [63].
Our variational perturbation method can be directly extended to deal with larger systems with more particle numbers though the integral coefficients, which are closely related to Tan's contacts, become more complicated with the increase of particle number. After constructing the variational basis, we can investigate the crossover from fewbody to many-body systems in the same scheme described in the present paper. However, when the particle number is large, it becomes a very difficult task to determine the integral coefficients via carrying numerical multiple integral. For the large-N system, it is more convenient to study the ground state properties by developing some methods based on the density functional theory [64][65][66][67]. It is also interesting to study the perturbation correction of the ground state energy for the large-N system in a harmonic trap by generalizing the method for the single-component bosonic system in Ref. [68]. Our method can be also applied to study other multi-component systems with more inner degrees of freedom, for example spinor quantum gases with S = 1 and S = 3/2. Based on the variational perturbation method, we can also study the quantum magnetism for strongly interacting multi-component quantum gases in the future work.