Impurity in a BoseEinstein condensate: study of the attractive
and repulsive branch
using quantum MonteCarlo methods
Abstract
We investigate the properties of an impurity immersed in a dilute Bose gas at zero temperature using quantum MonteCarlo methods. The interactions between bosons are modeled by a hard sphere potential with scattering length , whereas the interactions between the impurity and the bosons are modeled by a shortrange, squarewell potential where both the sign and the strength of the scattering length can be varied by adjusting the well depth. We characterize the attractive and the repulsive polaron branch by calculating the binding energy and the effective mass of the impurity. Furthermore, we investigate structural properties of the bath, such as the impurityboson contact parameter and the change of the density profile around the impurity. At the unitary limit of the impurityboson interaction, we find that the effective mass of the impurity remains smaller than twice its bare mass, while the binding energy scales with , where is the density of the bath and is the common mass of the impurity and the bosons in the bath. The implications for the phase diagram of binary BoseBose mixtures at small concentrations are also discussed.
I I. Introduction
The polaron problem is a paradigmatic topic in condensed matter physics: it concerns the effect of the quantum fluctuations of the surrounding medium on the properties of an impurity immersed in a bath. The first formulation of the problem is due to Landau and Pekar Landau46 () in their study of the motion of electrons in polar crystals. By using a variational approach valid in the limit of strong coupling it was shown that the electron becomes eventually trapped in the potential created by the selfinduced deformation of the lattice. Fröhlich Froehlich54 () proposed an effective Hamiltonian, which describes the coupling between charged impurities and longitudinal optical phonons of the lattice, providing the standard model description of the polaron problem. The groundstate energy of the Fröhlich Hamiltonian was calculated by Feynman Feynman55 () within a variational approach based on path integrals which yields an upper bound that nicely interpolates between the weakcoupling perturbative result and the strongcoupling LandauPekar prediction. After a large amount of theoretical work Mahan90 () spanning more than four decades, an exact solution of the model by means of diagrammatic Monte Carlo methods was finally presented in Refs. Prokofev98 (); Prokofev00 (). Remarkably, both the results obtained for the polaron groundstate energy and the effective mass are in very good agreement with Feynman’s findings.
Important generalizations of the Fröhlich Hamiltonian include the Holstein model on a lattice Devreese09 (), electrons coupled to acoustical phonons in a crystal Peeters85 () and, more recently, impurity particles immersed in a dilute BoseEinstein condensed (BEC) gas Cucchietti06 (); Kalas06 (); Jaksch08 (); Tempere09 ().
In the context of ultracold atoms the polaron concept has received large attention in its fermionic version, i.e. an impurity coupled to a Fermi sea Massignan14 (). Thanks to the use of a Feshbach resonance the wave scattering length between the impurity and the fermions of the bath can be tuned at will and experiments have probed various properties of attractive and repulsive Fermi polarons both in 3D Zwierlein09 (); Salomon09 (); Grimm12 () and in 2D Koehl12 (), ranging from the weak to the strongcoupling regime and including the polaronmolecule transition.
Bose polarons, which involve a bath consisting of a BEC and therefore are directly related to the original Fröhlich model, have also been realized and their dynamics experimentally investigated Catani12 (); Widera12 (); Oberthaler13 (). However, so far, there have been no studies exploiting Feshbach resonances to increase the strength of interspecies interactions, nor measuring basic polaron properties such as their binding energy, lifetime and effective mass. On the theoretical side, the selflocalization of Bose polarons was investigated using meanfield approaches Cucchietti06 (); Kalas06 (); Jaksch08 (); Santamore11 (); Timmermans13 () as well as Feynman’s variational method applied to the effective Hamiltonian describing the impurity Tempere09 (); Novikov10 (). Starting from the Fröhlich Hamiltonian other studies have focused on the calculation of the radio frequency response of the polaron Demler14 (), and of its binding energy and effective mass using renormalization group Demler141 () and diagrammatic Monte Carlo Devreese14 () methods. A more microscopic approach based on the Tmatrix approximation was used in Ref. Rath13 () where various quasiparticle properties are calculated for both attractive and repulsive Bose polarons close to a Feshbach resonance. Similar results are also obtained in Ref. DasSarma14 () by means of a variational ansatz for the wave function of the bathimpurity system. Finally, threebody correlations were explicitly included in the theoretical treatment both at the level of perturbation theory Bruun1 () as well as within a variational approach Bruun2 (), giving rise to a significant lowering of the binding energy of attractive polarons.
In this paper we address the problem of Bose polarons using a fully microscopic, non perturbative approach, consisting in the quantum Monte Carlo (QMC) method. This numerical technique can provide exact results for the groundstate energy and the effective mass of the impurity as a function of the parameters of the Hamiltonian describing the interspecies and intraspecies interaction potentials and the density of the bosonic bath. We model these interactions using a hardsphere potential for the interboson repulsion and both a purely repulsive hardsphere and an attractive squarewell potential for the impurityboson interaction. In particular, the latter model allows one to investigate situations where the impurityboson wave scattering length is either positive or negative giving rise to the groundstate attractive and excitedstate repulsive branches of the polaron. Our analysis is limited to the case where the mass of the impurity is equal to the one of the Bose particles in the medium, but generalizations to include different mass ratios can be easily implemented within the same method.
We investigate the properties of the Bose polaron both along the attractive and the repulsive branch. We find that for small values of the ratio of the impurityboson to the bosonboson scattering length our results for the binding energy and the effective mass are in good agreement with secondorder perturbation theory based on a Frölichlike Hamiltonian describing the coupling between the impurity and the bath. At the unitary point of resonant impurityboson scattering () the binding energy is found to scale with , where is the density of the bath and is the common mass of the particles and the impurity. We notice that this behavior is similar to the Fermi polaron case where the binding energy is proportional to the Fermi energy of the bath Massignan14 (). The effective mass ratio ranges from values close to one in the weakcoupling limit up to values which remain smaller than two close to the resonant point. We find no evidence of the selflocalization of the polaron, which in studies based on the Fröhlich model is signaled by an abrupt increase of the effective mass as the coupling strength exceeds a critical value. We believe that this wrong prediction has to do with the inadequacy of the effective Fröhlich Hamiltonian in the description of the pairing mechanism which takes place close to the resonance where the impurity and one boson from the bath can form a bound state.
We analyze the structural properties of the bosonic bath by calculating the contact parameter which characterizes the shortrange behavior of the impurityboson pair correlation function. The knowledge of how particles in the bath are distributed around the impurity enables us to evaluate the distortion of the density profile produced by the impurity. Within the attractive squarewell model, we find a pronounced peak in the density close to the impurity both on the attractive and on the repulsive branch. This peak is a result of the pairing induced by the impurityboson potential. It is a shortrange feature that can not be accounted for by the Fröhlich Hamiltonian which can only describe longrange distortions of the density profile.
An important point to analyze is related to the existence of fewbody bound states in vacuum, such as threebody Efimovlike states and deeper bound states with more than three particles. At the resonant point we calculate the energy of the deepest bound state with three and more particles (i.e. the impurity plus two or more bosons) finding evidence that such state exists only up to six particles (i.e. the impurity plus five bosons). Remarkably, the energy of these selfbound states is in absolute value much smaller than the polaron binding energy which involves the contribution from a large number of particles in the bath. One should notice that the values we obtained for the groundstate energy of the cluster states as well as the size of the largest cluster greatly depend on the details of the hardsphere bosonboson potential used in the simulations. However, we believe that the results for the polaron binding energy at unitarity are universal and only depend on the gas parameter of the bath and the mass ratio between the impurity and the bosons.
The structure of the paper is as follows. In Sec. II we first address the singlepolaron problem by introducing the model Hamiltonian (subsection IIA), and by reviewing the perturbation treatment leading to the Fröhlichtype Hamiltonian. Here, we also derive the results for the polaron binding energy and effective mass valid in the weakcoupling limit (subsection IIB). Finally, in subsection IIC, we briefly review the DMC method and we discuss the different trial wave functions used to describe the attractive and repulsive polaron branch. The results along the two branches concerning binding energy, effective mass, density profiles and contact parameter are presented in subsection IID. Furthermore, subsection IIE contains a discussion of these results specific of the resonant point for the impurityboson scattering. In Sec. III we report on calculations of the binding energy of fewbody states in vacuum at the unitary point and on the side of the resonance where a twobody bound state exists. In Sec. IV, we generalize the problem to many impurities obeying Bose statistics and we use DMC simulations to validate the perturbative equation of state in the limit of small concentrations bearing some consequences for the phase diagram of binary mixtures. Conclusions are finally drawn in Sec. V.
Ii II. Single impurity
ii.1 A. Model Hamiltonian
We consider a system of one impurity immersed in a dilute gas of Bose particles at described by the following Hamiltonian
(1)  
Here, the first two terms represent the kinetic and the interaction energy of the bosonic bath consisting of particles of mass and interacting through the twobody potential , which depends on the distance between a pair of bosons. Furthermore, is the kinetic energy of the impurity with mass denoted by the coordinate vector and is the bosonimpurity potential depending on the distance between the impurity and the th particle of the bath. The interboson potential is modeled by the hardsphere (HS) interaction
(2) 
where the diameter coincides with the wave scattering length. The impurityboson interaction, instead, is modeled by either a purely repulsive hardsphere potential
(3) 
or an attractive squarewell (SW) potential
(4) 
The latter is characterized by a range and a depth () chosen such as to yield the value of the scattering length. This is determined by the transcendental equation
(5) 
where in terms of the reduced mass . In the case of the repulsive potential the wave scattering length is always positive, whereas for the squarewell potential in Eq. (5) the value of can be either positive or negative depending on . In particular, we consider values in the range , corresponding to either no bound state () or one bound state () in the twobody sector. In this latter case the molecular binding energy is obtained from the equation
(6) 
where . We also notice that the value corresponds to the unitary limit of the impurityboson interaction where the scattering length diverges and the binding energy vanishes. In the present study we consider values of the range that are small compared to the interboson scattering length with ranging between 5 to . We expect that for such a shortrange potential the value of its scattering length is the only relevant parameter for all polaron properties.
We restrict the analysis of the Hamiltonian (1) to the case where the impurity and the bosons in the bath have the same mass: . The strength of the interboson interactions is determined by the gas parameter involving the bosonic density , whereas the intensity of the impurityboson coupling is given in terms of the ratio of the two scattering lengths.
ii.2 B. Perturbation theory
The problem of a single mobile impurity in a Bose gas can be thoroughly investigated using perturbation theory at least in the weakcoupling regime. The approach is based on the treatment of the bath within the Bogoliubov approximation of a dilute Bose gas described by the Hamiltonian
(7) 
Here, is the groundstate energy of the bosonic particles
(8) 
where is the interboson coupling constant. The operators , are the annihilation/creation operators of quasiparticles related to the bosonic particle operators , through the standard transformations
(9) 
with coefficients and . The elementary excitation energies are given by the Bogoliubov spectrum
(10) 
where is the dispersion of free particles and is the density of condensed particles.
At the meanfield level the interaction energy between the impurity, located in position , and the bath is described by the following expression
(11) 
involving the density of bosons and the interspecies coupling constant proportional to the bosonimpurity swave scattering length . In momentum space the above interaction Hamiltonian can be recast in the form
where the Bogoliubov approximation of the density fluctuations in terms of quasiparticle operators (9) has been used. By applying perturbation theory to the Hamiltonian , as it was done for example in Ref. Viverit02 () for a BoseFermi mixture, one finds the following result for the groundstate energy of the system of bosons plus one impurity
(13) 
valid to order in the bosonimpurity coupling strength Saam69 (); Tempere09 (). To the same order one can also calculate the effective mass of the impurity obtaining the result note1 (); Cucchietti06 ()
(14) 
Apart from the trivial firstorder contribution to the groundstate energy, the secondorder corrections in Eqs. (13)(14) scale in terms of the same dimensionless parameter , which should be much smaller than unity to ensure the validity of the perturbation approach. We also notice that the perturbative corrections to the energy and the effective mass in Eqs. (13)(14) diverge when the interboson scattering length tends to zero. This feature indicates the instability of the ideal Bose gas towards clusterization around the impurity and points out the crucial role played by the repulsive interaction between the bosons.
The Hamiltonian (II.2) has the general form of the Fröhlich polaron Hamiltonian describing an impurity coupled to a bath of non interacting bosonic quasiparticles Froehlich54 (). This analogy was first exploited in Ref. Cucchietti06 () where the intriguing problem of the “selflocalization” of the polaron was addressed in a fashion similar to the LandauPekar description of electrons in ionic crystals Landau46 (). In Ref. Tempere09 () the JensenFeynman variational scheme is applied to both the strong and the weakcoupling regime of the effective Fröhlich Hamiltonian. Both these studies predict that for the impurity “selflocalizes” inside the potential well produced by its distortion of the bosonic density. In terms of the effective mass of the impurity the selflocalization phenomenon corresponds to a large enhancement of the mass ratio, , which arises from the cloud of bosonic quasiparticles dressing the impurity. We would like to stress, however, that the Fröhlich Hamiltonian , where the bath and interaction term are given respectively by Eq. (7) and (II.2), is an effective lowenergy reduction of the microscopic Hamiltonian (1). Whether it captures the relevant physics of an impurity immersed in a Bose condensate when the coupling to the bath is strong is questionable and should be analyzed with care. This problem will be investigated in the remaining part of the article using quantum Monte Carlo methods that are particularly suitable to treat strongly correlated systems in a non perturbative manner.
ii.3 C. Quantum Monte Carlo method
We use the diffusion Monte Carlo (DMC) method which aims to solve the manybody Schrödinger equation in imaginary time for the distribution function , where is the wave function of the system of bosons plus one impurity and is a trial wave function of the particle coordinates used for importance sampling. The timedependent Schrödinger equation can be written as
(15)  
in terms of the so called local energy and quantum drift force . In the above equation plays the role of a diffusion constant and is a reference energy. The formal solution of the equation is given by
(16) 
where one introduces the Green’s function describing the time evolution governed by the Langevin operator . If the shorttime dependence of is known for sufficiently small , the asymptotic solution for large times, , can be obtained by iterating Eq. (16) for a large number of time steps . On general grounds, the initial time distribution is expanded in terms of the eigenfunctions of the system corresponding to the eigenenergies . For a system of bosons, provided the coefficient does not vanish, the solution at large times of Eq. (15) is given by the expression
(17) 
and is proportional to the nodeless ground state with energy . The value of the groundstate energy is determined from the condition of keeping the probability distribution stationary at large times or, more conveniently, from the average of the local energy
(18) 
Apart from statistical errors, the DMC method allows one to calculate the exact groundstate energy of a system of Bose particles. Importantly, the energy estimate obtained using the DMC technique with importance sampling is to a large extent independent of the detailed shape of the trial wave function as long as is positive definite. If the function changes sign in some regions of the configuration space, the DMC algorithm guided by this trial function yields, instead of the ground state, the lowestenergy eigenstate compatible with the nodal constraint fixed by .
The general form of the trial wave function used in the present study is given by
(19) 
where the functions and describe, respectively, interboson and impurityboson twobody correlations. The functional form of the interboson term is constructed from the twobody scattering solution of the hardsphere potential in Eq. (2)
(20) 
where the value of the wave vector is chosen such that the first derivative of the function vanishes at half of the size of the cubic simulation box: . This condition ensures that the Jastrow factor (20) is compatible with the periodic boundary conditions used in the simulation.
For the impurityboson correlation function we use instead different forms depending on the type of potential , hard sphere or square well, and on the polaron branch, repulsive or attractive.

