Heavy-light N + 1 clusters of two-dimensional fermions

We study binding of N identical heavy fermions by a light atom in two dimensions assuming zero-range attractive heavy-light interactions. By using the mean-field theory valid for large N we show that the N + 1 cluster is bound when the mass ratio exceeds 1.074 N 2 . The mean-field theory, being scale invariant in two dimensions, predicts only the shapes of the clusters leaving their sizes and energies undefined. By taking into account beyond-mean-field effects we find closed-form expressions for these quantities. We also discuss differences between the Thomas-Fermi and Hartree-Fock approaches for treating the heavy fermions


Introduction
Binding in the fermionic N + 1-body model with zero-range interactions is a fundamental problem, necessary for general understanding of fermionic mixtures with mass and population imbalance.Apart from the interspecies scattering length a > 0, which determines the size of the 1+1 cluster, the model is parametrized by the number of heavy fermions N , the mass ratio M/m, and the space dimension D. This parameter space has unexplored spots in spite of the constant interest to the problem from the nuclear-physics side and, more recently, from the ultra-cold-gas community.
In contrast to attractive bosons, which typically always bind, fermionic N + 1-clusters bind only above a critical mass ratio, such that the interspecies attraction overcomes the Fermi pressure.Previous studies have shown that binding of larger clusters requires higher mass ratios or lower dimension, which can be explained by the dependence of the Fermi pressure (kinetic energy of noninteracting fermions) on M and D. The current status of the fermionic N + 1 problem is as follows.
In three dimensions, the 2 + 1 trimer binds at M/m = 8.2 [1], the 3 + 1 tetramer at M/m = 8.9 [2,3], and the 4 + 1 pentamer at M/m = 9.7 [3].A qualitative picture which explains the relatively small spread in the critical mass ratios for these bound states, as well as their angular momenta and parities, is that the dimer acts as a p-waveattractive scattering center for heavy atoms and thus accomodates three orbitals with different projections of the angular momentum [3].These N + 1-body systems become Efimovian above certain mass-ratio thresholds in the vicinity of M/m ≈ 13 and require a three-body, four-body, and five-body parameter, respectively [3][4][5].
In one dimension there is no Efimov effect, no fermionic sign problem, and no shell effects.The 2 + 1 trimer exists for any M/m > 1 [6] and the mass-ratio thresholds of N + 1 clusters steadily grow with N [7].In the limit of large N their shapes and energies are well described by the mean-field theory; in this limit the critical mass ratios scale as M/m = π 2 N 3 /36 [7].
In two dimensions the atom-dimer attraction in the p-wave channel leads to a formation of the 2 + 1 trimer with unit angular momentum for M/m > 3.33 [8].The 3 + 1 tetramer emerges almost immediately, for M/m > 3.38 [9], pointing to even stronger shell effects than in three dimensions.Exact calculations demonstrate the presence of an excited tetramer for M/m > 5 [10] and a ground-state pentamer for M/m > 5.14 [9].Although there is no Efimov effect, the fermionic statistics and rapid growth of the configurational space with N makes the analysis of larger clusters technically difficult.
In this paper we address the large-N limit of the two-dimensional N + 1-body problem by applying the mean-field (MF) approximation together with the local-density Thomas-Fermi (TF) assumption for the energy of an ideal Fermi gas.We show that the N + 1 cluster emerges for M/m = 2N 2 /C, where C = 1.862.At this critical point the cluster shape is controlled by the nonlinear Schrödinger equation with attractive cubic nonlinearity like in the case of attractive two-dimensional bosons.We thus notice similarities of the critical N + 1 cluster with the Townes soliton [11], recently observed in ultracold bosonic atoms [12,13].For larger mass ratios our theory gives the shape of the cluster as a function of the parameter α = 4πN 2 /(M/m) < 2πC and the coupling constant for which the solution exists.However, being scale invariant this theory does not predict the size or energy of the cluster.We show that the leading beyond-MF correction for the N +1 cluster with N ≫ 1 is a local quantity and we use it to calculate the cluster energy and size, which both scale exponentially with N .We show that replacing the TF approximation by the full Hartree-Fock treatment is necessary for determining the preexponential factor.
2 Mean-field Thomas-Fermi approach for large N We consider a mass-imbalanced fermionic mixture governed by the Hamiltonian where φ † r and Ψ † r are the creation operators of light and heavy fermions, respectively, and we set h = 1.The short-range heavy-light interaction is characterized by the coupling constant g = 2π/[m r ln(2m r |E 1+1 |/κ 2 )] < 0, where m r = mM/(m + M ) is the reduced mass, E 1+1 is the dimer energy, and κ is the ultraviolet cut-off momentum assumed to be much larger than any other momentum scale in the problem.The dimer energy is related to the heavy-light scattering length by E 1+1 = −2e −2γ E /(m r a 2 ), where γ E is the Euler constant.
We write the MF energy functional for the N + 1 system as where ϕ(r) is the wave function of the light atom, the product N n(r) is the density profile of the heavy atoms, and we introduce two dimensionless parameters: α = 4πmN 2 /M and γ = 2mgN < 0. Equation ( 2) is valid for weak interactions, i.e., m r |g| ≪ 1.The term ∝ n 2 (r) is the kinetic energy density of an ideal Fermi gas taken in the TF local-density approximation valid for N ≫ 1, when n changes slowly on the mean interparticle distance.
To minimize Eq. ( 2) with the normalization constraints |ϕ(r)| 2 d 2 r = 1 and n(r)d 2 r = 1 we introduce the Lagrange multipliers ϵ and µ and minimize the grand potential Ω = E − [µN n(r)+ϵ|ϕ(r)| 2 ]d 2 r.The conditions δΩ/δϕ = 0 and δΩ/δn = 0 lead to the coupled equations where θ(x) = (x + |x|)/2.Equation (3) describes a light atom in an effective well formed by the MF attraction of the heavy fermions.Similarly, Eq. ( 4) is the TF density profile of the heavy fermions in the MF trap created by the light atom.The Lagrange multipliers ϵ and µ have physical meanings of the energy of the light atom and the chemical potential of the heavy atoms, respectively.For self-bound solutions of Eqs.(3) and (4) these quantities should both be negative since ϕ and n are not allowed to spread over the whole space.The binding threshold for the N + 1 cluster corresponds to µ = 0, when the heavy atom at the Fermi surface is nearly unbound.In this case Eq. ( 4) reduces to n(r) = −γ|ϕ(r)| 2 /α.The normalization constraints then imply −γ = α and Eq. ( 3) becomes the nonlinear Schrödinger equation with negative cubic nonlinearity like in the case of attractive two-dimensional bosons.The solution of this problem is the Townes soliton which exists only for a specific value of the coupling constant [11].In our case, this compatibility condition reads −γ = α = 2πC = 11.7 and the critical wave function is given by ϕ [11,14].The scale R is arbitrary and cannot be determined from the MF set of Eqs. ( 2)-( 4).Interestingly, the MF energy (2) vanishes for this solution independent of R [15][16][17][18].For bosons this problem is solved by the fact that the renormalized coupling constant g r depends logarithmically on R [14] leading to a shallow (beyond-MF) minimum in E(R) at a certain R. Here, for fermions, the stabilization mechanism is similar, but as we will show below, the beyond-MF contribution has a slightly different form and can be calculated in the local-density approximation.
By numerically solving Eqs. ( 3) and ( 4) above the critical mass ratio, i.e., for α < 2πC, we always find a radially symmetric real nodeless self-bound solution.More precisely, each α gives rise to a family of self-similar solutions which exist only for a certain γ = γ c (α).We parametrize this family by the length scale R = 1/ √ −2mϵ.Formally setting ϵ = ϵ 0 = −1/(2m), Eqs. ( 3) and ( 4) become adimensional and we denote their solution by ϕ 0 (ρ), n 0 (ρ), and µ 0 .Then, for any R > 0 the dimensional solution for the same α and γ All these solutions of Eqs. ( 3) and ( 4), for any α and any R, correspond to vanishing E. Physically, it follows from the fact that Eq. ( 2) scales with R as E ∝ R −2 .Then, if E ̸ = 0, the system would shrink or expand, contradicting the stationarity of the found solution.That the stationarity of a two-dimensional soliton with cubic (scale-invariant) nonlinearity is equivalent to E = 0 has been mathematically shown in Ref. [15] (see also Refs.[16][17][18]).In our case, to make sure that from Eqs. ( 3) and ( 4) indeed follows E = 0 we derivate Eq. ( 3) with respect to R using the scaling properties of ϕ(r) and n(r) mentioned above.We then eliminate ϵ from the result by employing the same Eq.( 3) again.We then multiply the resulting equation by ϕ and integrate it over space obtaining the equality [−2ϕ(r)∇ 2 ϕ(r) − γϕ 2 (r)rn ′ (r)]d 2 r = 0, which can further be reduced to the form E = 0 with the help of Eq. ( 4).
We emphasize that γ c is not an external parameter but a characteristic of the solution of the MF Eqs. ( 3) and ( 4).On the other hand, γ (or g) is what we are allowed to tune.Although for γ ̸ = γ c the MF solution is not stationary, it may become stationary once we take into account beyond-MF effects.The final result should not depend on the choice of g and κ as long as a is fixed (see Sec. 3).

