One-dimensional multicomponent Fermi gas in a trap: quantum Monte Carlo study

One-dimensional world is very unusual as there is an interplay between quantum statistics and geometry, and a strong short-range repulsion between atoms mimics Fermi exclusion principle, fermionizing the system. Instead, a system with a large number of components with a single atom in each, on the opposite acquires many bosonic properties. We study the ground-state properties a multi-component Fermi gas trapped in a harmonic trap by fixed-node diffusion Monte Carlo method. We investigate how the energetic properties (energy, contact) and correlation functions (density profile and momentum distribution) evolve as the number of components is changed. It is shown that the system fermionizes in the limit of strong interactions. Analytical expression are derived in the limit of weak interactions within the local density approximation for arbitrary number of components and for one plus one particle using an exact solution.


Introduction
Quantum one-dimensional systems can be realized in ultracold gases confined in cigar-shaped traps[1, 2, 3,4,5,6]. At ultracold temperatures atoms manifest different behavior, depending if they obey Fermi-Dirac or Bose-Einstein statistics. In terms of the wave function, a different symmetry is realized with the respect to an exchange of two particles, bosons having the symmetric wave function and fermions the antisymmetric one. A peculiarity of one dimension is that reduced geometry imposes certain limitations on probing the symmetry due to an exchange. The only way of exchanging two particles on a line is to move one particle though the other. This leads to important consequences when particles interact via an infinite repulsion. In this case, known as Tonks-Girardeau limit, systems might acquire some fermionic properties. The energy of a homogeneous Bose gas with contact interaction (Lieb-Liniger model) can be exactly found [7] using Bethe ansatz method and follows a crossover from mean-field Gross-Pitaevskii gas to Tonks-Girardeau gas, in which the energetic properties are the same as for ideal Fermions [3,8,9]. The equation of state can be probed [2,4,5] by exciting the breathing mode in a trapped gas, as the frequency of collective oscillations depends on the compressibility. This allowed to observe a smooth crossover from the mean-field Gross-Pitaevskii to the Tonks-Girardeau/ideal Fermi gas value [4]. For a large number of particles, the breathing mode frequency in the crossover is well described [10] within the local density approximation (LDA) which relies on the knowledge of the homogeneous equation of state [7]. Instead for a small number of particles, the LDA approach misses the reentrant behavior in the Gross-Pitaevskii -Gaussian crossover, which first was observed experimentally [4] and later quantitatively explained [11] using quantum Monte Carlo method.
The interplay between repulsive interactions and statistics was clearly demonstrated in a few-atom experiments in Selim Jochim's group [12] where it was shown that two component fermions with a single atom in each component fermionize when the interaction between them becomes infinite. The question becomes more elaborate when the number of components becomes large. A system of N c component fermions with a single atom in each component should have a similar energetic properties as a single component Bose gas consisting of N c atoms [13]. Recently a six-component mixture of one-dimensional fermions was realized in LENS group [14]. The measured frequency of the breathing mode approaches to that of a Bose system as N c is increased from 1 to 6. It was observed that the momentum distribution increases its width as the number of components is increased, keeping the number of atoms in each spin component fixed.
The Bethe ansatz theory is well suited for finding the energetic properties in a homogeneous geometry, predicting the equation of state for the case of N c = 2 components [15,16] and an arbitrary N c [17,18]. Unfortunately Bethe ansatz method is not applicable in presence of an external potential. Using LDA is generally good for the energy but misses the two-body correlations. Instead quantum Monte Carlo methods can be efficiently used to tackle the problem. Recently, lattice and path integral Monte Carlo algorithms were successfully used to study the properties of trapped bosons [19], trapped fermions with attraction [20], fermions in a box with periodic boundary [21]. The properties of a balanced two-component Fermi gas in a harmonic trap were studied by means of the coupled-cluster method [22].
In the following we study the energetic and structural properties of a trapped multicomponent system of one dimensional fermions.

