Self-binding of one-dimensional fermionic mixtures with zero-range interspecies attraction

For sufficiently large mass ratios the attractive exchange force caused by a single light atom interacting with a few heavy identical fermions can overcome their Fermi degeneracy pressure and bind them into an $N+1$ cluster. Here, by using a mean-field approach valid for large $N$, we find that $N+1$ clusters can attract each other and form a self-bound charge density wave, the properties of which we fully characterize. Our work shows that there are no fundamental obstacles for having self-bound states in fermionic mixtures with zero-range interactions.


Introduction
According to our current understanding of Nature, big composite self-bound objects (nuclei, atoms, molecules, liquids, solids, etc.) are formed due to attractive finite-range forces originating from exchanges of bosons (gluons, photons, etc.) However, one of the archetypal fermionic models, the two-component mass-balanced Fermi gas with zero-range attraction [1], exhibits no self-binding. Increasing the attraction in such a gas leads to the formation of dimers consisting of fermions of different sort. The dimers repel each other in any space dimension because of the Pauli exclusion principle for the constituent fermions [2][3][4][5][6][7][8][9]. The Pauli "repulsion" is also believed to be the main mechanism preventing binding of four neutrons, although the topic remains controversial since internucleon interactions are not zero range [10].
Introducing mass imbalance into the model can lead to binding of mesoscopic clusters of type N + 1, where the exchange force of a single light atom overcomes the degeneracy pressure of a small Fermi sea of N heavy atoms. Such clusters with N up to 5 have been studied by exact few-body techniques in all dimensions [11][12][13][14][15][16][17][18] and we have recently developed their mean-field theory in one dimension valid for large N [19].
Can a two-component Fermi mixture with zero-range interactions become self-bound in the thermodynamic limit? A good starting point to answer this question is to understand whether two clusters of the type N + 1 can stick together. This problem is nontrivial; the light atoms should be sufficiently light to ensure attraction for the heavies, but, on the other hand, their own degeneracy pressure (inversely proportional to the light mass) can hinder binding. No evidence of such binding has been reported. Here we can cite rather extensive studies of the fermionic 2+2 system and the dimer-dimer scattering problem in the mass-imbalanced case [20][21][22][23][24][25][26]. Although not fully comprehensive (i.e., not all dimensions and possible mass ratios covered), these studies are consistent with the scenario that the 2+2 fermionic system is either unbound or breaks into two repulsive dimers or into a heavy-heavy-light trimer plus a free light atom when the trimer gets below the two-dimer threshold. Naidon and co-workers [27,28] estimated that three-dimensional 2+1 trimers repel each other.
In this article we show that in one dimension N + 1 clusters can arrange themselves into a self-bound configuration, at least, for sufficiently large N . To this end we use the mean-field theory based on the Thomas-Fermi approximation for the heavy atoms, valid in the limit N ≫ 1 [19]. In this case, the system behavior is governed by a single parameter α = (π 2 /3)N 3 m/M . We find that two clusters bind for 2.3 < α < 9.4. Interestingly, instead of merging, they stay at a finite distance from each other and keep a doublepeak density profile (see the gray dotted curve in Fig. 2). Below α = 2.3 this state becomes metastable and the true ground state corresponds to two N +1 clusters at infinite separation. Our calculations show that three or more clusters can form a self-bound charge density wave with one light atom per period. We describe bulk properties of these states by a fully analytic weak-modulation theory, the Peierls instability [29] emerging as a complementary explanation of the modulation.