Squarewell potential: We use different functional forms of for the repulsive and attractive branch of the polaron.

Repulsive branch (): is constructed from the zeroenergy scattering solution of the potential (4), orthogonal to the bound state existing at for two particles:
(22) Here is a matching point and for the function goes to a constant reached at where . The coefficients , and ensure the continuity of and of its first derivative at the points and .

Attractive branch (): is constructed in the same way as for the repulsive branch of Eq. (22), with the only difference that the scattering length is negative.

Attractive branch (): We use the solution of the twobody bound state with energy given by Eq. (6)
(23) where . The coefficients , and again ensure the continuity of and of its first derivative at the potential range and at the matching point .
In all the above three cases, the values of the matching point and of the parameter are optimized by minimizing the variational energy. We notice that the function is positive definite along the attractive branch, whereas it changes sign at the value on the repulsive branch. For positive values of the nodal surface in the manybody trial wave function which originates from the choice (22) of the Jastrow correlation term allows one to discriminate between the groundstate attractive branch and the excitedstate repulsive branch.
Furthermore, we notice that the unitary limit, corresponding to , is reached following the attractive branch. This limit corresponds to a Jastrow term in the range and is obtained by approaching the resonance both from and from .
ii.4 D. Attractive and repulsive polaron branch
Simulations are carried out using periodic boundary conditions and the number of bosons in the bath is typically . Calculations with different numbers of particles up to are also performed in order to check that finitesize effects are below statistical uncertainty.