The model Hamiltonian and parameters
We consider a multicomponent one-dimensional Fermi gas at T = 0, trapped in a harmonic confinement of frequency ω. Inspired by the LENS experiment [14], we consider N c spin components of the same atomic species of mass m, with N p particles in each component, the total number of particles being equal to N = N c N p . The system Hamiltonian is given by where g is the coupling constant. In the LENS experiment g was not changed. However, in a more general case, its value can be fine-tuned by changing the magnetic field and exploiting the Feshbach and Olshanii [23] resonances.
In the LENS experiment [14] the number of particles in each spin component was around N p = 20. Such a number is rather small, so it is questionable if all quantities can be precisely described within the local density approximation (LDA). At the same time, this number is already larger than the system sizes which can be accessed with the direct diagonalization methods, as there the complexity grows exponentially with the system size. On the contrary, quantum Monte Carlo methods work very efficiently with the system sizes of interest.
The harmonic confinement defines a characteristic length scale, a osc = h/(mω) which we adopt as a unit of length. We use the inter level spacing of a free confinement, hω, as a unit of energy. Another length scale is associated with the coupling constant, g = −2h 2 /(ma s ), namely the s-wave scattering length a s . We remind that in onedimension, the s-wave scattering length has a different sign compared to the usual three-dimensional case, that is a s is negative for a repulsive interaction, g > 0.
Another peculiarity of a one-dimensional world is that the s-wave scattering length is inversely proportional to the coupling constant, for example, a s = 0 corresponds to an infinite value of g. The third characteristic length is the size of the system. Within the local density approximation it is the Thomas-Fermi size, R T F . While in homogenous system the system properties are governed by a single dimensionless parameter, namely the gas parameter na s , the presence of an external confinement requires, in general, an additional parameter. Within the local density approximation, which cannot describe the Friedel oscillations and is expected to be applicable for large system sizes, R RF a osc , the properties depend on the LDA parameter [24,25] At N a 2 s /a 2 osc → 0 the interaction is infinitely strong and the ground-state energy of N c -system becomes equal to the one of ideal one-component fermions. In the opposite limit of a 2 s /a 2 osc → ∞ the interaction vanishes and the system behaves as N c -component noninteracting Fermi gas.