Model and N + 1 cluster solution
We address the problem of N h fermions of mass M interacting with N l fermions of mass m through the mean-field density functional where g < 0 is the heavy-light interaction constant, ϕ i (x) are the wave functions of the light atoms, and n(x) is the density profile of the heavy atoms. The term ∝ n 3 (x) is the kinetic energy of an ideal Fermi gas taken in the Thomas-Fermi local-density approximation [19,30]. The model (1) thus requires weak interactions (a ≫ λ h , λ l ) and n(x) slowly varying on the scale ∼ λ h . Here, a = −(m + M )/(mM g) is the scattering length and λ h ∼ 1/n and λ l ∼ |ϕ i /∂ x ϕ i | are, respectively, the typical de Broglie wave lengths of the heavy and light atoms. The first line in Eq. (1) defines the total energy E, which we seek to minimize subject to constraints ϕ * i (x)ϕ j (x)dx = δ ij [31] and n(x)dx = N h . The normalization constraints are taken into account by introducing Lagrange multipliers ϵ i and µ.
One can show that up to an overall scaling factor, the behavior of the system satisfying Eq. (1) is governed by two dimensionless parameters. We choose the first to be N l and the second to be α = (π 2 /3)N 3 m/M , where N = N h /N l . Indeed, introducing the characteristic size λ = 1/(2m|g|N ), new coordinate u = x/λ, and rescaling whereμ = 2mN µλ 2 < 0,ε i = 2mϵ i λ 2 and the normalization constraints are now φ * i (u)φ j (u)du = δ ij and ñ(u)du = N l . We minimize Ω imposing that the variational derivatives of Eq. (2) with respect toφ * i andñ vanish. These conditions, respectively, lead to the equations The energy of the system then equals Let us briefly summarize the main results obtained for the case N l = 1 [19]. The N + 1 cluster exists for α < 12 and we show its energy, denoted by E N l =1 N +1 , as a function of α in the inset of Fig. 1. In the limit α → 0 the heavy atoms are much more localized than the light one, and the system can be described as a light atom bound by a point- With increasing α the heavy atoms get more freedom and their chemical potential grows till it reaches µ = 0 at α = 12. This corresponds to the right endpoints of the red dashed curves in Fig. 1 Beyond this point the droplet cannot accomodate more heavy atoms. The validity of the model (2) is verified by the following arguments. Note that Eqs. (3) and (4) are dimensionless and, for α ∼ 1 and N l ∼ 1, the spatial extent of the cluster is of order λ (in original units). The condition of the slowly varying density can be written as |∂ x n(x)| ≪ n 2 and translates to N ≫ 1. This inequality is equivalent to (m/M ) 1/3 ≪ 1 since we are mainly interested in α ∼ 1. We thus have a ≈ −1/mg, or, in rescaled units,ã = a/λ ≈ −2N , which is much larger thanλ h ∼ 1/N andλ l ∼ 1, ensuring weak interactions.

Binding of two or more N + 1 clusters
In the case N l = 2 we solve Eqs. (3) and (4) iteratively. Namely, we diagonalize Eq. (3) assuming a certain initialñ, substitute the obtainedφ i into Eq. (4), tuneμ to satisfy the normalization constraint forñ, and repeat the procedure till convergence. In this manner, we obtain three types of solutions. The first type is two isolated N + 1 clusters, the second is a bound state of two N + 1 clusters, and the third is a 2N + 1 cluster plus a free light atom. The first type is the ground state of the 2N + 2 system for 0.16 < α < 2.3.
N +1 shown in the inset. For small α, N l isolated N + 1 clusters can lose one of their light atoms and rearrange into N l − 1 isolated clusters with larger N . The final-to-initial energy ratio (with minus sign) is shown for N l = 2 (magenta dotted) and N l = 3 (green dash-dotted). The crossing region (black frame) will be shown in more detail in Fig. 3.
The second type is realized in the region α > 1.6. The corresponding energy per cluster (we denote it by E N l =2 N +1 ) is shown in Fig. 1 as the gray dotted curve. For our purposes it is convenient to normalize all energies in Fig. 1 to the energy of an isolated N + 1 cluster, shown as the horizontal red dashed line. For α > 2.3 the bound-state configuration is the ground state, and in the region 1.6 < α < 2.3 it is only dynamically stable (thermodynamically it prefers to break into isolated clusters). This state is characterized by a density profile with two maxima separated by a finite distance (see Fig. 2). Below α = 1.6 the metastable state disappears and the clusters unbind. To qualitatively understand this phenomenon, we have performed a variational analysis by takingñ(u) as a sum of two Gaussians with variable widthσ, placed at distanceξ from each other. The minimization of Eq. (5) with respect toσ then gives the energy as a function ofξ. The curves E(ξ) obtained in this manner can feature a (meta)stable minimum at finiteξ. They are very similar to what we obtain for infinite N l (see Sec. 4 and Fig. 4). The left panel in Fig. 4 corresponds to the critical α where the metastable minimum disappears. This is the point where the minimum and maximum of E(ξ) merge, creating nonanalytical singularities in both the optimalξ and the energy. The same nonanalytic behavior is observed in our exact (not variational) numerics. In fact, the critical points in Fig. 1 are determined by gradually decreasing α and monitoring the distance between the peaks and the energy of the systems.
We find that the curves E(ξ) obtained by the variational procedure are characterized by a repulsive tail at largeξ. That the clusters repel each other at large separations is due to the exchange of the identical light fermions. The mechanism can be understood in the Born-Oppenheimer approximation and is similar to the long-distance repulsion between two heteronuclear mass-imbalanced dimers [23]. The minimum at finiteξ appears because heavy atoms distribute their density in an optimal manner to provide binding in spite of the degeneracy pressure. This phenomenon, which is the main result of this work, is not obvious and rather subtle (note relatively low binding energies). As we increase α beyond 9.4, similarly to the case N l = 1, the chemical potential µ crosses zero and becomes positive. This simply means that, in free space, the cluster will eject the excess of heavy atoms, effectively decreasing α to subcritical values. Still speaking about solutions of the second type we find that the phenomenon of selfbinding persists for N l > 2. Solving Eqs. (3) and (4) for N l up to 5 we observe that clusters tend to form a regular chain or polymer, with the number of density peaks equal to N l . To avoid cluttering in Figs. 1 and 2 we show only the case N l = 3 (dash-dotted) in addition to N l = 1, 2 already discussed. The three-cluster bound state shows qualitatively the same behavior as the two-cluster state. In brief, it is stable above and metastable below α = 2.4. Below α = 2.2 the system energetically prefers to break into three isolated N +1 clusters. An interesting feature of this system [32] is that in the region 2.2 < α < 2.4 (more precisely 2.21 < α < 2.38) its true ground state is an N l = 2 chain with N h = 2N ′ plus a single (3N − 2N ′ ) + 1 cluster, where N ′ is determined by the minimization of . We observe similar behavior for higher N l . The curves corresponding to the binding energies (per cluster) of the chains with different N l bundle together in the region α ≈ 2.4 ± 0.1. In Fig. 3 we plot the results for N l =1, 2, 3, 4, 5, and ∞. The determination of the ground state for a given N l in this region is complicated since a longer chain can break into shorter chains with generally different α. We note, however, that in this region all partitions of a chain are almost degenerate and correspond to (meta)stable states separated by energy barriers similar to the one shown in the middle panel of Fig. 4. For larger α, sufficiently far from the crossing region, clusters do prefer to merge into a single chain since longer chains feature higher binding energy per cluster.
The third type of solutions of Eqs. (3) and (4) realizes for small α when N l isolated N +1 clusters eject one light atom forming N l −1 isolated N ′ +1 clusters with N ′ = N N l /(N l −1), the new configuration becoming energetically favorable for (N l − 1)E N l =1 N ′ +1 < N l E N l =1 N +1 . In Fig. 1 we show the quantity (N l − 1)E N l =1 N ′ +1 /N l |E N l =1 N +1 | for N l = 2 (magenta dotted) and N l = 3 (green dash-dotted). These curves cross the line -1 at α = 0. 16 Figure 3: Zoom into the framed region in Fig. 1. Here, in addition, we show the energies for N l = 4 (magenta dotted) and N l = 5 (green dash-dotted). The curves show that for a given N l the longest possible chain or N l isolated clusters are not the only possible ground states. For instance, for N l = 3, the energetically optimal configuration can be a cluster with one light atom isolated from a cluster with 2 light atoms (both clusters having generally different values of α, see text).

