Quantum Coherent States of Interacting Bose-Fermi Mixtures in One Dimension

We study two-component atomic gas mixtures in one dimension involving both bosons and fermions. When the inter-species interaction is attractive, we report a rich variety of coherent ground-state phases that vary with the intrinsic and relative strength of the interactions. We avoid any artifacts of lattice discretization by developing a novel implementation of a continuous matrix product state ansatz for mixtures and priorly demonstrate the validity of our approach on the integrable point that exists for mixtures with equal masses and interactions (Lai-Yang model) where we find that the ansatz correctly and systematically converges towards the exact results.

Since the experimental achievement of quantum degeneracy in weakly interacting atomic gases, first for bosons [1] and soon after for fermions and boson-fermion mixtures [2,3], there has been a tremendous interest in the study of such systems. They provide an ideally clean and controllable experimental platform for fundamental studies of unitary quantum gases and for their use in applications ranging from quantum simulation to metrology. For the case of mixtures, the new experimental platform reignited a long-standing theoretical interest connecting back to 4 He- 3 He solutions [4] and related to superfluidity and phase separation; but now combined with new practical considerations, since mixtures turned out to be also a useful stepping stone for the evaporative cooling of fermions with otherwise vanishing s-wave cross sections (a method known as sympathetic cooling [3,4]).
While the initial theoretical studies of gas mixtures (with mixed statistics) considered three-dimensional (3D) gas clouds [5], soon the interest expanded to include one-dimensional confinement (via elongated traps or optical lattices) [6][7][8]. Those early studies unveiled a rich Luttinger-liquid phenomenology with instabilities towards collapse, demixing and pairing depending on the sign and strength of the interactions [8]. Subsequent studies paralleled the even greater theoretical activity attracted by two-component Fermi mixtures [9], that accompanied the experiments as they also evolved from 3Dtrapped clouds [10] to the realm of 1D confinement [11]. A similar evolution in the experimental study of Bose-Fermi mixtures has the potential to be equally rich and interesting.
The physical modeling of (one-dimensional) Bose-Fermi mixtures of atomic gases with short-range (contact) interactions is well approximated by the generalized , where the fieldsψ † α (x) create bosons or fermions with mass m α (for the fermionicity α ∈ {b = 0, f =1}, respectively) and obey the appropriate (anti)-commutation They interact with a contactpotential strength g αβ = 2 ω ⊥ a αβ , where a αβ is the corresponding scattering length between the two species and ω ⊥ is the common transverse trapping frequency [13]. Since fermions avoid each other, for contact interactions (with a range smaller than the exchange-correlation hole) the strength becomes arbitrary and we can set g ff = 0. This model has been studied (away from integrable points) using mean-field and field-theoretical approaches (such as bosonization, in the hard-core-boson limit) [6,8,14]. However, most of the further studies of onedimensional Bose-Fermi mixtures [7,15] relied on lattice discretizations and mappings to the Bose-Fermi Hubbard model; plus the use of wave-function approximations such as the BCS mean-field, Gutzwiller and Jastrow Ansätze [16], or numerical methods such as the matrix-productbased density matrix renormalization group (DMRG) [17].
During the last decade, a generalized coherent-state version of the matrix-product variational Ansatz for quantum fields in a one-dimensional continuum was successfully put forward [18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][36][37]. This collective work explored the physics of a large number of systems of single or multi-component gases with a given type of statistics. We present here an extension of those developments that demonstrates an efficient implementation of the continuum matrix-product-state (cMPS) Ansatz for a twocomponent gas with mixed statistics. We apply it to the study of the Lai-Yang model, first validating the method against the Bethe Ansatz at the (non-trivial) integrable point and then studying the nature of pairing tendencies in certain regions of the non-integrable regime.
A cMPS for a binary (or multicomponent) gas mixture, for a system of length L and periodic boundary conditions, has the general form [21] [38] (4) where the complex-valued local matrices Q(x) and R α (x) arXiv:2110.13899v2 [cond-mat.quant-gas] 24 May 2022 act on an auxiliary space of dimension D, called the bond dimension. In addition, α is the index for the atom type (b or f),Î is the identity operator on the Fock space, Tr aux is a trace over the auxiliary space, Pe is a path-ordered exponential, and |0 is the bare vacuum state annihilated by all the field operators. Henceforth, we will consider flat-bottom trapping potentials and adopt a phase-modulated uniform Ansatz [27] given by Q(x) = Q and R α (x) = e iqαx R α , where Q and R α are three position-independent matrices, and q α are two additional real-valued variational parameters (continuum in the large-L limit).
For regularity, the R α have to obey the same local algebra as the matching field operators, [R α , R β ] (αβ) = 0. In particular, nilpotency in the fermion sector, R 2 f = 0, has to be implemented exactly for improved numerical stability [27]. Thus, ignoring zero eigenvalues, R f has to have an unrescaled Jordan form [39] composed entirely of 2×2 blocks with zero Jordan eigenvalue (i.e., proportional to the spin-rising Pauli matrix, σ + ); notice this restricts D to take only even values. The inclusion of scaling parameters in the Jordan blocks provides better convergence for mixtures with (very) different densities of each species. Explicitly, one can write R f = σ + ⊗ Γ and enforce commutation with its boson counterpart by we do by solving for D b and keeping the other matrices arbitrary for variational optimization [40] [41].
The cMPS norm is given by χ|χ = Tr e T L where That guarantees a cMPS has unit norm in the thermodynamic limit [18,22]. We can then have a right-identitynormalization of T by taking A is an arbitrary anti-Hermitian matrix [42]. We find that a numerically accurate normalization of T is important for the convergence properties of the algorithm.
The α-atom number density can be readily computed as and similarly for the corresponding kinetic-energy density by using Interaction terms and correlators admit similar (longer) expressions in which the R α matrices replace the field operators [18,22,27].
In order to demonstrate the accuracy of the cMPS Ansatz, we focus first on the Lai-Yang integrable point [12,[43][44][45], which corresponds to equal masses (m b = m f = m) and interactions (g bb = g bf = g fb = g > 0 [46]) between bosons and fermions. In experiments, the

E/E [%]
Lai-Yang (nb = nf) Lieb-Liniger Variational cMPS ground-state energy, e(γ) = 2mE/(n 3 2 L), for the Lai-Yang model with different bond dimensions as a function of the fermion-density fraction and compared with the Bethe Ansatz result. Inset: cMPS convergence with bond dimension of the relative energy error, ∆E = Ecmps − E, as a percentage of the Bethe-ansatz value, for n f /n = 0.5 and 0 (the latter is the Lieb-Liniger limit that we computed with a bosons-only cMPS that has more variational parameters for a given D and yields lower energies).
first condition can be approximately fulfilled by considering different isotopes of the same atom, and the second one by tuning the interactions via shape and Feshbach resonances (deviations break integrability, but can still be modeled within cMPS [32]). There have been several studies of this special case using the Bethe Ansatz and framed in the modern context of cold-atom gases [47][48][49][50]. They found that, despite earlier suggestions, the demixing tendency is absent for the homogeneous system (but one can still have trap-induced Bose-Fermi separation; cf. Ref. 51). We will not explore those questions here; rather, we will focus on comparing variational groundstate energies with the exact results and assessing the convergence properties of cMPS; see Fig. 1. We will work with units in which = 1 and m = 1 [52], will take n = n b + n f = 1/4 per unit length (in some arbitrary units), and the interaction strength will be given by setting the dimensionless parameter γ = g/n = 8 [53].
We observed good systematic convergence of the cMPS energies as a function of bond dimension-as expected, since the Ansatz captures progressively more and more multi-particle entanglement. For two-component Bose-Fermi mixtures, convergence in D is slower than in the Lieb-Liniger case, but faster than for the Gaudin-Yang model (the two-fermion Ansatz can be constructed as an outer nesting layer on the one given here and is further restricted to D being a multiple of 4;cf. Ref. 27).
We minimized the energy of the system directly in the L → ∞ limit and tested several standard localoptimization algorithms. We found that a principal-axis minimization was usually the best strategy [54], although other ones, such as the Nelder-Mead simplex method, also converged well. The global optimization can be done using simulated annealing, but we found that repeated random starts and parametric variations provided a simpler strategy and gave comparable results for not-toolarge bond dimensions. We targeted fixed particle densities using (augmented) Lagrange multipliers and penalties. The single-optimization running times on regular hardware ranged from a few seconds for D = 2 to several hours for D = 8.
Having established the viability of the cMPS Ansatz for the study of mixtures, we next move on to address the case of regimes with attractive interactions between bosons and fermions, g bf = g fb = −g < 0. All the while, the interactions between bosons shall remain repulsive, g bb /g = G > 0, in order to guarantee the thermodynamic stability of the system by preventing bosonic collapse. The large-G limit corresponds to the Tonks-Girardeau regime [55] in which the effective exclusion statistics of the bosons gradually approximates that of fermions and the low-energy spectrum of the model converges to that of a Gaudin-Yang gas. One would thus expect a tendency to (algebraic) superconducting order as in the two-fermion case; however, while some of this intuition is correct, the exchange statistics of the bosons remains bosonic (they tend to hard-core bosons) and the emerging physical picture is considerably more subtle, as shall be established below.
Our implementation of the cMPS Ansatz constitutes a sui generis matrix representation of generalized coherent states, and as such it is particularly well suited to capture the emergence of ground states with spontaneously broken U (1) symmetries and off-diagonal (quasi) long-range order. For instance, bosonic field operators can acquire a non-zero vacuum expectation value (vev), ψ b (x) , that signals the occurrence of (quasi) Bose-Einstein condensation (BEC), a state with macroscopic quantum phase coherence [56]. We shall refer to the absolute values of such appropriately normalized vevs as coherence order parameters (cf. Ref. 57). They take values in the unit interval and a comprehensive set of them is displayed in Fig. 2 for a wide range of values of γ = g/n and G.
Due to their statistics, fermions cannot condense as single atoms, but they can pair up and condense as a molecular BEC (localized pairs), or a more BCS-like state (extended pairs), or any intermediate scenario (BEC-BCS crossover). Such pair-coherent states (also possible for bosons) are signaled by the vevs ψ α (x)ψ α (x + δx) , where x is arbitrary due to translation invariance and δx is chosen to maximize the amplitude (and corresponds to the most likely atom-atom distance in the pair, which has to be nonzero for spinless fermions due to exclusion but turns out to be zero for bosons). Naïvely, one would want to consider also mixed-species pairs, but their combined Fermi-Dirac statistics precludes condensation

FIG. 2. Normalized coherence order parameters for a Lai-
Yang gas with repulsive inter-boson and attractive bosonfermion interactions. The calculations were done using D = 4 cMPS states and the typical error is expected to be in the 5-10% range. Contrasting behaviors are clearly seen between the strongly (γ 1) and weakly (γ 1) interacting regimes.
(even in the Tonks-Girardeau regime). Rather, Bose-Fermi coherence manifests in composite order parameters (cf. Refs. 15,58). After exploring different correlators, we are led to define the local bf-molecule field ψ bf (x) ≡ψ b (x)ψ f (x) and consider its pair-of-pairs (PoP) vev defined as above but with α → bf.
Schematic representation of the different quasilong-range-order ground-state phases present when the interactions among bosons are repulsive while bosons and fermions attract each other, as a function of interaction strengths. The dominant orders are indicated in each case.
The combined consideration of all the order parameters defined above, yields a comprehensive picture of a ground-state phase diagram with four starkly different sectors and interpolating crossover regions between them. This information is schematically summarized in Fig. 3 to guide the discussion. Rather than comparing algebraic exponents as in bosonization-based studies, cMPS forces the breaking of U (1) symmetries and allows the direct comparison of coexisting coherence-order-parameter amplitudes (normalized so that the dominant phases can be identified by their larger relative magnitudes in Fig. 2). The simplest case is the weakly interacting limit deep into the third quadrant of Fig. 3, in which ψ b is the dominant order (approximately saturating the bound) signaling boson condensation. Notice in the top panel of Fig. 2 that ψ bψb is also comparably large, as expected for true higher-order coherence [57], while the other order parameters involving fermions are highly suppressed in comparison. Increasing G and moving into the fourth quadrant, the amplitude of the ψ b order parameter is partially suppressed and the coherence reduces to first-order only. On the other hand, the ψ fψf order is enhanced in the process, with extended BCS-like pairing, and we find a scenario of combined coherence (compatible with the polaronic picture of earlier studies [58]).
If, instead, we keep G small and increase γ, going into the second quadrant, the single-boson coherence is suppressed to very low values while the two-boson one remains close to maximal, and in addition the two-fermion coherence is greatly enhanced and more BEC-like than in the fourth quadrant (see the middle panel of Fig. 2). This is interpreted as Bose-Fermi mutualism, with the bosons providing the glue for the fermionic pairing while, simultaneously, the fermions do the same for the bosonic pairing. The two order parameters are not mixed (they involve different degrees of freedom), but neither condensate would be able to exist without the other. Finally, when G and γ are both large, we go into the first quadrant of Fig. 3. We find that all of the purely bosonic or purely fermionic orders are suppressed to zero (see the top and middle panels of Fig. 2), while only the mixed ψ bfψbf four-particle coherence remains present in the system (bottom panel of the same figure). The latter corresponds to tightly bound Bose-Fermi molecules that condense in loosely bound pairs. Moreover, this PoP order is also enhanced in the intermediate-coupling region at the center of the phase diagram, where it coexists with a strength similar to that of the dominant orders from the other quadrants in a large mixed-coherence crossover region.
We have focused on the case of balanced populations of equal-mass Bose and Fermi atoms, but those conditions might well not be the most easily achievable in experiments. We also considered the case of a 20% density imbalance (both ways) and found that the phase diagram remains the same asymptotically while the location of the crossover boundaries shifts with the density ratio. One would expect similar results as a function of mass imbalance and cMPS would be an ideal tool for that study (based on our past experience with Gaudin-Yang systems [32]); however, since the space of parameters is large, we leave a focused study for the future, to be guided by the parametric choices dictated by experimental considerations. Similarly to the case for twofermion gases, a convenient setup might be a quasi onedimensional array of tubes created by a transverse optical lattice. This would turn the ordering tendencies of the system into true long-range order by exploiting the 1D-3D crossover and can have a number of additional experimental benefits [11,15,59]. On the other hand, the inclusion of a (weak) longitudinal optical lattice would break translation invariance and open up the possibility of density-wave orders that do not exist in the continuum (as found in studies based on the Bose-Fermi Hubbard model). The cMPS-ansatz implementation can be extended to that case as well, as has been done already for Lieb-Liniger gases [33]. This carries the advantage of treating translation invariance, or lack thereof, without the additional umklapp scattering introduced by the lattice discretization needed for the use of standard MPS methods (which, moreover, require the truncation of the local bosonic Hilbert spaces at each lattice site and introduce biases in the capture of BEC order). Finally, our implementation of the particle statistics turns out to be quite elegant and simpler than in the MPS/DMRG set-ting, (cf. Ref. 17 and the large body of subsequent work). This work was funded by the NSF (Grant No. PHY-1708049) and at times used computing resources from the Ohio Supercomputer Center (OSC) [60]. In its initial stages, the work benefited from an extended visit to ITAMP (Harvard-Smithsonian Center for Astrophysics). We also acknowledge the hospitality of the IIP-UFRN (Brazil).