Binding energy
We determine the polaron binding energy by calculating the energy difference
(24) 
where is the groundstate energy of the system of bosons alone and is the energy of the system of bosons plus the impurity in the same volume .
Two different branches are obtained depending on the impurityboson interaction potential and on the choice of the Jastrow term in the trial function (19). The attractive branch, which corresponds to being the ground state of the composite system, is simulated using the squarewell potential (4) and the positive definite function described in Sec. II C. Along the repulsive branch is still the ground state of the hardsphere potential (3), but it corresponds to an excited state of the potential (4) which we calculate by imposing the nodal constraint given by Eq. (22) on the trial function.
The results for both branches are shown in Fig. 1. The gas parameter of the bosonic bath is here corresponding to a dilute gas whose groundstate energy is found to be very close to the result (8) of second order perturbation theory. The results reported in Fig. 1 are obtained both with the hardsphere potential (3) and with the squarewell potential (4) where we used two different values of the ratio of the bosonboson scattering length to the potential range. The figure clearly indicates that the energies scale with the ratio and that the details of the impurityboson potential are irrelevant. For the repulsive branch we find a remarkably good agreement with the perturbation result in Eq. (13) up to values of . For larger values of the impurityboson scattering length calculations using the squarewell potential get increasingly difficult because of large fluctuations arising from the nodal constraint imposed on the Jastrow correlation term which is no longer adequate to define the excited state of the polaron. It is worth stressing at this point that the results for the repulsive branch strongly depend on the choice made for the nodal constraint which provides the correct description of the repulsive polaron only in the limit . The attractive branch can instead be followed, starting from small negative values of , down to the unitary point () and when approaching this point the polaron binding energy shows large deviations from the perturbation expansion (13). We find that at the unitary point , resulting in a binding energy of the impurity much larger than the chemical potential of bosons in the bath. The results corresponding to positive values of along the attractive branch are reported in the inset of Fig. 1). We find that lies always significantly below the twobody binding energy . One expects that, by increasing on the positive side of the resonance, the polaron binding energy eventually approaches the energy of the deepest cluster state (see Sec. III).
In Fig. 2 we show a more detailed comparison of the binding energy of the two polaron branches with the perturbation expansion (13). This is carried out by subtracting from the values of the meanfield contribution , i.e. the first term in bracket in Eq. (13), and by comparing with the second order contribution of the perturbation expansion. We notice that remains positive in the limit , in contrast with the predictions of Ref. Cucchietti06 () and Tempere09 () where this quantity should turn negative for as a consequence of the selflocalzation of the polaron.