limit, where it is obtained from the condition ∂[E N l =1
N +1 /N ]/∂N = 0. In principle, starting from a long chain of isolated clusters and trying to follow its ground state by decreasing α below 0.47, the system will lose light atoms one by one till it eventually ends up in the state where all N h heavy atoms are bound by a single remaining light atom. We should note, however, that these transitions are associated with a global redistribution of the heavy particles such that the system will likely get stuck in metastable states with "wrong" N l . This is because these transitions are not associated with max[ϵ i ] crossing zero (isolated clusters individually never lose their light atoms).

Infinite chain analysis
We now go back to the second type of solutions and discuss bulk properties of self-bound chains in more detail. We assume thatñ(u) is periodic with the modulation lengthξ and that light atoms are filling the first Brillouin zone of the lattice (we have one light atom per modulation length). We aim to calculate the energy per cluster, which we denote E N l =∞ N +1 , as a function ofξ. In principle, Eqs. (3)(4)(5) are suitable for the task. In this casẽ ϕ i become Bloch functions and i is the real Bloch wave vector in the first Brillouin zone, i.e., i ∈ (−π/ξ, π/ξ]. It is however convenient to rescale the coordinate again, introducinḡ u = u/ξ. Making related changes and rescalings (we mark new rescaled quantities by a bar), we arrive at the following formulation of the problem. The energy per cluster is given by where the spectrumε p is determined by the equation with the periodic boundary condition χ p (0) = χ p (1). The function χ p is the periodic part of the Bloch wave function corresponding to the wave vector p ∈ (−π, π]. The densityn is given byn whereμ =μξ 2 , and the normalization conditions read 1 0 |χ p (ū)| 2 dū = 1 and 1 0n (ū)dū = 1.
A few examples of these curves are shown in Fig. 4. We find that there is always a local minimum atξ = ∞ corresponding to isolated noninteracting clusters We also see that E N l =∞ N +1 (ξ) decreases withξ for smallξ. This is understandable as we are dealing with two Fermi seas at high densities ∝ 1/ξ where the interaction is asymptotically negligible. The function E N l =∞ N +1 (ξ) has a local (for 1.8 < α < 2.4) or global (for α > 2.4) minimum at finiteξ. This minimum can correspond to the self-bound solution, as this is the point of zero pressure. The energy per cluster in this state is shown in the main panel of Fig. 1  (black solid curve). However, we should specify that the minimum of E N l =∞ N +1 (ξ) does not necessarily mean that this state is self-bound in free space. We have yet to check the conditions µ < 0 (otherwise the chain will lose heavy atoms) and max[ε p ] =ε π < 0 (otherwise it will lose light atoms). In fact, the right end point of the black solid curve in Fig. 1 corresponds toμ = 0. The behavior of the self-bound chain near this point is thus similar to what happens in the cases N l = 1 and N l = 2 already discussed. By contrast, the left end point does correspond to the disappearance of the minimum (see left panel in Fig. 4). There, the chain breaks into free N + 1 clusters. We should mention that by assuming one light atom per period we disregard possible breaking of a long chain into smaller clusters sufficiently close to α = 2.4 as we have pointed out in Sec. 3.
Although there is no apparent small parameter that allows us to solve Eqs. (7) and (8) perturbatively, the assumption of weak modulation, which provides a completely analytic description of the system, turns out to work extremely well for all values of α andξ relevant for analyzing the self-bound regime. The weak-modulation theory is based on the ansatz n(ū) = 1 + (δ/ξ) cos(2πū), where δ is assumed to be small. One then calculates E N l =∞ N +1 (ξ) given by Eq. (6) up to terms ∝ δ 2 , and minimizes it with respect to δ.
We calculate the spectrumε p following the standard weak-modulation approach [33]. Namely, substituting the expansion χ p = ∞ j=−∞ β j e i2πjū into Eq. (7), we get the set of equations for all integer j. Sinceε p =ε −p we can consider only the positive half of the first Brillouin zone, i.e., 0 < p < π. One can then check that the solution of Eqs. (10) for the lowest band is characterized by the following hierarchy. The coefficients β 0 and β −1 are of order one or smaller, β 1 and β −2 are ∼ δ or smaller, β 2 and β −3 are ∼ δ 2 , etc. Therefore, up to the second order in δ we can write β 1 = (δ/2)β 0 /[(p + 2π) 2 − p 2 ], where we use Eq. (10) with j = 1 and neglect the small difference (at most ∝ δ) betweenε p and the unpertubed energy p 2 −ξ. In a similar way, Eq. (10) with j = −2 gives Substituting these β 1 and β −2 into Eq. (10) with j = 0 and j = −1 we obtain Solving this linear system gives the spectrumε p with the desired accuracy.
To integrateε p we divide the p > 0 part of the Brillouin zone into two regions. The first is 0 < p < π − √ δ. Here we just use the Taylor expansion ofε p up to terms ∝ δ 2 . In the remaining interval π − √ δ < p < π we cannot Taylor expandε p as this would lead to a divergent integral. However, replacing p by π in the terms proportional to δ 2 in Eqs. (11) and (12) (one can check that this is a legal approximation in the considered integration interval) givesε p in the form suitable for analytic integration. The result is π −πε p dp 2π Finally, using 1 0n 3 (ū)dū = 1 + (3/2)δ 2 /ξ 2 , the minimization of Eq. (6) with respect to δ gives δ = 16π 2 e −24π 2 α/ξ 2 (14) and To complete the theory, we note that the chemical potential can be determined by raising Eq. (8) to the second power and by integrating the result overū. We obtainμ = 3α[1 + δ 2 /(2ξ 2 )] −ξ.
In Fig. 4 we show E N l =∞ N +1 (ξ) in units of the single-cluster energy |E N l =1 N +1 | as a function of the modulation period for α = 1.85, 2.4 and 4.0. The black solid curves are determined by exactly solving Eqs. (7) and (8) and the gray dotted curves are given by Eq. (15). The weak-modulation approximation is very precise (much less than a percent deviation from the exact numerics) up to the minima for all considered α. An appreciable difference can be seen only at largeξ, far from the minima.

Conclusion
In conclusion, we show that two N + 1 clusters formed in a two-component fermionic mixture can attract each other by a peculiar potential with a minimum at a finite intercluster separation. This attraction persists in the thermodynamic limit such that the Fermi-Fermi mixture can become self-bound forming a polymer of N + 1 clusters. One can also think of this state as a self-bound homogeneous liquid, undergone the Peierls charge density wave instability [29,34]. Since our theory is valid for N ≫ 1 [or (M/m) 1/3 ≫ 1] we cannot determine the smallest N [or M/m] at which N + 1 clusters bind. This problem should be tackled by other methods such as, for instance, exact diagonalization, quantum Monte-Carlo, or density matrix renormalization group. Our theory, taken at its face value, predicts binding of 4+1 clusters in the fermionic 173 Yb-6 Li mixture (M/m = 28.75), which can in principle be checked in current experiments [35,36].