Two impurities in a Bose-Einstein condensate: from Yukawa to Efimov attracted polarons

The well-known Yukawa and Efimov potentials are two different mediated interaction potentials. The first one arises in quantum field theory from the exchange of virtual particles. The second one is mediated by a real particle resonantly interacting with two other particles. This Letter shows how two impurities immersed in a Bose-Einstein condensate can exhibit both phenomena. For a weak impurity-boson attractive interactionattraction with the condensate, the two impurities form two polarons that interact through a weak Yukawa attraction mediated by virtual excitations. For a resonant attraction with the condensate, the exchanged excitation becomes a real boson and the mediated interaction changes to a strong Efimov attraction that can bind the two polarons. The resulting bipolarons turn into in-medium Efimov trimers made of the two impurities and one boson. Evidence of this physics could be seen in ultracold mixtures of atoms.

A particle interacting with a surrounding medium can form a polaron, i.e. it becomes dressed by a cloud of excitations of the medium that alters its properties. This general concept, introduced by Landau and Pekar [1] to describe electrons coupled to the vibrations of a lattice in solids, has proved useful to understand a variety of physical systems such as semi-conductors and superconductors [2]. In the last few years, polarons with arbitrarily strong interactions with the medium could be investigated experimentally using ultra-cold atoms [3][4][5][6]. These experiments have realised Fermi polarons (impurities embedded in a Fermi sea) by mixing different kinds of fermionic atoms and tuning their interaction by a Feshbach resonance. Recently, two experimental works [7,8] have reported the observation of Bose polarons (impurities embedded in a Bose gas) by using bosonic ultra-cold atoms. While the properties of a single Bose polaron are interesting and theoretically challenging [9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24], it is also of fundamental interest to understand the interaction between Bose polarons induced by their medium [24]. An exchange of bosonic excitations is known to induce a Yukawa potential between two polarons. This occurs in mixtures of bosonic and fermionic helium liquids [25].
A similar phenomenon appears in high-energy physics, where the nuclear force is mediated by mesons [26].
On the other hand, at the few-body level, it is known that for sufficiently strong interactions, an effective three-body force called the Efmov attraction can bind three particles into one of infinitely many three-body bound states, known as Efimov trimers [27][28][29]. The Efimov attraction can be understood as an interaction between two particles mediated by a third particle. It scales as the inverse square of the distance between particles, conferring discrete scale invariance to the system. Efimov trimers and their singular properties have been observed in ultra-cold atom experiments in the last few years [30], triggering the question of the influence of the surrounding medium on these trimer states. In a condensate of heavy bosons strongly interacting with light impurities, there is a strong Efimov attraction that can form Efimov trimers of two bosons and one impurity. Theoretical studies [23,31] have shown how a single polaron can turn into such an Efimov trimer, and a similar effect was found for an impurity in a two-component Fermi superfluid [32,33]. Reference [23] found that the in-medium Efimov trimer is stabilised by the surrounding condensate. In the opposite limit of heavy impurities in a condensate of light bosons, these trimers are very weak and the Efimov attraction favours instead the formation of trimers of two impurities and one boson. This indicates that in this case two polarons may turn into one such Efimov trimer. A study [34] has suggested that such an Efimov trimer would be weakened by the surrounding condensate. However, the theory could not completely describe the interaction at large distance between the two impurities. The precise effect of a surrounding Bose-Einstein condensate on Efimov trimers and the mediated interaction thus remain to be clarified.
Motivated by these theoretical questions and the recent experiments with ultra-cold atoms, this work presents a minimal description of two impurities in a Bose-Einstein condensate that bridges the perturbative regime of weakly attracted polarons and the nonperturbative regime corresponding to a bound Efimov trimer immersed in a Bose-Einstein condensate. In particular, the effective interaction between impurities is shown to go from the Yukawa type, mediated by virtual bosonic excitations, to the Efimov type, mediated by a real boson. This description is based on the method of Refs. [23,35] which uses a variational wave function for the impurities and the excitations of the medium. Here, the excitations are the Bogoliubov quasiparticles of the condensate. In the following, at most one excitation will be considered, which is the minimal requirement to reproduce the expected Efimov three-body physics. First, the mediated interaction between the two impurities will be derived, and then the energy spectrum of the system will be presented and discussed.
The Bose-Einstein condensate is assumed to be a homogeneous gas of bosons of mass m interacting via a weak pairwise interaction U B , whereas the interaction U between an impurity and a boson may be arbitrarily strong. No direct interaction between the impurities is considered. The impurities are assumed to be identical bosons of mass M . The Hamiltonian thus reads in arXiv:1607.04507v3 [cond-mat.quant-gas] 13 Apr 2018 second quantisation: where V is the system's volume, k = 2 k 2 2m and b k are the kinetic energy and annihilation operator for a boson with momentum k, and ε k = 2 k 2 2M and c k are the kinetic energy and annihilation operator for an impurity with momentum k. Since the bosons are weakly interacting, the first line of Eq. (1) can be approximately diagonalised where the operator β k annihilates a quasi-particle with momentum k, , where n 0 = N 0 /V is the condensate density. For convenience, the origin of energy is set to the condensate ground-state energy E 0 .
The total wave function |Ψ of the system can be expanded exactly as a superposition of any number of excitations on top of the Bose-Einstein consdensate |Φ and the two impurities. Truncating this expansion to at most one excitation gives the following ansatz, Applying the variational principle δΨ|H − E|Ψ = 0 to the Hamiltonian of Eq. (1) with the ansatz of Eq. (3), where α q and α q,q are varied independently, gives a set of two coupled equations: is the total density of bosons. Here, a B = m 4π 2 U B (0) is the boson scattering length in the Born approximation.
Let us first consider a weak interaction U , i.e. that can be treated perturbatively. This imposes that the Born expansion of the scattering length a = a 0 + a 1 + . . . converges rapidly, and a 0 = 2µ 4π 2 U (0) is much larger than is the boson-impurity reduced mass. In this case, one can neglect the sum in Eq. (5), as it contributes to higher orders in U . Let us now consider the limit M → ∞ of heavy impurities separated by a vector r, and perform a Fourier transform with respect to q, the conjugate momentum of r. Eliminating the second equation into the first, one obtains where E = E − 2nU (0). The solution E(r) of this equation as a function of r gives the effective potential between the two impurities in the Born-Oppenheimer approximation. Equation (6) shows that it decays as a Yukawa potential (see Appendix A.1), where ξ = (8πn 0 a B ) −1/2 is the condensate coherence length, and E(∞) = 8π 2 2µ na 0 + n 0 a 1 + √ 2a 2 0 /ξ is the asymptotic energy of the separated impurities, which is essentially twice the mean-field energy E M F = 4π 2 2m na of a single impurity [22]. This confirms well-known results for small scattering lengths [25,26,36]. Note that, due its variational nature, the potential E(r) is unreliable for r ξ, nevertheless it can be shown that the form of Eq. (7) remains correct for r → ∞ (see Appendix A. 1

and A.2).
To investigate the non-perturbative regime, let us now consider the limit of a contact interaction. It corresponds to a constant interaction in momentum space, i.e. U (k) = g < 0, up to some arbitrarily large momentum cutoff Λ. The scattering length a of this interaction is given by the relation which is used to renormalise all final results, i.e. express them in terms of the scattering length a instead of g. Using this interaction in Eqs. (4)(5), one encounters the terms F q = g 1 V p u p α q,p−q and G q = g 1 V p v p α q,p−q . Although F q remains finite when Λ → ∞, since the sum in its expression diverges as g −1 for a fixed value of a, the term G q vanishes, since the sum in its expression does not diverge. In the end, one finds the following equation (see Appendix B.1), (10) As previously, one can find the mediated interaction in the Born-Oppenheimer limit M → ∞ by fixing the distance r between the two impurities and performing the Fourier transform of these equations. One obtains This equation differs from Eq. (4) of Ref. [34] by its nonzero right-hand side and the coefficient u 2 k = 1. Let us consider its solution for weak (1/a → −∞), unitary (1/a = 0), and strong (1/a → +∞) boson-impurity interactions.
For weak boson-impurity interactions (a < 0 and |a| a B ), the solution E(r) of Eq. (11) is of the Yukawa type (see Appendix B.2), where It can be seen that this potential is slightly different from the perturbative result of Eq. (7). This is attributed to the limitation of the variational ansatz, which only provides an upper bound of the exact potential (see Appendix B.2).
Let us now consider the cases 1/a → 0 and 1/a → +∞.
It follows that in the large scattering length limit 1/a → 0, the mediated interaction E(r) has the form (see Ap- and L = (8πn 0 ) −1/3 . One recognises at short distances the 1/r 2 Efimov attraction (in the Born-Oppenheimer limit [29]) between two impurities mediated by a boson. The Efimov attraction can support an infinite number of bound states. However, here it is truncated at distances on the order of the mean boson spacing L = (8πn 0 ) −1/3 , and asymptotes to the energy E(∞) ≈ − 2 2mL 2 . As a result, the infinite number of possible trimer states in vacuum is reduced to a finite number, such that only those trimers whose energy is lower or comparable to E(∞) survive in the presence of the condensate.
Finally, as the boson-impurity interaction is strengthened towards small positive scattering length a, each polaron is expected to turn into a dimer of energy E d = − 2 2ma 2 , as each impurity should strongly bind with a nearby boson. However, in the present theory, the asymptotic energy E(∞) of the two separated impurities goes to E d instead of 2E d , as can be seen from Eq. (13) for r → ∞ and 1/a → +∞. The reason is that the ansatz of Eq. (3) includes only one bosonic excitation, and as a result only one impurity can bind with that excitation. A more quantitative treatment of the scattering threshold of the impurities in this regime would thus require at least two bosonic excitations.
Let us now turn to the energy spectrum of the system. In the contact model, the Efimov attraction exists at infinitely small distances, as seen in Eq. (14a), leading to the so-called Thomas collapse [28,37]. Some additional short-range scale is necessary to cure this problem and set the three-body observables [29]. This can be done at the two-body level by keeping a finite momentum cutoff Λ for the sum in Eq. (9). Alternatively, one may introduce a three-body force. The simplest way to introduce such a force is to set a momentum cutoff Λ 3 on the second Jacobi momentum, i.e. the argument of F in Eq. (9). In atomic gases, this three-body parameter is related to the van der Waals length of the atoms [29,[38][39][40]. Figure 1 represents the exact energy spectrum of the system for a mass ratio M/m = 19 (such as caesium-133 atoms in a lithium-7 condensate), as a function of 1/a, and calculated numerically from Eq. (9) with a three-body cutoff Λ 3 .
For any boson-impurity interaction, the spectrum shows a continuum corresponding to scattering states of two attractive polarons. Its threshold, shown by the dotted curve in Fig. 1, corresponds to the asymptotic limit of the mediated interaction, which for the large mass ratio used here is well approximated by the Born-Oppenheimer threshold given by the solution of Eq. (13) for r → ∞. As noted before, the threshold corresponds to the mean-field energy 2E MF of two polarons for small a < 0, and (unphysically) asymptotes to the energy E d of a single dimer for small a > 0. The spectrum also features discrete bound states for sufficiently strong boson-impurity interaction. This is expected since the mediated interaction becomes strong enough to bind the two polarons into bipolarons as it gradually turns from a weak Yukawa potential into a strong Efimov attraction. As the interaction further increases, the bipolarons (shown as solid curves in Fig. 1) turn into Efimov trimers made of two impurities and one boson. As anticipated from the Born-Oppenheimer potential between the two impurities -see Eq. (14) -only the trimers whose energy is lower than the polaron scattering threshold survive in the presence of the condensate. Near unitarity (1/a = 0), the bipolaron energies are pushed down from the trimer energies in vacuum due the attractive effect of the surrounding bosons, but the binding energies relative to the polaron scattering threshold are smaller than in vacuum. Near unitarity, the trimers are therefore weakened by the condensate. However, interestingly, the bipolarons exist for weaker boson-impurity interactions than the vacuum trimers. The scattering lengths at which the bipolarons appear (indicated by arrows in Fig. 1) are indeed reduced in magnitude with respect to vacuum. In this sense, the condensate favours the appearance of the trimers. This is especially true when the polaron scattering threshold at unitarity (∼ − 2 /(2µL 2 )) is comparable to the energy of an Efimov trimer in vacuum, as shown by the second bipolaron in Fig. 1. It should be mentioned that the boson-boson interaction has the opposite effect of weakening the bipolarons, but this effect remains small in the assumed dilute regime na 3 B

1.
Finally, it is important to specify the range of validity of the present treatement. As noted above, for sufficiently large attraction between the impurities and the bosons, more than one bosonic excitation are needed. Since the bipolarons are seen in Fig. 1 to correlate to the  , where Λ3 is a three-body cutoff inversely proportional to the range of the boson-impurity interaction. The spectrum is obtained from Eq. (9), which is expected to be quantitative for 1/a 0. The shaded area represents the scattering continuum of the two attractive polarons and its threshold is shown as a dotted curve. The solid curves correspond to the bound states (bipolarons). The points where they appear from the polaron scattering threshold are indicated by the vertical arrows. For reference, the black dot-dashed curve shows the boson-impurity dimer energy in vacuum, and the red dashed curves correspond to the boson-impurity-impurity trimer energies in vacuum.

Two-polaron spectrum
vacuum trimer states of two impurities and one boson, it is natural to expect that an additional bosonic excitation would correlate them to tetramer states of two impurities and two bosons. Such tetramers do exist for weak interactions between the bosons and would significantly affect the two-polaron spectrum beyond the results presented here. However, for a moderate interaction between the bosons a B Λ −1 3 that is typical for ultracold atoms, the tetramers are found to be suppressed for 1/a 0 [41]. It is therefore expected that the present theory is quantitative in this regime.
In summary, a simple variational ansatz has been used to investigate the problem of two impurities in a Bose-Einstein condensate. The ansatz bridges the well-known perturbative regime to the non-perturbative regime, where the Bose-mediated interaction takes the form of the Efimov attraction. It shows that the two polarons formed by the two impurities merge into one or several Efimov trimers for sufficiently strong interaction. The stability of these bipolarons under the influence of the condensate has also been revealed. Although their binding energy is reduced near unitarity with respect to that of trimers in vacuum, they exist for smaller interaction as the density of the condensate is increased. In a mixture of resonantly interacting ultra-cold atoms, this would appear as a boson-density-dependent shift of the three-body loss peaks associated with the appearance of Efimov trimers. The direct effect of the mediated interaction between impurities could be observed as an impurity-density-dependent mean-field shift (estimated on the order of a few percents) in the single-polaron energy spectrum.
The author thanks Nguyen Thanh Phuc, Takumi Doi, and Tetsuo Hatsuda for helpful discussions. This work was partially supported by the RIKEN Incentive Research, iTHES, and iTHEMS projects.. The Yukawa potential of Eq. (7) is obtained as follows. For sufficiently small scattering length, the term E in the denominator of Eq. (6) may be neglected as it contributes to higher orders, and the term (u k − v k ) 2 = k /E k . It follows that: Using k = 2 k 2 /(2m) and U B (0) = 4π 2 a B /m, one gets The term E (∞) does not depend on r and can be written as a follows: The first sum in the above expression converges due to the decay of U (k) at large k, and is related to the second term a 1 = − 2m in the Born expansion of the scattering length a. In the second sum, one may take 2m a 0 since the momentum range of U is typically much larger than ξ −1 . One then obtains The last sum in Eq. (17) goes to zero as r → ∞. Its asymptotic behaviour at large r may be obtained from the low-momentum contribution in the sum. In this limit, U (k) may be approximated by U (0), i.e.

Shortcoming of the variational ansatz
We should note that while Eq. (19) gives the form of E (r) in the limit of small interaction, the convergence to this potential is not uniform. For small but finite attraction U , there is indeed a distance beyond which the term E in the denominator of Eq. (6) may not be neglected. Taking into account this term, one can find the true asymptotic form to be , and the facts that F (0) = 0 and F (2) One can show that this 1/r 4 repulsion occurs for r ξ 2m 2 E (∞) −1 . However, this asymptotic behaviour has no physical reality as it is an artifact of the variational ansatz of Eq. (3) for large distances. Indeed, the variational ansatz only gives an upper bound of the exact potential and its scattering threshold. Even if the variational threshold is only slightly above the exact one, the way the variational potential asymptotes to this threshold may be completely different from the way the exact potential asymptotes to the exact threshold. As a result, the analytic form of the asymptote may be wrong, as seen here.

Exact derivation
It turns out that the perturbative result Eq. (19) holds exactly beyond the variational ansatz of Eq. (3) and to infinite distances. To show this, let us first treat the boson-impurity as a perturbation to first order in the Hamiltonian. We obtain the following Frölich-like Hamiltonian: where g k = U (k) √ n 0 (u k −v k ). Taking the limit of static impurities (M → ∞), we get ε k → 0 and write p c † p+k c p = d 3 Rn c (R)e ik·R , where n c (R) is the density of the impurities. For two impurities separated by r, we have n c (R) = δ 3 (R) + δ 3 (R + r). This gives the r-dependent Hamiltonian, We now note that the Hamiltonian can be diagonalised exactly by introducing the operatorβ k such that One can check that this operator satisfies the bosonic commutation relations [β k ,β q ] = 0 and [β k ,β † q ] = δ k,q , and the Hamiltonian expressed in terms of this operator reads The ground state is therefore given by which is exactly the same as Eq. (15). It follows that the results Eqs. (18) and (19) are exact in the perturbative limit, beyond the variational ansatz Eq. (3).

Derivation of the equation
The following provides the derivation of Eq. (9). Starting from the general equations Eqs. (4)(5), one performs the changes α q−k,−q = α −q,q−k in Eq. (4) and α q−p−k,k−q = α k−q,q−p−k in Eq. (5), using the bosonic exchange symmetry α q,q = α q ,q . Then, one sets the potential U (k) = g for k < Λ, U (k) = 0 for k ≥ Λ. This yields In the first equation, one can change the term α −q,q−k into α −q,k+q by performing the change of variable k → −k. The equation then reads, where the terms F q and G q are defined by Next, the change of variable p → −p is performed in Eq. (21). For sufficiently large Λ, the sum |p−k|<Λ p can be approximated by p<Λ p , so that Eq. (21) can be expressed in terms of F and G: Using this equation to express α q,k−q in Eqs. (23) and (24), one finds Owing to the renormalisation relation Eq. (8), for a fixed scattering length a, the term 1/g in the left-hand side of Eqs. (26) and (27) diverges as Λ for very large Λ. In Eq. (26), this divergence in the left-hand side is cancelled by another divergent term in the right-hand side. In contrast, in Eq. (27), the right-hand side does not diverge for large Λ. Indeed, for large k, the denominator in Eq. (27) is ∼ k 2 , and the numerator involves the terms v k u k ∼ k −2 and v 2 k ∼ k −4 . The term in the sum of Eq. (27) thus decay as k −4 or faster, and the sum is therefore convergent. One concludes that G may be neglected for sufficiently large Λ.
There only remain two equations, Eq. (22) and (26), which for G = 0 read One can further rewrite Eq. (29) as The sum in the last term is convergent, since for large k, u k (u k − v k ) ∼ 1 and α q−k k −2 according to Eq. (28). Since it is multiplied by the vanishing factor g, it can therefore be neglected. On the other hand, the sum in the second term of Eq. (30) diverges as −1/g -as seen from Eq. (8). It therefore cancels the factor g for large enough Λ. Finally, the divergence of the sum in the left-hand side of Eq. (30) is cancelled by the term 1/g, which can be done explicitly by using the renormalisation relation Eq. (8). Finally, Eq. (30) simplifies to The energy E = E − 2ng may be replaced by E since g vanishes for large Λ. As a last step, one can perform the changes q → −q and k → −k in Eq. (31), and observe that F −q satisfies the same equation as F q . One can thus set F −q = F q in Eq. (28). Combining this equation with Eq. (31) finally yields Eqs. (9), where the limit Λ → ∞ is taken.
It is worthwhile to note that the same equation can be obtained in a different way from a two-channel contact model, whose range parameter is set to zero.

Derivation of the mediated potentials
The following provides the derivation of the Born-Oppenheimer potentials Eq. (12) and Eq. (14) in the limit of small scattering length (1/a → −∞) and unitarity (1/a → 0).
The Born-Oppenheimer equation (11) may be written as follows, which can be further expressed as Small scattering length In the limit of small scattering length 1/a → −∞, the energy E goes to zero. Therefore, in this limit one can neglect the term E(r) in the denominator in the above equation. The integral can then be calculated analytically, yielding the following explicit expression for E(r), where F (r) = [1 + e − √ 2r/ξ + 2I 0 ( √ 2r/ξ) − 2L 0 ( √ 2r/ξ)]/(4r), and I 0 is the modified Bessel function of the first kind, and L 0 denotes the modified Struve function. One can check that F (r) ≈ exp(− 4+π 2 √ 2π r ξ )/r for r ξ and F (r) ≈ 1/(4r) for r ξ. However, similary to the derivation of Appendix A, the convergence to Eq. (35) for small negative a is not uniform. At very large distance, one may not neglect the term E(r) in the denominator, and one obtains instead F (r) = 1 2 E(∞)| , giving the asymptotic behaviour, which can be shown to hold in the present situation for r ξ −1 /| 2m 2 E(∞)|. As noted in Appendix A, because of the variational nature of the calculation, the form of E(r) at large distance is not physical, so we restrict our consideration to r ξ. Treating a as a small perturbation, one obtains which yields Eqs. (12). As noted in the main text, this potential is close but somewhat different from the perturbative result of Eq. (7). If one slightly worsens the variational potential by replacing u k by u k − v k in Eq. (33), one obtains a form that is closer to Eq. (7), Unitary limit For larger scattering lengths, the term E becomes larger than 2 ξ −2 2m , so that the integral in Eq. (34) is mostly determined by E(r) and we can treat ξ −2 as a perturbation. To first order in ξ −2 ,we get: where we set E = − 2 κ 2 /(2m). The integration can then be performed, yielding E(r) = 8π 2 n 0 2m or equivalently, from which Eq. (13) is obtained. Now we take the unitary limit 1/a → 0 in Eq. (41). For small r, the term 1/r dominates over the other terms except κ so that the equation reduces to −κr + e −κr = 0. (42) It can be checked that the missing terms can be neglected when r L (and r ξ which is readily satisfied since L ξ). The solution of Eq. (42) is κr = W (1), where W is the Lambert function, leading to the result of Eq. (14a). For large r, one can first neglect the r-dependent terms since they vanish. One obtains the threshold value κ ∞ satisfying the equation This cubic equation in κ ∞ admits the solution with X = 1 2 L −3 + L −6 + 4 27 ( 1 2 ξ −2 ) 3 = L −3 + O(ξ −6 ). Treating the r-dependent terms of Eq. (41) as a perturbation, one can set κ(r) = κ ∞ (1 + ε(r)) in this equations, which yields ε(r) = 3κ ∞ + ξ −2 2κ∞ −1 1 r − ξ −2 2κ∞ e −κ∞r . This gives, to first order in ξ −2 , the result of Eq. (14b). The smallness of ε(r) 1 requires that r κ −1 ∞ ∼ L. Again, we should note that according to Eq. (36) one finds that for much larger distances (in this case for r ξ) the potential approaches its threshold as 1/r 2 , although this is an unphysical artifact of the variational ansatz.

Numerical solution
In the s-wave channel (F q = F q ), where the Efimov attraction takes place, the equation (9) simplifies as follows, where a three-body momentum cutoff Λ 3 has been imposed on the argument of F . Making the substitution 1 V k ≈ (2π) −3 d 3 k and using the explicit forms of T k (E), E k and ε k gives