Effective mass
The effective mass of a distinguishable particle can be determined in a DMC simulation by calculating its diffusion constant in imaginary time Boninsegni95 (); Boronat99 (). The main assumption is that the energy of the system with the impurity having momentum can be written in the form
(25) 
in terms of the impurity binding energy and effective mass . The ratio of the bare to the effective mass of the particle is then given by
(26) 
where is the diffusion constant of a free particle and is the mean square displacement of the impurity in imaginary time. One can determine the value of from the large time slope of as a function of the imaginary time . The results are shown in Fig. 3 for the attractive and repulsive branch of the polaron. We notice that far away from the resonant point the increase of the effective mass agrees with the prediction of perturbation theory. On approaching the resonance, remains finite reaching values . This result is again in contrast with the selflocalization picture of Refs. Cucchietti06 (); Tempere09 () which predicted a large increase of the effective mass with increasing coupling strength. Along the attractive branch we calculated only up to the unitary point, we expect that following this branch on the positive side of the resonance the value of should continue to increase.

Density profiles
Another important output of our QMC simulations, useful to understand the changes induced in the bosonic bath by the impurity, is the pair correlation function giving the probability of finding a bosonic particle at a distance from the impurity. At large distances , whereas its shortrange behavior is determined by the impurityboson potential . The density profile of the particles of the bath surrounding the impurity can be calculated using the following integral of the pair correlation function
(27) 
which approaches the bulk value far away from the impurity. As a technical remark, we compute the function by carrying out both a variational and a diffusion Monte Carlo calculation, which provide the estimates and respectively, and by using the extrapolation formula Kolorenc11 ().
The ratio of the local to the bulk density is shown in Figs. 46 for three values of , both on the attractive and the repulsive branch of the polaron. The calculations are carried out at the bath density , using the squarewell potential (4) with , and distances are reported in units of the healing length . In all cases exhibits a pronounced peak at the position of the impurity caused by the attractive potential well. If (attractive branch) the local density of particles decreases monotonously, reaching the bulk value when . Instead, if (repulsive branch) the density goes through a minimum before reaching the bulk value. The position of the minimum lies in the region and decreases with increasing . In this case, the density depletion occurring at a large distance from the impurity arises from the effective repulsive interaction associated with the positive value of .
The average number of particles of the bath surrounding the impurity is obtained from the integral and is shown in the inset of Figs. 46 as a function of the distance . We notice that the number starts to grow faster with increasing for positive values of , consistently with the larger peak of at very short distances. In particular, for and 20 (see Figs. 45), rapidly reaches the value of one particle already at distances .
We also notice that any perturbative approach based on the Fröhlichtype Hamiltonian (II.2) can only be meaningfully applied if the density perturbation induced by the impurity, , satisfies the condition Demler141 (). These approaches are therefore limited to values of such that and can never describe correctly the structural properties of the bath at short distances from the impurity.