Beyond-mean-field correction
We see that the MF analysis cannot predict the sizes and the binding energies of the clusters, although it does predict their shapes (up to the rescaling) and determines the threshold mass ratio M/m = 2N 2 /C.Since the beyond-MF correction is not scale invariant, it introduces preferred length and energy scales, which can be understood from the following arguments.In two dimensions the second-order correction to the energy of two atoms interacting via a delta potential is logarithmically diverging at high momenta.Therefore, the beyond-MF correction to Eq. ( 2) is dominated by the renormalization of the two-body coupling constant, logarithmic in κ.It is thus convenient to express the beyond-MF-corrected energy by writing Eq. ( 2) with g replaced by g r = g + δg, where One can check that the renormalized coupling constant g r = g + δg is cut-off independent up to the second-order terms in the small parameter m r |g| ≪ 1.This renormalization removes the cut-off dependence from the energy to this order.The physical (i.e., cut-off independent) part of the beyond-MF contribution is absorbed into the length scale ξ, which is a functional of the fields ϕ and n, in general nonlocal.Qualitatively, 1/ξ is the characteristic momentum governing the many-body or few-body problem at hand.For two atoms in a box ξ is proportional to the box size.For a weakly interacting uniform Bose gas ξ is proportional to the healing length and this result can also be applied in the inhomogeneous case, if the density varies slowly on the scale ξ (see, for instance, [19]).On the other hand, for attractive bosons the local-density approximation does not work since ξ is proportional to the soliton size.However, this very fact that ξ ∝ R leads to important predictions for the energy and size scalings of bosonic solitons [14].
In our fermionic N +1 case the typical second-order process contributing to the beyond-MF term is a virtual excitation of the light atom creating a particle-hole excitation in the Fermi sea of heavy atoms.The typical momentum transfer is on the order of the Fermi momentum, which means that ξ is comparable to the mean interparticle separation for the heavy atoms, which scales as R/ √ N .Therefore, the beyond-MF correction in our case is local and can be obtained by analyzing the homogeneous problem.We just need to know the second-order ground-state energy shift for a single light atom immersed in a uniform Fermi sea of heavy atoms with Fermi momentum p F .Having found no answer in the literature we briefly outline this calculation.
Normalizing the single-particle states per unit surface we write the second-order energy correction as [20] ∆E (2) The integration domain in Eq. ( 7) is defined by the inequalities p 2 < p F , p < κ, and |p 2 − p| > p F , corresponding to the following virtual process.The unperturbed state is the impurity at rest and a Fermi sea filled up to p F .The virtually excited state is the light atom at momentum p, a heavy hole at momentum p 2 , and a heavy atom at momentum p 2 − p.We find it convenient to expand p 2 into a vector p ∥ parallel to p and a vector p ⊥ perpendicular to p.The integral over the angle of p gives 2π.We then integrate Eq. ( 7) over p ⊥ , then over p, and finally over p ∥ .In this manner, neglecting finite-range corrections p 2 F o(p F /κ), we obtain where Since we are interested in the regime M/m ∼ N 2 ≫ 1, Eq. ( 9) further simplifies to ξ = e 1/2 /p F which confirms the qualitative guess ξ ∼ 1/p F and shows that the large mass ratio does not dramatically influence this estimate.We can now proceed with the localdensity approximation.Substituting p F = 4πN n(r) into Eqs.( 8) and ( 9), keeping only the leading-order terms at large M/m, and going back to notations of Eq. ( 2) we write the beyond-MF correction to the cluster energy in the form Note that E BMF ∝ N −1 ln N is smaller than any of the three terms in Eq. ( 2), which are ∼ 1 (for a cluster of unit size).Therefore, Eq. ( 10) cannot strongly influence the shape of the cluster, but it can remove the degeneracy related to arbitrariness of R. Substituting ϕ(r) = ϕ 0 (r/R)/R and n(r) = n 0 (r/R)/R 2 into Eqs.( 2) and (10) and assuming γ = γ c + O(1/N ) we obtain up to the terms of order 1/N where . Minimization of Eq. ( 11) then gives and The parameters g and κ drop out from Eqs. ( 12) and ( 13) consistent with the fact that the length scale of the problem is given only by the scattering length and the energy scale by the energy of the 1 + 1 molecule.The preexponential-factor accuracy in Eqs. ( 12) and ( 13) follows from the beyond-MF accuracy of Eq. ( 11) and from the fact that the weak-interaction parameter m r |g| ≪ 1 is equivalent to 1/N ≪ 1 in the self-bound regime since γ ∼ 1.We note, however, that the TF approximation for the kinetic energy of the heavy fermions is guaranteed only to the leading order in 1/N .To estimate the error we pass to the Hartree-Fock (HF) description.