Methods
We resort to the fixed-node diffusion Monte Carlo (FN-DMC) technique to find the ground-state properties of the system. The proper fermionic symmetry of the wave function are imposed by using an antisymmetric trial wave function. For a given nodal surface the FN-DMC provides the rigorous upper bound to the ground state energy. If the nodal surface of the trial wave function is exact, the FN-DMC obtains the statistically exact ground-state properties of the system. Importantly, the nodal surface in one dimension is known exactly as the fermionic wave function must vanish when any two fermions approach each other.
We chose the trial wave function as a product of determinants S of a singlecomponent Fermi gas and correlation terms which ensure Bethe-Peierls boundary condition [26] between different species α and β, ∂Ψ/∂( Here R is the multidimensional vector which contains the coordinates x α i of all particles and R α contains all particles of component α. The Slater determinant in a harmonic trap is constructed from Hermite polynomials and has a special structure of the van der Monde determinant which can be explicitly evaluated [27] into a N p × N p pair product: with γ = 1/2 value corresponding to the ground-state wave function of a singlecomponent Fermi gas. Slater determinants (4) are antisymmetric with respect to exchange of same-component fermions and define the nodal surface of the total wave function (3). We consider γ as a variational parameter and optimize it value by minimizing the variational energy. The nodal surface of the Ψ T (R) is not affected by γ, but a proper choice of the variational parameter significantly reduces the statistical noise in the ground-state energy and in the correlation functions. For comparison, we also consider a system of N trapped bosons interacting via a contact potential. In this case the trial wave function is chosen the following form:

System energy and contact
It is of a large interest to understand how the interplay between interactions and the quantum statistics affects the energetic properties of the system. The release energy can be measured by suddenly removing the trapping potential and measuring the kinetic energy of the spreading gas. Furthermore, the short-range interactions result in a very intrinsic relation between the interaction energy and the correlation functions, a phenomenon which is truly unique for dilute ultracold gases and is absent in condensed matter. Namely, it was shown in 2008 by Shina Tan [28,29,30] that in three-dimensional Fermi gases many short-range properties of correlation function are related to a universal number, commonly referred as Tan's contact C, which, in turn, is related to the equation of state as In one dimension a similar relation was anticipated in 2003 by Maxim Olshanii and Vanja Dunko [31] (on 1D, see also Ref. [32]). Generally, the contact provides a number of universal relations which connect the short-range correlations to the thermodynamics of the whole many-body system. The universality of the contact was experimentally verified in Bose [33] and Fermi [34,35] gases. In order to obtain the FN-DMC value of contact C we calculate derivative (6) of the ground-state energy E 0 using finite difference. The resulting dependence of the contact on parameter N a 2 s /a 2 osc is reported in Fig. 1. It is convenient to begin the analysis with the simplest case of a two-component system, N c = 2. For a single particle in each component, N p = 1 (green solid squares), the many-body problem reduces to an exactly solvable two-body problem [36] (green solid line). A very good agreement between quantum Monte Carlo data and the analytic result verifies the used numerical approach. For a single particle in each component, the quantum statistics is not yet important. Instead, for N p = 5 (purple down triangles) and N p = 20 (yellow up triangles), particles in each component obey the Fermi-Dirac statistics. Our numerical finding is that the rescaled contact C/N 5/2 has a negligible dependence on the number of particles in the component N p for the considered values of N a 2 s /a 2 osc . From a practical point of view this allows to use the analytic two-particle result as a good approximation to the contact in a two-component system.
For fermions with three components, N c = 3 (pink left triangles), and six components, N c = 6 (red right triangles), the absolute value of the contact is increased. At the same time the qualitative dependence on N a 2 s /a 2 osc remains the same. Inspired by a good agreement between the contact in a two component fermions and two bosons, we also perform the simulations for the case of one-component bosons. In this case Ψ T (R) = Ψ bos J 1 (R)Ψ G (R), where in the first term a b = a. The DMC results for 3 and 6 bosons are shown in Fig. 1 (black open circles and squares, correspondingly). One can see that also in the case of a multi-component fermions the contact is similar to that of a gas of N c bosons in the considered range of parameters.  The interaction energy, E int = g Nc α<β Np i,j=1 δ(x α i − x β j ) , can be calculated by using Hellmann-Feynman theorem as E int = −a s dE 0 /da s . In other words, there is a close relation between the contact (6) and the interaction energy for a contact potential,

bosons 6 bosons
We show the interaction energy in Fig. 2. The interaction energy vanishes in the limit of non-interacting species, g → 0, which corresponds to N a s /a ho → ∞. In the opposite limit of very strong interspecies interactions, g → +∞, there is a node in the wave function when particles approach each other (refer to Eq. (3) with a s = 0). As a result the interaction energy vanishes also in the limit N a s /a ho → 0. The maximum in the interaction energy is observed for N a s /a ho ≈ 1.
Similarly to the contact, we find that the most important effect is the number of components N c and not the number of particles in each of them, at least for the considered parameter range.

Correlation functions
In one-dimensional geometry there exists an unusual relation between the interactions and statistics, which is absent in higher dimensions.
By tuning Busch et al. the interspecies interaction strength it is possible to change the correlations from uncorrelated "ideal boson-like" case to the "ideal fermion-like" case.
According to Girardeau's mapping [37], infinite repulsion between one-dimensional particles mimics the Pauli principle and the absolute value of the wave function is the same. This holds true both for a system of one-component bosons (Tonks-Girardeau gas) but also a system of two-component Fermions. For the case N c = N p = 1 the fermionization of two distinguishable fermions was experimentally demonstrated in Selim Jochim's group [12]. From the Girardeau's mapping it immediately follows that the diagonal properties, which depend on the square of the absolute value of the wave function, are the same in both system. Instead, off-diagonal properties are different. From that it is highly instructive to analyze the evolution of both one-and two-body correlation functions. We consider the density profile n(x) and the density-density correlation function g(x), which obey the Girardeau's mapping, and experimentally relevant momentum distribution n(k), for which the mapping does not apply.

Density profile
We calculate the density profile n(x) in a two-component Fermi gas, N c = 2, with N p = 5 particles in each component, as a function of the interaction parameter N a 2 s /a 2 osc and show it in Fig. 3). For no interaction between two species, N a 2 s /a 2 osc → ∞, the density profile is given by that of a spin-polarized Fermi gas with N p particles (dashed line in Fig. 3). In this case, the interaction energy E int = 0 (refer to Fig. 2). As the interaction strength is increased, the interaction energy increases and the system size becomes larger. Stronger repulsion pushes the particles to the edges of the trap, making the density in the center drop down. For infinite interaction strength, N a 2 s /a 2 osc → 0, Girardeau's arguments predict that wave function can be mapped to that of an ideal single-component Fermi gas of N c N p particles. Indeed, we see that in this limit the density profile is that of an ideal Fermi gas with 2N p particles (solid line). In between we see a smooth crossover from one regime to the other. In other words we observe a fermionization of two fermion species, similarly to experimentally observed N p = 1 case [12] (see also Ref. [38] for Bose-Fermi mixtures).  It is important to note that the Friedel oscillations, clearly seen in Fig. 3 are completely missed within local density approximation and cannot be predicted using Bethe solutions for homogeneous gas. Furthermore, here we observe an interesting parity effect. When interspecies interaction is absent, there is a maximum at n(x = 0). For infinite interspecies interaction, there is a minimum at n(x = 0). In addition we see a doubling in the number of peaks which sometimes is interpreted as a transition from 2k F -Friedel to 4k F -Wigner oscillations [39].
On the other hand, even if LDA misses the shell structure, it is still capable to predict the overall shape and the cloud size. As well a polytropic approximation to the equation of state permits to understand some basic features analytically. In the non-interacting limit, N a 2 s /a 2 ho → ∞, the density profile is a semicircle, n(x) ∝ (µ − mω 2 x 2 /2) with µ = N ph ω. In the strongly-interacting limit, N a 2 s /a 2 ho → 0, the density profile is again a semicircle but of a larger size, µ = 2N ph ω. In the limit of large number of components, N = N c and N p = 1, the density profile is a semicircle for infinite interaction, µ = Nhω, is an inverted parabola for weak interactions, n(x) ∝ (µ − mω 2 x 2 /2); and is a free-harmonic oscillator Gaussian profile for zero interactions (not captured by LDA).