Contact parameter
The impurityboson contact parameter can be determined from the behavior of the pair correlation function in the range of distances , but still much larger than the typical radius of the impurityboson interaction. We define the dimensionless contact parameter as
(28) 
where the limit should be intended in the sense specified above. The results for the contact parameter, obtained at using the SW potential with , are shown in Fig. 7 for both the attractive and the repulsive branch. For small values of the impurityboson pair correlation function is well approximated by the simple expression , determined solely by twobody physics, yielding the estimate for the contact parameter. The derivative of the polaron binding energy with respect to the inverse scattering length should also be related to Werner12 (). From the behavior in the regime one finds
(29) 
This result is shown in Fig. 7 together with the contact extracted from the pair correlation function. Good agreement is found along the attractive branch, whereas the two estimates of on the repulsive branch are compatible only in the weakcoupling limit. The disagreement between the contact parameter obtained from the equation of state and from the pair correlation function indicates that our choice of the trial wave function does not provide a fully satisfactory description of the repulsive polaron in the region where becomes very large.
ii.5 E. Resonant interaction
In this section we focus on the properties of the Bose polaron when the interaction between the impurity and the bath is resonant, i.e. . In Fig. 8 we show the binding energy of the polaron calculated at resonance as a function of the gas parameter of the bath. The results show that scales with the energy and that, once expressed in these units, it depends weakly on the gas parameter over many orders of magnitude. As decreases the value of also decreases, reaching at the very small density . We can not establish weather the binding energy continues to decrease for even smaller densities, signaling the instability of the non interacting gas in the presence of an impurity with attractive interaction, or it reaches a constant value in agreement with the findings of the fieldtheoretical calculation in Ref. Rath13 (). Remarkably, the binding energy of a Fermi polaron resonantly interacting with the bath is given by , where is here the density of the Fermi sea Lobo06 (), and differs approximately by a factor of two compared to the results in Fig. 8 for the smallest values of .
The effective mass as a function of the gas parameter is shown in Fig. 9. Also in this case we find a small variation of following a change of over orders of magnitude. The largest effective mass, , is achieved at the smallest densities of the bath.
In Fig. 10 we show the density profile of the bath surrounding the impurity obtained using the pair correlation function and Eq. (27). The behavior is qualitatively similar to the one reported in Fig. 6 and corresponding to along the attractive polaron branch. By decreasing the value of the bath gas parameter we find that the density peak around the impurity sharpens and the size of the deformation in units of the healing length shrinks. Finally, in the inset of Fig. 10, we show the value of the contact parameter , determined from the shortrange behavior of the pair correlation function, for different values of . Also for this quantity we observe a weak dependence on the value of the bath gas parameter.
Iii III. Fewbody physics
In this Section we consider the problem of the existence of bound states in vacuum consisting of the impurity and a number of bosons. Of course, such bound states can only occur in the case of the squarewell model for the impurityboson potential. This potential supports a twobody molecular state having energy , given by Eq. (6), for all positive values of the interspecies scattering length .
The search for the ground state of clusters with particles is carried out using the DMC method based on the following trial wave function
(30) 
The above wave function differs from the one of Eq. (19), used in simulations of homogeneous configurations, by the exponential term which depends on the hyperradius of the cluster
(31) 
where is the coordinate of the center of mass. The Jastrow correlation terms in Eq. (30) are similar to the ones of Eqs. (20) and (23), respectively for the bosonboson and the impurityboson function. Since periodic boundary conditions are absent here, the length scale in Eq. (20) is replaced by the large distance . Moreover, the boundary condition on the derivative of is relaxed with the choice , holding for , with the constants and determined in the same way as in Eq. (23). Free parameters that are optimized using a variational procedure are the matching point and the coefficients and . In particular, the latter fixes the size of the cluster in terms of its hyperradius.
Calculations are performed in the reference frame where , in order to eliminate the contribution from the center of mass motion. Furthermore, we consider only the resonant point where and , and the point on the positive side of the resonance where . In Fig. 11 we show the results for the groundstate energy of the cluster with particles as a function of the number of bosons. At unitarity the twobody binding energy, corresponding to in Fig. 11, is identically zero, whereas the threebody Efimov state () is found to feature an extremely shallow groundstate energy: . This result is consistent with the prediction for the lowest Efimov state in terms of the threebody length and the Efimov parameter Petrov10 (). In the case of equal masses for the impurity and the bosons, the value of is very small, , resulting in which is of the same order as our estimate if . For increasing the cluster groundstate energy decreases markedly up to , while clusters with are undoubtedly unbound. These findings are compatible with the results at (see inset of Fig. 11), where the binding energy appears not to decrease further already for .
It is important to stress that the largest size of bound clusters, as well as the precise value of their groundstate energies, depend on the details of the interboson and impurityboson interactions. However, we believe that the qualitative behavior emerging from our simulations should hold for any shortrange interaction with scattering length and , respectively. In particular, we notice that the polaron binding energy shown in Fig. 1 is more than a factor larger than the deepest cluster state at unitarity and remains larger also at . A possible reason for the irrelevance of cluster states at unitarity is the feature of equal masses for the impurity and the bosons, which makes Efimov states extremely shallow.
Iv IV. Many impurities
Let us now analyze the case of a small concentration of impurities immersed in a BEC at . For this problem the statistics of the impurities is important and in the present study we consider only impurities which obey Bose statistics. A binary BoseFermi mixture with a small concentration of bosons in a Fermi sea and featuring resonant BoseFermi interactions has been investigated using QMC methods in Ref. Bertaina13 ().
A collection of impurities immersed in a gas of particles is described by the Hamiltonian
The impurityimpurity potential is modeled by the same hardsphere interaction, including the same scattering length , which characterizes the coupling between the bosons of the bath: . We also assume that the masses of the two types of particles are the same (), and the impurityboson interaction is as described in Sec. IIA.
The perturbation treatment of a binary mixture of Bose condensates at has been carried out in Ref. Balabanyan86 () using an extension of the standard Bogoliubov approach. The result for the groundstate energy in the low concentration limit, , is obtained as follows
(33)  
up to quadratic contributions in the ratio of scattering lengths and in the impurity concentration. Here is the energy (8) of the bath without impurities. Furthermore, one should notice that the term linear in the concentration coincides with the single polaron energy of Eq. (13), whereas the term proportional to describes the repulsive interaction between polarons.
In QMC simulations we calculate the groundstate energy of the mixture of bosons plus impurities making use of the following trial wave function
(34) 
where the same Jastrow factor of Eq. (21) accounts for the repulsive correlations between the particles of the bath and between the impurities while the impurityboson term is as described in Sec. IIC. We notice that the above wave function is symmetric under the exchange of the impurity coordinates fulfilling Bose statistics.
We perform QMC calculations using particles in the bath and a varying number of impurities in a cubic box with periodic boundary conditions, aiming to simulate a homogeneous binary mixture characterized by a small concentration, , of the minority component. In particular, we calculate the shift between the groundstate energy of the mixture and of the bath without impurities:
(35) 
The results for different values of the ratio along the repulsive branch are reported in Fig. 12 for the bath density . The energy shift is in remarkable agreement with the perturbation result (33) for all values of up to . Only at the largest values of and of the concentration significant deviations from Eq. (33) are visible. In particular, the interaction between polarons appears to be overestimated by the term proportional to in Eq. (33). In Fig. 13 we show the results of as a function of the concentration for on both the repulsive and the attractive branch. Also in this case, the comparison with Eq. (33) shows a very good agreement.
The above results validate the expression (33) for the equation of state of a binary BoseBose mixture in the regime of parameters: and . Such a validation is important in order to establish the stability conditions and the phase diagram of the mixture in a quantitatively reliable way. From the analysis of the compressibility matrix , where and , one finds that the homogeneous mixture is stable if the ratio of scattering length satisfies the condition
(36) 
This relation holds in the limit and, compared to the meanfield result given by Balabanyan86 (), includes the leading order correction in the small parameter . We notice that the results reported in Figs. 12 and 13 lie outside the stability range (36) of the homogeneous binary mixture. The spinodal instability arising from a vanishing compressibility is associated to systems approaching the thermodynamic limit and is usually prevented in simulations of finitesize systems.
Another possible state of the binary mixture corresponds to a complete phase separation between the bosons and the impurities. In this case the groundstate energy can be written as
(37)  
in terms of the densities and of the two species and their relative volumes fulfilling the condition . The stability of this state requires that and that the energy cost to add one impurity to each of the two uniform phases is positive: . From Eqs. (8) and (13), giving respectively the energy of the bath without impurities and the excess energy of a single polaron, we get that the phase separated state is stable if
(38)  
(39) 
One can easily show that, if , the energy (37) of the phase separated state lies below the energy (33) of the homogeneous binary mixture for any value . On the contrary, no gaslike phase appears to be stable outside the regions of Eqs. (38) and (36) when .
V V. Conclusions
We investigated the properties of an impurity immersed in a BoseEinstein condensate at . This Bose polaron study has been carried out using a fully microscopic approach and QMC simulation methods in analogy with previous investigations of the more thoroughly expounded, both theoretically and experimentally, Fermi polaron problem. The main results concern the binding energy of the impurity and its effective mass along the attractive and repulsive polaron branch, explored by changing the scattering length of the impurityboson interaction potential. These results are expected to be universal, for a given bosonboson and impurityboson scattering length and for equal masses of the two components. At the resonant point of the impurityboson interaction the polaron binding energy scales with the equivalent of the Fermi energy of the bath in analogy with the behavior found for Fermi polarons.
The measurement of the polaron binding energy should be accessible in experiments using radiofrequency spectroscopy while its effective mass affects the dispersion of collective excitations in highly imbalanced twocomponent mixtures. We hope that such experiments, which were so successful in the investigation of the properties of the Fermi polaron, will be carried out also for its Bose counterpart providing a deeper knowledge of this clean and simply stated, but highly non trivial manybody problem.
Acknowledgments
We gratefully acknowledge useful discussions with J. Tempere, W. Zwerger, R. Grimm, M. Cetina and D. S. Petrov. This work has been supported by ERC through the QGBE grant.
References
 (1) L. D. Landau and S. I. Pekar, Zh. Eksp. Teor. Fiz. 16, 341 (1946).
 (2) H. Fröhlich, Adv. in Phys. 3, 325 (1954).
 (3) R. P. Feynman, Phys. Rev. 97, 660 (1955).
 (4) See, e.g., G. D. Mahan, ManyParticle Physics (Plenum, New York, 1990), Chap. 6, 2nd ed.
 (5) N. V. Prokof’ev and B. V. Svistunov, Phys. Rev. Lett. 81, 2514 (1998).
 (6) A. S. Mishchenko, N. V. Prokof’ev, A. Sakamoto, and B. V. Svistunov, Phys. Rev. B 62, 6317 (2000).
 (7) For a recent review see, e.g., J. T. Devreese and A. S. Alexandrov, Rep. Prog. Phys. 72, 066501 (2009).
 (8) F. M. Peeters and J. T. Devreese, Phys. Rev. B 32, 3515 (1985).
 (9) F. M. Cucchietti and E. Timmermans, Phys. Rev. Lett. 96, 210401 (2006).
 (10) R. M. Kalas and D. Blume, Phys. Rev. A 73, 043608 (2006).
 (11) M. Bruderer, W. Bao, and D. Jaksch, Eur. Phys. Lett. 82, 30004 (2008).
 (12) J.Tempere, W. Casteels, M. K.Oberthaler, S.Knoop, E. Timmermans, and J. T. Devreese, Phys. Rev. B 80, 184504 (2009).
 (13) P. Massignan, M. Zaccanti, and G. E. Bruun, Rep. Prog. Phys. 77, 034401 (2014).
 (14) A. Schirotzek, C.H. Wu, A. Sommer, and M. W. Zwierlein, Phys. Rev. Lett. 102, 230402 (2009).
 (15) S. Nascimbène, N. Navon, K. J. Jiang, L. Tarruell, M. Teichmann, J. McKeever, F. Chevy, and C. Salomon, Phys. Rev. Lett. 103, 170402 (2009).
 (16) C. Kohstall, M. Zaccanti, M. Jag, A. Trenkwalder, P. Massignan, G. M. Bruun, F. Schreck, and R. Grimm, Nature 485, 615 (2012).
 (17) M. Koschorreck, D. Pertot, E. Vogt, B. Fröhlich, M. Feld, and M. Köhl, Nature 485, 619 (2012).
 (18) J. Catani, G. Lamporesi, D. Naik, M. Gring, M. Inguscio, F. Minardi, A. Kantian, and T. Giamarchi, Phys. Rev. A 85, 023623 (2012).
 (19) N. Spethmann, F. Kindermann, S. John, C. Weber, D. Meschede, and A. Widera, Phys. Rev. Lett. 109, 235301 (2012).
 (20) R. Scelle, T. Rentrop, A. Trautmann, T. Schuster, and M. K. Oberthaler, Phys. Rev. Lett. 111, 070401 (2013).
 (21) D. H. Santamore and E. Timmermans, New J. Phys. 13, 103029 (2011).
 (22) A. A. Bilinova, M. G. Boshier, and E. Timmermans, Phys. Rev. A 88, 053610 (2013).
 (23) A. Novikov and M. Ovchinnikov, J. Phys. B: At. Mol. Opt. Phys. 43, 105301 (2010).
 (24) A. Shashi, F. Grusdt, D. A. Abanin, and E. Demler, Phys. Rev. A 89, 053617 (2014).
 (25) F. Grusdt, Y. E. Shchadilova, A. N. Rubtsov, and E. Demler, preprint arXiv:1410.2203.
 (26) J. Vlietinck, W. Casteels, K. Van Houcke, J. Tempere, J. Ryckebusch, and J. T. Devreese, preprint arXiv:1406.6506.
 (27) S. P. Rath and R. Schmidt, Phys. Rev. A 88, 053632 (2013).
 (28) Weiran Li and S. Das Sarma, Phys. Rev. A 90, 013618 (2014).
 (29) R. S. Christensen, J. Levinsen, and G. M. Bruun, preprint arXiv:1503.06979.
 (30) J. Levinsen, M. M. Parish, and G. M. Bruun, preprint arXiv:1505.04530.
 (31) L. Viverit and S. Giorgini, Phys. Rev. A 66, 063604 (2002).
 (32) W. F. Saam, Ann. Phys. 53, 239 (1969).
 (33) This result was first derived in Ref. Saam69 () for equal scattering lengths and as a function of the mass ratio .
 (34) M. Boninsegni and D. M. Ceperley, Phys. Rev. Lett. 74, 2288 (1995).
 (35) J. Boronat and J. Casulleras Phys. Rev. B 59, 8844 (1999).
 (36) For more details see, e. g., J. Kolorenc̆ and L. Mitas, Rep. Prog. Phys. 74, 026502 (2011).
 (37) F. Werner and Y. Castin, Phys. Rev. A 86, 053633 (2012).
 (38) C. Lobo, A. Recati, S. Giorgini, and S. Stringari, Phys. Rev. Lett. 97, 200403 (2006).
 (39) K. Helfrich, H.W. Hammer, and D. S. Petrov, Phys. Rev. A 81, 042715 (2010).
 (40) G. Bertaina, E. Fratini, S. Giorgini, and P. Pieri, Phys. Rev. Lett. 110, 115303 (2013).
 (41) G. O. Balabanyan, Teoret. Mat. Fiz. 66, 121 (1986).