Hartree-Fock approach
The Hartree-Fock approach for solving the N + 1 cluster problem consists in introducing N orthonormal orbitals Ψ ν (r) and minimizing Eq. ( 2), in which the TF approximation αn 2 (r)/(4m) is replaced by N ν=1 |∇Ψ ν | 2 /(2M ).With these notations, the minimization of the energy functional with respect to ϕ gives Eq. (3) with and the minimization with respect to the orbitals Ψ ν (r) leads to where ω ν are Lagrange multipliers corresponding to the normalization constraints |Ψ ν | 2 d 2 r = 1.The functions ϕ and n determined by the TF and HF approaches are different, but one can easily check that their scaling properties are the same.In the HF method the length scale can thus also get fixed by setting 2mϵ 0 = −1 and denoting the corresponding solutions by ϕ 0 and n 0 .In addition, we assume cylindrical symmetry by imposing ϕ 0 (r) = ϕ 0 (r).Equation ( 15) then splits into one-dimensional Schrödinger equations for functions ψ l,ν (r), such that Ψ ν (r) = ψ l,ν (r)e ilφ , with l being the integer angular momentum.States with l ̸ = 0 are doubly degenerate corresponding to ±l.To find the ground state, we diagonalize Eq. ( 15) and we select the N states with the lowest energy among all possible channels.We then plug these states into Eq.( 14), substitute n 0 into Eq.( 3), and find γ for which Eq. ( 3) has a ground state corresponding to 2mϵ 0 = −1.The function ϕ 0 is then substituted back into Eq.( 15) and the process is repeated until convergence.This procedure results in the critical γ HF c which, in contrast to γ TF c (we use superscripts to specify the method), does not only depend on α but also on N .We show results up to N = 256.The black straight lines indicate the slope N −1 .We believe that stronger fluctuations for larger α are due to the fact that the uppermost filled heavy-atom orbitals are closer to the dissociation threshold and the system is thus more sensitive to changes in N .Accordingly, we find that the case of larger α and N requires finer and larger spatial grids for accurate calculations.
The fact that |γ HF c − γ TF c | scales as N −1 on average shows phenomenologically that to keep up with the claimed accuracy (up to the preexponential factor) for the cluster energy we should use γ HF c instead of γ c in Eqs.(12) and (13).The best TF-based prediction for the cluster energy is therefore E N +1 = −a −2 e −4πN/γc−2 ln N +O(N 0 ) .Advantages of the TF approximation is that it has α as the only input parameter and that limiting cases can be worked out analytically.By contrast, the higher accuracy of the HF method comes at the price of doing separate calculations for each N and M/m.However, the HF method predicts the structure of the cluster, its angular momentum and parity.It can also naturally handle excited states.
5 Hartree-Fock method applied to small clusters Although the HF method is valid for N ≫ 1, it is tempting to benchmark its performance for the few lowest-order clusters, for which exact results are known [8][9][10].In addition, the obtained solutions can also be used as guiding functions for the fixed-node Monte-Carlo scheme (see, for instance, [21]).We find that the HF method describes these small clusters rather well and we expect the accuracy to further improve with increasing N .
In Fig. 3 we plot the energies E HF N +1 in units of the exact dimer energy as a function of the mass ratio.The different curves stand for the 1 + 1 cluster (dotted black, in the HF description the heavy atom occupies the lowest s-wave orbital), 2 + 1 trimer (solid black, occupied are the lowest s-wave and one of the two degenerate l = ±1 orbitals), 3 + 1 ground tetramer (dotted red, occupied are the lowest s-wave and both lowest p-wave orbitals), 4 + 1 pentamer (dash-dotted blue, occupied the lowest and the first excited swave and both p-wave orbitals), and the excited 3 + 1 tetramer found in Ref. [10] (dashed brown, occupied are the lowest and the first excited s-wave and one of the lowest pwave orbitals).One can see that the HF approach reproduces the structure of the levels rather well (cf.[9]), although the artifacts of the approach are also visible.For instance, the trimer and the tetramer emerge immediately with finite binding energies, which is a consequence of the nonlinearity of the equations.The threshold behavior depends on the angular momentum of the orbital or, more precisely, on the convergence properties of the corresponding normalization integral.Since at zero energy the orbitals with angular momentum l behave as ψ l,ν ∝ r −|l| , the normalization integral for s-wave orbitals diverges.Therefore, a heavy atom in the s-wave orbital is in the halo state and does not influence the core.The crossing is therefore smooth (see the crossings of the pentamer and the excited tetramer).By contrast, for |l| > 1 the zero-energy orbital function is normalizable that the newly bound heavy atom is inside the core right at the threshold.This creates artifacts due to the nonlinearity.We find that the case |l| = 1, in spite of the logarithmic divergence of the normalization integral, is also prone to this nonlinear effect.
The cylindrical-symmetry assumption has to be carefully checked, but it requires a more involved two-dimensional analysis.We leave this task as well as the investigation of higher-order clusters to the future.

Conclusion
A two-dimensional fermionic N + 1 cluster binds for sufficiently large M/m.The MF theory valid for large N predicts the threshold value M/m = 2N 2 /C = 1.074N 2 and the cluster shape at this point and for larger M/m.The beyond-MF analysis based on the local-density approximation gives closed-form expressions for the size and energy of the cluster.The accuracy and practical relevance of the obtained results can be increased by switching to the Hartree-Fock form of the MF density functional.Finally, our findings have implications for ultracold fermionic mixtures.We can think of strongly mass-imbalanced mixtures of 6 Li with 173 Yb [22,23] or with other heavy Lanthanides such as Dy or Er.

Figure 2
Figure2shows the convergence of the HF value γ HF c (α, N ) towards the TF value γ TF c (α) at large N for α = 2 and α = 8.We show results up to N = 256.The black straight lines indicate the slope N −1 .We believe that stronger fluctuations for larger α are due to the fact that the uppermost filled heavy-atom orbitals are closer to the dissociation threshold and the system is thus more