Pair correlation function
The pair correlation function g(x) is proportional to the probability of finding a pair of particles separated by a distance x. It is a diagonal quantity and is a subject to Girardeau's mapping, so we expect to see fermionization for strong repulsive interactions. In a trapped system the long-range decay zero comes from vanishing particle density at the edges of the trap and is eventually a one-body effect. From the quantum-correlation perspective the most interesting behaviour is that of short and intermediate distances.
For the intra-component correlation function the Pauli exclusion principle implies that g ↑↑ (x = 0) = 0. Instead, for inter-component correlation function g ↑↓ (x) the value at the contact, x = 0, is related to C using the following relations: In Fig. 4 we show the inter-component pair-correlation function g ↑↓ (x) in a twocomponent system. One can observe how the decrease of the interaction parameter N a 2 s /a 2 osc (which means the increase of the interaction strength) the value g ↑↓ (x = 0) decreases until it reaches the zero value corresponding to a fully fermionized system. The horizontal arrows in Fig. 4 shows g(0) at given N a 2 s /a 2 osc calculated from the FN-DMC energy through the the contact, see Eq. 6. We observe a good agreement between these two methods of the calculation of the contact.

Momentum distribution
An important experimentally observable quantity is the momentum distribution n(k). Being an off-diagonal quantity, it will not fermionize in the limit of infinitelystrong repulsion and will be essentially different from n(k) of an ideal Fermi gas. Also, the use of the Bethe ansatz methods is very limited in calculation of n(k).
We show the momentum distribution of a two component system in Fig. 5. One observes that with the decrease of the parameter N a 2 s /a 2 ho (increasing the interactions) the amplitude of n(k) increases and the distribution becomes narrower. This latter effect is opposite to the particle distribution becoming wider in the real space, see To see that we plot the same results on a double logarithmic scale in Fig. 5(b). The black lines at large k show the function C/k 4 , where C is the value of the contact obtained from the numerical derivative of the energy. We find that asymptotic regime is reached for momenta larger than ka osc ≈ 5.

Conclusions
To conclude, we studied the ground-state properties of a one-dimensional multicomponent Fermi gas in presence of a harmonic confinement (N c components with N p particles each). Interspecies interaction is considered to be of a contact form, We used fixed-node diffusion Monte Carlo (FN-DMC) method to calculate the energy and correlation functions. A unusual feature of one-dimensional geometry is that there is an interplay between quantum statistics and contact-interaction. According to Girardeau's mapping [37], the wave function in the limit of infinite repulsion, g → +∞, can be mapped to the wave function of an ideal Fermi gas. In this regime, the multicomponent gas fermionizes having the same energy and diagonal properties as a single-component ideal Fermi gas of N = N c N p particles. The fermionization is observed in the energy, the density profile, the pair correlation function but not in the momentum distribution.
We calculate the Tan's contact C for different number of particles N p in each of N c components. We find that for a given number of components N c , the C/N 5/2 has an almost universal dependence on parameter N a 2 s /a 2 osc for the considered interaction and number of particles. Importantly, we discover that in our case the contact can be well approximated by considering only one particle in each component N p = 1. For a two-component Fermi gas N c = 2 in harmonic confinement this allows to use analytic result of Busch et al. [36]. Also we check numerically that for N c = 3 and N c = 6 components the value of the Tan's contact is close to that of a Bose system with 3 and 6 particles. We verify that the Tan's contact, calculated from the derivative of the groundstate energy (6), as well provides the description of the interaction energy (7), value of the pair-correlation function at zero (8), and large-momentum tail of the momentum distribution (9).

Acknowledgements
NM would like to acknowledge the Nanosciences Foundation of Grenoble for financial support. GEA acknowledges partial financial support from the MICINN (Spain) Grant No. FIS2014-56257-C2-1-P. The Barcelona Supercomputing Center (The Spanish National Supercomputing Center -Centro Nacional de Supercomputación) is acknowledged for the provided computational facilities.
Appendix: contact for weak interaction, LDA The limit of weak interaction can be analyzed assuming ideal Fermi gas density profile for each component and taking into account the interaction potential perturbatively. For energetic calculations it is sufficient to assume the LDA shape, which for a single component provides the correct energy even for N p = 1 The interaction energy is (11) Considering an unperturbed LDA density profile (10) we get from (11) the interaction energy and the contact We note that Eqs. (12)(13) are derived assuming weak interaction, a s → −∞ and a large number of particles N p in each component. A priori, it is not obvious if N p = 1, E int /(hω) = −0.77a osc /a s , expression (12) can be used at all. Nevertheless, a comparison with the exact expression (16) is less than 5% of the exact coefficient in the asymptotic expansion. The major difference between multiple-components and two-particle cases is the presence of a combinatoric N C (N C − 2)/2 term in Eq. (12).
Appendix: contact for weak interaction, two particles The problem of two particles in a harmonic trap can be solved analytically [36]. The harmonic confinement permits to separate the problem into center of mass (CM) and relative motion (rel). The ground-state energy of the CM motion is E CM =hω/2. The energy of the relative motion E rel is defined as a solution to the following equation (harmonic oscillator units are used) For the same value of the s-wave scattering length a s Eq. (14) permits multiple solutions for the energy, which correspond to multiple level structure. For a weakly interacting case, a s → −∞, the energy of the relative motion is close tohω/2 value and Gamma functions in Eq. (14) can be expanded into series around this point. The resulting total energy E = E rel + E CM is