Phase diagram for strong-coupling Bose polarons

Important properties of complex quantum many-body systems and their phase diagrams can often already be inferred from the impurity limit. The Bose polaron problem describing an impurity atom immersed in a Bose-Einstein condensate is a paradigmatic example. One of the most interesting features of this model is the competition between the emergent impurity-mediated attraction between the bosons and their intrinsic repulsion. The arising higher-order correlations make the physics rich and interesting, but also complex to describe theoretically. To tackle this challenge, we develop a quantum chemistry-inspired computational technique and compare two state-of-the-art variational methods that fully include both the boson-impurity and boson-boson interactions on a non-perturbative level. For a sweep of the boson-impurity interaction strength, we find two regimes of qualitatively different behaviour. If the impurity-mediated interactions overcome the repulsion between the bosons, the polaron becomes unstable due to the formation of large bound clusters. If instead the interboson interactions dominate, the impurity will experience a crossover from a polaron into a small molecule. We achieve a unified understanding incorporating both of these regimes and the transition between them. We show that both the instability and crossover regime can be studied in realistic cold-atom experiments. Moreover, we develop a simple analytical model that allows us to interpret these phenomena in the typical Landau framework of first-order phase transitions that turn second-order at a critical endpoint, revealing a deep connection of the Bose polaron model to both few- and many-body physics.


Introduction
Understanding the properties of complex quantum many-body systems is one of the major challenges in modern physics [1,2].Gases of ultracold atoms provide a unique playground to study such systems in a controlled environment [3][4][5][6][7].From a theoretical perspective, a special feature of these systems is the practically complete understanding of the microscopic interactions between the atoms [8].Moreover, the universality of ultracold scattering enables the translation of these detailed models into simple effective Hamiltonians, which can still be used to describe experiments on a quantitative level.This has led to great synergy between theory and experiments and a significantly improved understanding of paradigmatic systems such as the Fermi-Hubbard model [9][10][11][12][13][14][15][16], the unitary Bose-Einstein condensate (BEC) [17][18][19][20][21][22][23][24], and the crossover from a BEC to a BCS (Bardeen-Cooper-Schrieffer) state of paired fermions [25][26][27][28][29][30].Nevertheless, even though the theoretical models and Hamiltonians seem simple at first sight, important aspects of these systems are still not understood.
Particularly interesting systems are ultracold mixtures of different species of bosons and/or fermions.For example, boson-mediated interactions between fermions are a classic and important mechanism of superconductivity [31].With cold atoms, also the opposite scenario of fermion-mediated interactions in a BEC has been studied [32,33].Mixtures of bosons on the other hand, can form quantum droplets [34][35][36][37][38][39].Other than being of fundamental interest, ultracold mixtures are also a common starting point for creating ultracold molecules [40][41][42][43][44]. or crossover (dashed) as a function of the interaction strength.In the instability regime, the polaron state becomes unstable at the scattering length indicated by the dot.At this point the polaron decays into a cluster.The energy of the cluster before the instability is drawn with the dash-dotted line (only shown in the regime where it is lower than the polaron energy).The bold and dash-dotted lines are computed using a Gaussian-state Ansatz, and the dashed line using a double-excitation Ansatz (see Sec.A good understanding of ultracold mixtures can therefore also be harnessed to create ultracold molecules more efficiently [45].
The rich properties and the variety of phenomena that make these ultracold mixtures so interesting have the immediate consequence that the development of a unified theoretical picture or a global phase diagram is extremely challenging [46].A successful approach to achieve essential insights has been to start from the quantum impurity limit, of a single particle of one species immersed in a sea of the other.Often, key signatures of the many-body properties of the full mixture can already be recognized in this well-defined limit.
A good example is the problem of the Fermi polaron [47][48][49][50][51][52], an impurity immersed in a Fermi sea.As a function of interaction strength, this model shows a transition from a polaronic to a molecular ground-state, which is a precursor of the BCS-BEC crossover.Furthermore, this model shows an interesting connection to the elusive Fulde-Ferell-Larkin-Ovchinnikov phase of superconductivity [53][54][55][56].
The Bose polaron problem [57][58][59][60][61][62] of an impurity immersed in a BEC has proven to be more challenging.We focus on the case of a charge-neutral impurity, but also the case of an ionic impurity has drawn recent theoretical interest [63][64][65].Despite much effort, a unified understanding is lacking of what happens to the impurity as a function of interaction strength.The complications which arise here are the large number of excitations around the impurity, the importance of higher-order correlations, and the interactions between bosons from the BEC.In three dimensions, 1 no theoretical approach has so far captured all of these aspects simultaneously.
In Fig. 1a) we illustrate how the character of the attractive Bose polaron is predicted to change as the boson-impurity interactions are swept across a Feshbach resonance.Specifically, we show the energy of the polaron as a function of the inverse scattering length.State-of-theart theories envision two different scenarios.In the first scenario, the polaron experiences a smooth crossover into a small molecule, such as a dimer, trimer, or tetramer [67][68][69].In the second scenario, an instability is predicted, in the form of a collapse of the BEC triggered by attractive impurity-mediated interactions [70,71].The result is the decay of the polaron into a multi-particle bound cluster.These scenarios are qualitatively different, and so far they have not been captured together in a single theoretical framework.
In this work we achieve a unified picture of the properties of the attractive Bose polaron as a function of the BEC density, impurity mass, the boson-impurity and boson-boson interaction strengths.We show in Fig. 1b), that, depending on the system, both of the scenarios in Fig. 1a) can occur.The largest part of the parameter space corresponds to the crossover regime.The crossover occurs if the impurity is relatively heavy or similar in mass to the bosons of the BEC, or if the interboson scattering length a B is significantly larger than the van der Waals length.This is the regime where most experiments would naturally be realized.However, the instability predicted in Refs.[70,71] persists in presence of modest interboson repulsion for light impurities, in a regime which is also experimentally achievable.
To arrive at the unified picture of Bose polarons, it is crucial to explicitly incorporate the interactions of the bosons in the BEC into the model, and to compare the variational Gaussianstate and double-excitation approaches.Importantly, both of these methods also incorporate the impurity-mediated interactions between bosons and the Efimov effect.To facilitate efficient numerical implementation, we parameterize the variational wave functions in terms of quantum chemistry-inspired Gaussian basis sets.
Finally, using a simple Gaussian wave function, we capture the qualitative behavior of the Bose polaron in an analytical model.This links the phenomena we see in the Bose polaron model to the paradigmatic Landau model of first-and second-order phase transitions [72].In fact, we can interpret the instability-to-crossover physics as an analog of the typical liquid-togas transition, but appearing at zero temperature.
The structure of this paper is as follows.In section 2 we introduce the model and the theoretical methods, and in section 3 we give an overview of the theory background.Then we demonstrate how we include the renormalized interboson interactions with Gaussian states in section 4. In section 5, we describe our numerical results showing the transition from the instability to the crossover regime.The analytical model that captures this behavior is presented in section 6.Finally, we conclude our work in section 7 and we provide an outlook on future directions and interesting avenues to pursue.

Hamiltonian and variational methods
We consider the problem of a mobile impurity of mass M in a homogeneous BEC of bosons with mass m and chemical potential µ B .We denote the impurity-boson interaction potential by V I B and the boson-boson interaction potential by V BB .The interactions are treated within a single-channel model, assuming the boson-impurity scattering length is tuned close to a broad Feshbach resonance.We treat the impurity in first quantization with quadrature operators R and P, and the bosons in second quantization with creation and annihilation operators b † k and bk , respectively.We set ħ h = 1.This gives the following Hamiltonian Here the unitary Ûn 0 displaces the background condensate described by a coherent-state and performs the Lee-Low-Pines transformation to transform to the reference frame of the impurity [73,74].The operator Â(x ) is dependent on variational parameters x , which are optimized to minimize the energy; |0〉 is the bosonic vacuum state with the impurity at the center of the comoving frame.In the following we set the total momentum of the system to zero, which yields the following transformed Hamiltonian.
The second term in the transformed Hamiltonian originates from the Lee-Low-Pines transformation applied to the impurity momentum operator.It is this term which gives rise to impurity-mediated interactions between the bosons [70,71], crucial for the description of the Efimov effect [75].
In this work we will compare two types of variational Ansätze.The first is a Gaussian-state (5) Here the subclass where ξ = 0 is referred to as a coherent-state (CS) Ansatz.Importantly, with Gaussian states one can compute expectation values with Wick's theorem [76], simplifying calculations.
The second wave function is a double-excitation (DE) Ansatz, given by For the double-excitation Ansatz, the case where α = 0 is often called the Chevy Ansatz, after Chevy, who first introduced an Ansatz of this kind for the Fermi polaron problem [47].

Basis set and computations
We parameterize our wave functions in terms of a Gaussian basis set.This approach is inspired by quantum chemistry where the use of Gaussian basis functions is common practice [77].
Concretely, in the Gaussian-state case, this parameterization corresponds to where the functions χ l m (σ, k) are spherical Gaussian basis functions, with spherical harmonics Y l m (θ , φ).In Appendix A we show how to compute expectation values with the Gaussian states and Gaussian basis functions.Since the polaron cloud has a smooth shape and is localized around the impurity, this approach requires much fewer variational parameters than in our previous approach in Refs.[70,71], where the wave functions were parameterized by simply discretizing φ and ξ in a spherical wave basis.Gaussian basis functions are chosen over other types of basis functions which might more closely resemble shape of the polaron cloud, because integrals over Gaussian functions give simple analytical expressions.In particular, the matrix elements of the interboson interactions generally take a complicated form, whereas for Gaussian basis functions they can still be computed analytically, at least for Gaussian potentials.
For the calculations one can either choose to keep the exponents of the Gaussian basis functions fixed, or to also treat them as variational parameters.Here we leave them fixed.The size of the smallest σ is determined by the range of the potential, and the size of the largest σ by the extent of the polaron cloud.Since these length scales are orders of magnitude different, we choose the values of σ to be spaced logarithmically.The spacing of σ is then chosen by ensuring convergence of the parameters of interest.For varying calculations we typically use between five and twenty values of σ per angular momentum mode, depending on the variational method, the observable, and the desired convergence.
In the Gaussian-state case, we optimize the variational parameters φ and ξ by using imaginary time evolution, for which the equations are derived in Appendix B. Since the problem is spherically symmetric, we can restrict ourselves to only the zero angular momentum modes for φ.For ξ, the two created particles always need to have opposite angular momenta.We solve these equations numerically with a solver based on backward-differentiation formulas [78,79], which greatly outperforms standard Runge-Kutta methods for this problem due to the stiffness of the differential equations.The stiffness originates from the interplay of the vastly different length scales of the range of the potential and the healing length of the BEC.We find that using the Gaussian basis set, qualitative and near-quantitative results can already be retrieved with a relatively small number of parameters.However, in the regimes with the strongest correlations the stiffness of the non-linear equations of motion can lead to problems reaching strict convergence when increasing the number of parameters. 2or the double-excitation Ansatz, most computation steps proceed analogously.Here one can also derive equations of motion for imaginary time evolution.However, opposed to the Gaussian-state case, the equations of motion are linear and can be solved much more efficiently by direct diagonalization.Here stiffness is thus less of an issue, and more rigorous convergence can be reached.

Interaction potentials
We model the impurity-boson and boson-boson interactions using Gaussian potentials Here L g and L U set the ranges of the potentials and g and U set the coupling strengths.The matrix elements over these interaction potentials and the spherical Gaussian basis functions are computed analytically and given in Appendix D. We fix the coupling strengths g and U to give us the desired scattering lengths a and a B , respectively.The scattering lengths corresponding to the Gaussian potentials can simply be determined by solving the two-body problem, or using the simple formulas from Ref. [80].There is no unique choice of U and g, but we take U > 0 and for g we take the smallest negative value that reproduces the desired scattering length.
The range of the boson-impurity Gaussian potential can be related to the range of the typical cold-atom van der Waals potential via the effective range r e f f .We do this at unitarity a → ∞.There, r e f f ≈ 1.4L g for the Gaussian potential, whereas for a van der Waals potential r e f f ≈ 2.8L vd w [8,81], meaning that L g ≈ 2L vdw .For modest positive scattering lengths, the results are less universal, and the repulsive interboson Gaussian potential we use here can generally not reproduce the effective range of a van der Waals potential.By fixing the absolute range relative to L g , however, we believe we still obtain representative results.Note that having finite-range interactions is crucial for the description of the Efimov effect, since the range of the interactions sets the three-body parameter and therefore the scattering length of the first Efimov resonance a − [75,[82][83][84].

Current theoretical status
By now it has been firmly established that the Efimov effect plays an important role in the Bosepolaron problem [67][68][69][70][71]85].The Efimov effect has many interesting features [75,86], but most important for this work is that it is a cooperative binding effect.One speaks of cooperative binding if the ability of particles to bind increases with the number of particles in a bound state.In the case of the Efimov effect for example, three-body bound states can be bound even if the potential is too shallow to support two-body binding.In fact, the cooperative binding of the Efimov effect also persists for more than three particles, both in the homonuclear [87,88] and the heteronuclear case [69,70,89].
In the quantum impurity case, the Efimov effect arises from an effective impurity-mediated interaction between the bosons.The cooperative binding effect comes into play in the following way: The more bosons bind to the impurity, the smaller the impurity kinetic energy per boson becomes.As a result, the bosons can use the attractive boson-impurity interaction more efficiently.In the Hamiltonian of Eq. ( 4), the second term is responsible for these impuritymediated interactions, and it originates from carrying out the Lee-Low-Pines transformation.Since the prefactor of this term is 1/M , one immediately notices that for the heteronuclear Efimov effect, the mass of the impurity is of crucial importance.The heavier the impurity is, the weaker the mediated interactions are and the more the Efimov effect is hence suppressed.
The exact role the Efimov effect plays in the Bose-polaron framework is still open to debate.Both the crossover from a polaron into an Efimov state [67-69] and a complete polaronic instability caused by the Efimov effect [70,71] have been predicted.
When a variational Ansatz of the double (or triple) excitation type is used [67, 69, 90] a crossover behavior is found by construction.This is similar when the virial expansion is performed, when truncating on the level of few excitations [68,85].In these cases the number of excitations is simply not large enough to describe the formation of bound clusters containing many particles.These results were corroborated with quantum Monte Carlo studies, for equal mass or heavy impurities [90,91].In these studies either interboson repulsion or an effective three-body repulsion originating from the use of a two-channel model [67,69,90] prevent the build-up of many excitations on the polaron.
If an arbitrary number of excitations is allowed, a divergence of the polaron energy and number of particles in the polaron cloud is predicted in the limit of non-interacting bosons [74].This persists if the interboson repulsion is treated on the level of the Bogoliubov approximation.In this case a coherent-state description yields a polaron energy given by the simple formula [74] Here a 0 is positive and arises from the finite size of the polaron cloud resulting from the Bogoliubov dispersion.In Ref. [92] a renormalization group approach was used, predicting a divergence of the polaron energy at attractive scattering lengths.In Refs.[70,71] we extended the above coherent-state approach to include pairwise interboson correlations, and three-body correlations involving the impurity.With this Gaussianstate method, we showed that for light impurities and weakly interacting bosons large Efimov clusters can be found at energies much lower than the polaron energy, rendering the polaron a metastable local minimum in the energy landscape.At a given critical interaction strength, the polaron state loses its stability and decays into the clusters.Contrary to the coherent-state or renormalization-group approaches, the polaron energy does not "smoothly" diverge as 1  a −1 −x , but a stable polaron just abruptly ceases to exist at a certain point, where the BEC locally collapses onto the impurity.The onset of this collapse was shown to be tied to a many-body shifted Efimov resonance, highlighting the importance of Efimov physics for this polaronic instability.
Even though the divergence of the particle number and energy should not occur in presence of realistic interboson repulsion, taking the limit of the interboson repulsion going to zero is not unphysical, and the Bose polaron phenomenology should therefore smoothly connect to this non-interacting limit.Hence, the natural question arises what happens to the qualitative picture of the polaronic instability for bosons with repulsive interactions.To answer this question, an approach is needed which allows for an arbitrary number of particles, but which also includes the interboson repulsion beyond the Bogoliubov level.
Coherent-state approaches treating the interboson repulsion without the Bogoliubov approximation were introduced in Refs.[93][94][95].The Born approximation for the description of the interboson repulsion is still needed though, and the resulting approach is then equivalent to using the Gross-Pitaevskii equation (GPE).However, the GPE-approach does not include the interboson correlations required to see the Efimov effect.Furthermore, whether this type of Ansatz is still applicable at strong interactions is questionable, since in this regime it might not capture well the atomic nature of cold gases [90].In presence of significant repulsion between the bosons, it can be favorable for exactly one or two bosons to be around the impurity.This cannot be captured by a coherent-or Gaussian-state approach as these describe superposition states of different particle numbers, without full control over the weights of every contribution.
In this work we address these issues by comparing the Gaussian-state and double-excitation methods on equal footing, including fully the interboson repulsion beyond the Bogoliubov and the Born approximation.In this way we aim to develop a unified understanding of the Bose polaron problem and to connect the ideas and phenomenology found from all these different approaches.

Treatment of the interboson repulsion 4.1 Interboson repulsion energy functional
As a first step, we discuss how we describe the interboson interactions with our Gaussian-state Ansatz.Since we approximate in Eq. ( 2) the background BEC with a coherent-state, no interboson correlations are included, and the interactions are treated on the level of the Born approximation.If we now include interboson correlations close to the impurity, the interboson interactions will be renormalized, unphysically resulting in a weaker effective interboson repulsion close to the impurity than in the background BEC.
A natural way to overcome this issue is to also treat the background BEC on the level of a Gaussian state.In this case the Gaussian part of the state would effectively perform a Bogoliubov rotation.However, this would severely complicate the structure of the cubic and quartic terms of the polaron Hamiltonian in Eq. ( 4).
Instead, we choose a hybrid approach.Far from the impurity we take a coherent-state wave function and describe the interactions within the Born approximation, whereas close to the impurity we keep the bare coupling, which we fully renormalize with our Gaussian-state wave function.Concretely, we achieve this by removing from the energy functional the terms responsible for renormalizing the interactions in the background BEC, which appear in the expectation values of the quadratic and cubic terms of the interboson repulsion term in Eq. (4).In the remaining quadratic and cubic terms, we replace the coupling U by U B = 8a B m πL U which gives the same scattering length on the level of the Born approximation.The quartic term is treated fully within our Gaussian-state approach.This term is the most important for the strong-coupling physics close to the impurity, since it describes the repulsion between the excitations from the BEC.The precise energy functional we use and a more detailed explanation are given in Appendix C.

Infinitely heavy impurity
To test this approach we consider first the case of an infinite-mass impurity.This is a wellstudied case [90,[93][94][95], for which there are no impurity-mediated interactions.For negative scattering lengths up to unitarity, the Born approximation for the interboson interactions is expected to hold well [94].As a result, a coherent-state approach should describe the repulsion well when the interactions are treated using the Born approximation.
In Fig. 2 we compare our Gaussian-state result using the hybrid Born description with coherent-state results using varying interboson coupling constants (see Table 1).In the CS1 approach we fully take the Born approximation, for CS2 we take the bare coupling for the quartic term, and for CS3 we take the bare coupling in all the terms.Note that the CS1 approach is equivalent to using the GPE.We set L g = L U = a B and consider a high density BEC, n 0 = 10 −5 L −3 g (corresponding to a density of ∼ 10 −14 cm −3 ).We plot our results in units of the characteristic wave vector k n = (6π 2 n 0 ) 1/3 .In Fig. 2a) we plot the energies from these methods as a function of the inverse scattering length.In Fig. 2b) we plot the density Table 1: Explanation of lines and methods used for Fig. 2.Here U 2 and U 3 stand for the coupling constants in the quadratic and cubic terms of the interboson repulsion in Hamiltonian (4) and U 4 stands for the quartic term.The coupling constant U is the bare coupling and U B = 8a B m πL U , is the coupling giving the same scattering length on the level of the Born approximation.For more details on the modified Gaussian approach, see App. C. The method CS1 is equivalent to the GPE approach.
For the characterization of the computational approaches, see Tab. 1 and the main text.The legend for the line style and color in a) and b) also apply to c),d) and e). of bosons as a function of the distance from the impurity.In Fig. 2c), d) and e) we show the number of excitations surrounding the impurity in a shell at a certain distance R: for scattering lengths (ak n ) −1 = 2, (ak n ) −1 = 0, and (ak n ) −1 = −2, respectively.In the weak-coupling regime, i.e., in the left of Fig. 2a), the quadratic interboson repulsion term is most important, and all approaches treating this term on the same footing (GS, CS1 and CS2) agree with each other within a percent.The result from CS3 already gives a difference in energy, and furthermore, in Fig. 2e) we see that this approach underestimates the extent of the polaron cloud.This is because the healing length of the BEC, which sets the extent of the polaron cloud, is too small with the unrenormalized coupling constant.
Going to stronger coupling, the short-range repulsion becomes more important.While the GS and CS1 results keep agreeing with each other within 3% up to unitarity, (ak n ) −1 = 0, the CS2 approach starts to strongly deviate from these results at (ak n ) −1 ≈ −1.Indeed, Fig. 2d) shows that also the differences in the wave functions become larger.
For positive scattering lengths larger than (ak n ) −1 ≈ 1 the Born approximation appears to break down 3 and the Gaussian-state result starts to deviate from the CS1 result.Another interpretation is the increased relevance of quantum fluctuations.In this regime, the interboson repulsion close to the impurity is more important than the repulsion in the long-ranged polaron cloud, since the density close to the impurity is more than a factor 1000 higher than the background density, see Fig. 2b).

Discussion
We now have the confidence that the Gaussian state with the hybrid Born approach properly renormalizes the interboson interactions.The Gaussian-state results namely agree with the GPE or CS1 results in the parameter regime where this approach is valid.Without renormalization of the repulsion, the results would clearly be different, as shown via the CS2 and CS3 results.Furthermore, in the region where the Born approximation is expected to break down, indeed the Gaussian-state results are different from the coherent-state results.
Properly including the interboson interactions is more complicated for the single-and double-excitation Ansätze, since these are not mean-field approaches.Restricting the number of excitations namely implicitly leads to many-body correlations.To make the most fair comparison in the next section, we as far as possible treat the double-excitation Ansatz calculation on the same footing as the Gaussian-state calculation.The explicit form of the energy functional we use for the double-excitation approach is given in Appendix C. For the boundstate regime, where the double-excitation Ansatz is at its best, the repulsion at short range is most important, and this is described properly.

Results: Polaronic instability or smooth crossover?
We now move on to the case of a finite-mass impurity, where aside from the interboson and boson-impurity interactions, there are also mediated interactions between the bosons that are generated by the impurity.We will show that this new ingredient drastically changes the behavior of the polaron.It can lead to a polaronic instability [70,71] and induce physics akin to that of first-order phase transitions.
The reason why the polaron can become unstable is simple.If the attractive impuritymediated interactions overcome the interboson repulsion, the net interactions between the bosons close to the impurity are attractive.Attractive BECs are known to collapse [96].The polaronic instability is therefore nothing else than a local, impurity-induced collapse of a part of the BEC.The mechanism of the instability is described in more detail in Refs.[70,71] as well as in Sec.6, where we show how this physics can be qualitatively captured in a simplified analytical model.
In the present section we explore the occurrence of the instability as a function of the model parameters.We start discussing the regime of light impurities, where this instability was predicted to arise [70,71].These previous studies established the presence of the instability in absence of interboson interactions, but no convincing claims could be made in presence of significant interboson repulsion.Here, we fully include a varying interboson repulsion and investigate its impact. 4Then we will study what happens as the mass ratio in the system is changed.Throughout the whole section we make a comparison of the Gaussian-state and double-excitation methods.

Light impurities and the polaronic instability
For concreteness, we consider a 6 Li-impurity in a BEC of 133 Cs.We vary both the impurityboson and the boson-boson scattering lengths.We take the ranges of the potentials, L g and L U , to be related via Here, L vd w,CsCs(LiCs) is the Van der Waals length for the Cs-Cs (Li-Cs) potential.This results in a value of L U ≈ 2.3L g .
First we test what our variational approaches predict for the interboson scattering lengths a B = 1.5L g and 2.7L g , at a density of 10 −5 L −3 g .In a typical scenario this density approximately corresponds to 10 14 cm −3 .We use the Gaussian-state and double-excitation approaches as discussed in the previous section.Moreover, we study a coherent-state approach using the Born approximation in all terms of the energy functional (CS1 in the nomenclature of the previous section).The results are shown in Fig. 3.Here we plot the energies of the various methods as a function of the inverse impurity-boson scattering length, in units of L g .In Fig. 3a) the parameters lie in the instability regime and in Fig. 3b) in the crossover regime.One immediately sees that the curves in Fig. 3 are qualitatively different from those in Fig. 2a), and that the predictions of the three methods differ by orders of magnitude.
In Fig. 3a), the three approaches strongly start to differ around the Efimov resonance, since the three-body correlations introduced by the Efimov effect are treated in a widely varying manner.Before this point, on the far left of the figure, the boson-impurity coupling is weak, and all curves coincide.The coherent-state does not capture three-body correlations at all, and continues describing the mean-field polaron.The double-excitation Ansatz predicts a relatively sharp crossover into the trimer state which appears at the Efimov resonance.Finally, the Gaussian-state Ansatz predicts a polaronic instability, marked in Fig. 3a) by the red dots, where the mean-field polaron ceases to be stable and decays into a many-particle cluster.The red dashed line indicates the energy of the many-body cluster before the polaron becomes unstable.In this regime the polaron is metastable [70,71] and not the ground-state.Since the polaronic instability corresponds to a many-body shifted Efimov resonance [70,71], it happens close to where the trimer crosses the continuum.Compared to the Efimov energy scale, the density is still relatively low and the resonance is therefore not shifted substantially (see also Fig. 6).Since all three approaches are variational, the lowest energy state best describes the ground-state, and therefore the Gaussian-state approach is most appropriate.
In the crossover regime, such as in Fig. 3b), the story is different.Here the interboson repulsion is much larger and overcomes the mediated interactions, pushing the Efimov resonance Energy(M 1 Figure 3: Polaron energies compared to cluster energies as a function of scattering length for a mass ratio M /m = 6/133, n 0 = 10 −5 L −3 g and interboson scattering lengths of a) a B = 1.5L g and b) 2.7L g , and L U = 2.3L g .The three lines correspond to the coherent-state, Gaussian-state and double-excitation Ansätze.In figure a) we observe the polaronic instability for the Gaussian Ansatz, where the energy jumps from a polaron state to a cluster state, indicated with the dots and the dotted line.The cluster state before the instability is indicated with the red dashed line, and is only shown for energies lower than the polaron energy.towards unitarity.We therefore see that the mean-field regime, where the curves coincide, extends to stronger interactions.Around unitarity, now the double-excitation Ansatz gives the best description.In this regime, no large clusters can be formed, and the polaron experiences a crossover into a trimer-like state.As evident, this behavior cannot be captured well by the Gaussian-or coherent-state methods, whose curves lie relatively close together.Only when the attractive impurity-boson interaction is increased further, more particles can bind to the impurity, and the Gaussian-and coherent-state methods outperform the double-excitation approach, as signified by their resulting energies.Note, however, that this occurs relatively far outside the universal regime.
The character of the clusters that form in the instability regime in Fig. 3a) is quite distinct from the polaron.Example wave functions of the clusters are shown in Fig. 4 for increasing interboson repulsion in absence of a background BEC (n 0 = 0) and for unitary boson-impurity interactions.We see in Fig. 4 that many particles come together within the range of the potential.A polaron cloud can also host many particles, but for a polaron state most particles are far away from the impurity, at a distance set by the healing length (see Fig. 2c-e)).
For increasing repulsion, the number of particles in the cluster rapidly decreases.At some point the bound state contains only a few particles.Here the Gaussian-state description (as the coherent-state approach) fails since it is bound to represent a superposition of states with different particle numbers, with limited control over the weights of their contributions.This is especially detrimental for the description of a bound state with exactly one or two particles.Therefore, a double-excitation Ansatz is more suitable in this regime.
As a technical side note, numerically, we find not just one, but in fact two types of stable cluster states from the Gaussian-state approach.The clusters shown so far contain both a coherent and a Gaussian contribution to the wave function, but another local minimum on the variational landscape arises when there is solely a Gaussian contribution to the wave function (φ = 0).This second type of clusters is generally higher in energy, but the real-space wave functions such as in Fig. 4 are qualitatively similar.There is one regime where the second type of cluster is lower in energy than the first type: when the particle number in the cluster goes to zero.In this case the Gaussian state just describes a superposition of the free impurity and a trimer.

Emergent "phase diagram"
We now discuss the behavior of the polaron as predicted from the variational methods, for varying densities of the background BEC.In Fig. 5 we show the energy of the polaron as a function of the inverse scattering length for several densities (indicated by the colors).We have also added the energies of the dimer (black dashed), trimer (black dashed with triangles) and the two types of cluster discussed before (thin and thick black solid) as a function of the inverse scattering length in absence of a background BEC.The different panels correspond to Gaussian-state results [Fig.5a) and b)] and double-excitation results [Fig.5c) and d)] for different values of the interboson repulsion.Note that the scale of the x-axis is different in the left and right panels.
For weak impurity-boson attraction the polaron is in the mean-field regime, and the energy should depend approximately linearly on the density [see Eq. ( 12)].We observe on the left side of all four panels that this expected result is indeed recovered, since the colored lines are spaced linearly on the logarithmic energy grid and the corresponding densities are also spaced logarithmically.The energy scale of the polaron is much smaller than the energy scale of the bound states, except at the highest densities.
In Fig. 5a) we again note the presence of the polaronic instability, but interestingly, this instability disappears for large densities.The transition from a polaronic instability to a crossover can therefore happen both as a function of density, and as a function of the interboson repulsion.The point where the instability disappears is the point where the gap between the polaron energy and the cluster energy closes.In the case of Fig. 5a), this gap is closed by increasing the polaron energy through an increase in the density.The character of the cluster is largely unaffected by the increase of the density, meaning that the formed clusters can still contain many particles even in the crossover regime.
In Fig. 5b), the repulsion is large enough that there is a crossover for all densities.For all but the largest densities, the Gaussian state does not describe this crossover well around unitary interactions [compare to Fig. 5d)] since the Gaussian-state energy actually lies above the trimer energy.
In the right two panels [Fig.5c) and d)], we see that the double-excitation Ansatz behaves qualitatively the same for all densities and for both values of the interboson repulsion: it gives a smooth crossover between the polaron and the trimer state, irrespective of the presence of larger clusters, which cannot be captured by this approach.It is interesting to notice that for increasing density, the crossover happens over a broader range of interaction strengths, and Gaussian state a B = 1.5L g 0.0 0.5 that the interaction strength where the polaron energy merges with the trimer line becomes larger.This is the opposite trend compared to the instability in Fig. 5a), where the polaronic instability happens for decreasing interaction strengths as the background density is increased.This difference can be explained as follows.For the Gaussian state, a larger particle number in the polaron cloud means that it is easier to decay into a cluster.Hence, at larger densities smaller interaction strengths are needed.For the double-excitation Ansatz, there is a competition between the polaron energy, which increases with the density, and the trimer energy, which is not strongly affected by the density.Therefore, the regime of the crossover is shifted towards larger interaction strengths as the density increases.Next, we examine in more detail the density dependence of the critical scattering length of the polaronic instability.In Fig. 5a) we can already see that the point of instability shifts to smaller scattering length as the density increases.In Fig. 6a), this critical scattering length is plotted as a function of the density for various values of the interboson repulsion, increasing from top to bottom.
For low densities, on the left hand side of Fig. 6a), the lines of critical scattering length approach the value of the Efimov scattering length [70,71].This scattering length already gives a dependence on the interboson repulsion, reflecting that also the value of the Efimov scattering length depends strongly on the interboson repulsion.
When moving to higher density, the critical interaction strength (at fixed repulsion) decreases and at some critical density the gap between the polaron and cluster energies closes.This is also visualised in Fig. 6b), where the underlying energy landscape (computed using Gaussian states) is shown as a colormap for a B = 1.2L g .Here the instability is marked with the red line.The point where the gap between polaron and cluster closes, marks the end point of the transition from the polaron to the cluster states.Beyond this density a crossover occurs.Looking again at Fig. 6a), we find that the density where the instability terminates increases rapidly with the interboson repulsion.This is because the cluster binding energy decreases with the repulsion.Therefore the gap closes already at smaller polaron energies and thus smaller densities.Aside from the red line marking the instability, Fig. 6b) also shows (white dotted line) the scattering length where the cluster energy first crosses the polaron energy.This line therefore marks the point where the ground-state of the model truly switches character.In between the white dotted and red lines we find the regime of metastability, where the polaron is not the ground-state, but where there is no decay process included in our Ansatz from the polaron state into the cluster state.In Fig 3 this is the region before the instability where the dashed line lies below the solid line.Note that the white dotted line is almost horizontal in this figure.That is because the cluster energies are on a different scale than the polaron energy, and changing the polaron energy via the density therefore does not strongly affect the crossing point of the polaron and cluster energies.
In Fig. 6c) we show a similar graph, but here we fix the density n 0 L 3 g = 10 −4.5 and vary instead the interboson repulsion.While as a function of the density, the polaron energy changes and the cluster energy is approximately constant, the opposite is true when varying the interboson repulsion, which has a much larger impact on the cluster state than the polaron.Therefore in this case the gap between the polaron and cluster energies is closed by decreasing the cluster energy.However, qualitatively the picture is the same.There is still a transition terminating in a critical point, beyond which there is a crossover.As we saw in Fig. 6a), the point of instability shifts to larger interaction strengths for increasing repulsion, following the trend of the position of the Efimov resonance.
In contrast to Fig. 6b), where the white dotted line is almost flat, here the white dotted line varies much more with the interboson repulsion than the line of instability.This is because the line of instability is determined by the smallest cluster into which the polaron can initially decay, and the line of metastability by the many-particle cluster which is lowest in energy.Since the more particles there are in the cluster, the more important their repulsion, it is not surprising that the metastability line depends more strongly on the repulsion than the instability line.
Altogether, we see that a remarkable picture emerges of a first-order transition between a polaron and a cluster, which terminates at a critical point.This is strongly reminiscent of classical first-order phase transitions, such as the phase transition of condensation of a gas into a liquid.We discuss this analogy in more detail in Sec.6 where we develop an analytical model for the Bose polaron showing the same qualitative behavior.

Mass dependence
So far we have considered the scenario of a light impurity in a BEC, where the impuritymediated interactions are particularly strong.Now we explore in more detail how the phenomena we observe manifest themselves for a wider range of impurity masses.
In Fig. 7, diagrams are shown that indicate for which values of the mass ratio and interboson repulsion the polaronic instability appears.The color code gives the corresponding critical scattering length.The background density n 0 is fixed within both panels, and given by a) 10 −6 L −3 g and b) 10 −4 L −3 g (approximately 10 13 and 10 15 cm −3 respectively).In Fig. 1b) another such colormap is shown for the density 10 −5 L −3 g .Note that on the x-axis the interboson scattering length is given in units of L U .We have chosen these units, because varying L U while keeping a B /L U fixed only gives rise to minor changes in the plots.This indicates that also the range of the interboson repulsion matters, and not just its scattering length.This is not surprising, since the range of the potentials is known to be important for the Efimov effect.
We see that the region of instability appears for light impurities and weak interboson repulsion, in the lower left of both panels.Here, the impurity-mediated interactions are the strongest and dominate over the interboson repulsion.For equal mass or heavy impurities the regime where an instability appears shrinks drastically.Furthermore, no instability appears if a B is significantly larger than L U .
Comparing Figs.7a) and b), we observe that the regime of instability shrinks as the density is increased.This is because the larger the density, the larger the polaron energy, and hence the smaller the energy gap between polaron and cluster.From the change of color one sees that the critical scattering length at which the instability occurs also becomes smaller (as also visible in Figs.5a) and 6).
g .Here L g = L U and the density is given by n 0 = 10 −5 L −3 g .

Comparison of Gaussian-state and double-excitation Ansatz
Having discussed the qualitative distinction between the regimes of the instability and the crossover, we now systematically compare the Gaussian-state and double-excitation methods.
The main aim of this comparison is to demonstrate in which parameter regime which method is best to use.In Fig. 8 we show the polaron energies from the Gaussian-state Ansatz (first column), the double-excitation Ansatz (second column), and their ratio, as a function of the mass ratio and the interboson repulsion.The different rows correspond to different bosonimpurity scattering lengths.Comparing to Fig. 7, the y-axis extends to larger mass ratios.
In the lower left of the plots, again the region of polaronic instability appears.Here the energy found from Gaussian states is obviously lower than the energy from the double-excitation Ansatz, as evident from the right column and the results we have shown before.The region in parameter space where a many-particle cluster is found with Gaussian states is more or less similar to the regime where a trimer is formed with much lower energy than the energy scale of the mean-field polaron, as seen from the double-excitation results.
For system parameters in the crossover regime, the Gaussian-state and double-excitation results appear more similar, because here at least the energy scale is the same: the scale of the mean-field polaron energy, set by the density of the BEC.For heavy impurities and a = −0.05L−1 g , in the blue region of Fig. 8a) and b), the system is truly in the polaron regime, and bound state physics does not play a crucial role.Here the ground-state energy (which is already rescaled by the reduced mass) is only mildly dependent on the mass ratio and the interboson repulsion.Despite being in the polaron regime, in the far left and far right of the plot there is a significant energy difference of up to 30 % between the Gaussian-state and double-excitation results.For low interboson repulsion the Gaussian-state performs better, because here many excitations can come close to the impurity.In contrast, in the regime of large repulsion the double-excitation Ansatz works better.Here the number of excitations is limited, and having more correlations between these few excitations and the bath is therefore more effective.Note that this argument even holds for the infinitely heavy impurity.Thus, our results show that merely the Born approximation being satisfied is not sufficient to show that a coherent or Gaussian state accurately describes the ground-state.This has the important implication that also the Gross-Pitaevskii equation, which is equivalent to a coherent-state approach, loses its validity.This confirms the analysis of Ref. [90].
When moving from the upper row of Fig. 8 to the lower rows, i.e., to stronger bosonimpurity interaction strengths, we see that the differences between the Gaussian-state and double-excitation results become progressively larger.This is because here bound-state physics starts to play a larger role.In the right-hand side of the plots a bound state containing only one or two particles is more favorable than a true polaron state.This is reflected in the much lower energy of the double-excitation result compared to the Gaussian-state result.However, towards the left-hand side of these plots the number of particles close to the impurity grows significantly, rendering, in turn, the double-excitation Ansatz insufficient.
In the strong-coupling regime where Gaussian states and the double-excitation Ansatz give similar energies [white regions in figures f) and i)], we expect that actually neither of them work very well.While the double-excitation Ansatz does not have enough excitations, the Gaussian-state Ansatz does not give enough independent control over the different particle number sectors.We identify this as the most challenging regime, which requires going beyond the double-excitation Ansatz [69] or using Quantum Monte Carlo.
Quantum Monte Carlo calculations have previously been performed mostly for equal mass or heavy impurities in presence of significant repulsion [90,91,97], and here good agreement was found with results from the double [90]-or triple-excitation [69] Ansatz.In this parameter regime these variational approaches outperform the Gaussian-state approach.In the regime of significantly lighter impurities and weak repulsion, where Gaussian states perform better, so far no Quantum Monte Carlo calculations were reported.This would be a very interesting avenue of research.

Experimental implementation
Now we turn to a discussion of the feasibility of realizing the physics predicted in this work in experiments.First, let us discuss the density range of the BEC.The density range of 10 −6 -10 −4 L −3 g , corresponds via L g ≈ 2L vd w and a typical value of L vdw ≈ 50a 0 [8] to densities of 7 • 10 12 -10 14 cm −3 .Such densities are readily available experimentally.
Next, we turn to the interboson repulsion.With the usual magnetic-field controlled Feshbach resonances, the boson-impurity and boson-boson scattering lengths can not be tuned independently.This limits experiments to the background scattering lengths in the BEC at the positions of the boson-impurity resonances.Far away from Feshbach resonances and for a stable BEC, typically a B ∼ L vd w .As a result, most experiments without large mass imbalance naturally operate relatively deep in the crossover regime (see Figs. 1 and 7).
There are other mechanisms to change the scattering lengths, which could be used in combination with the magnetic field approach to be able to tune multiple scattering lengths simultaneously [98].Optical control of the scattering length has already been observed [99,100], and Feshbach resonances tuned via radiofrequency or microwave fields [98,101] have been proposed theoretically.Even though these schemes usually lead to losses, in the Bose polaron context it might help that the interboson scattering length needs to be decreased instead of increased, and that experiments would only need to run for a short time.
From the Bose polaron experiments carried out so far, the experiments using 39 K [59, 62] incidentally have a very small interboson scattering length of a B ≈ 9a 0 .While this is most likely not in the regime of instability, it is relatively close, and therefore the formation of large clusters could be expected.However, to populate such deeply bound states one needs to go beyond the standard injection spectroscopy, because starting with a non-interacting state will give only very small overlap with such deep, many-particle bound states.
Concerning the realization of varying mass ratios between the impurity and the bosons, there are setups realizing a large mass imbalance with light impurities.Prominent examples are the Li-Cs mixtures used to study heteronuclear Efimov physics [102][103][104][105]. Unfortunately, for 133 Cs the background scattering length is usually very large, so this would most likely prevent the presence of the polaronic instability.
Another issue which arises in systems of large mass imbalance is very fast three-body loss at strong coupling.This loss arises precisely from the same mediated interactions between bosons which also lead to the polaronic instability.While this highlights another interesting connection between the field of polarons and conventional few-body physics, this may also make preparing such mixtures in a stable way at ultracold temperatures difficult.
Altogether, the physics which would be observed could differ tremendously between the different mixtures and even for different Feshbach resonances in the same mixture, see below.This gives the exciting opportunity to explore the different types of behavior.However, this also implies that caution is required when comparing different experiments, since seemingly small differences in the experimental parameters might lead to drastically different behavior.This gives a unique opportunity to explore the rich interplay between few-and many-body physics.

Discussion
We now discuss some limitations of the theoretical approaches employed in this work.

Description of the clusters and the instability
The Gaussian-state approach includes only pairwise interboson correlations and it is therefore important to discuss what would happen if higher-order correlations would be included.The polaronic instability is a cascade process, and these higher-order correlations will affect both the onset and the endpoint of the cascade.For the Gaussian-state Ansatz the onset of the instability can be viewed as a many-body shifted three-body Efimov resonance [70,71].The polaron first decays into a few-body bound state before decaying into larger clusters.If up to n-body correlations would be included then the instability would instead occur at a manybody shifted n-body Efimov resonance.However, the relevant timescale for these higher-order processes might be too slow to observe.If the order of the included correlations is as high as the number of particles of the cluster which first appears from the continuum, then one would find no instability any more, but a sharp crossover.However, for large clusters we expect that the avoided crossing in the spectrum between the polaron and the cluster would be extremely narrow.
The endpoint of the polaronic instability is a deeply bound cluster, and the structure and energies of such deeply bound clusters are not well described by Gaussian states.There are several reasons for that.First, in absence of a background condensate, the true cluster states are particle number eigenstates and Gaussian states are not.In fact, the spread in the particle number of a Gaussian state is of the order of the number of particles.Second, higher-order correlations will most likely become important in the microscopic description of many-body clusters.Even though three-body correlations between the bosons and the impurity are included, this will generally not be sufficient for bound states of more than three particles.Third, the properties of deep many-particle bound states can not be expected to be universal.Therefore the properties of these states will depend on the details of the interactions potentials and the use of simple model potentials is not warranted.
However, the detailed structures of these deeply bound clusters are unlikely to be important for experimental observables.Experimentally, once such a deeply bound cluster is formed, this would immediately lead to fast recombination losses.Furthermore, in the density regime reachable in cold-atom experiments, the qualitative mechanism of the polaronic instability should not be very sensitive to the detailed structure of the underlying clusters.In fact, unlike the wave function of the deeply bound clusters, the mechanism of the instability itself and the mediated interactions should be universal.The mediated interactions namely originate from the Lee-Low-Pines term in the Hamiltonian, which is independent of any potential.

The double-excitation Ansatz
In the crossover regime of Figs. 1 and 7, far enough away from the polaronic instability, we believe that the ground-state energy is well described by the double-excitation Ansatz.Possible corrections can be accounted for by an extension to the triple-excitation Ansatz [69].Therefore, we believe that regarding the ground-state energies, the picture sketched for equal mass and heavy impurities is accurate in Refs.[69,90].
However, whether the wave function of the resulting state is also well described, remains an open point of discussion [74].For example, in the regime where the ground-state of the system is a trimer, the double-excitation Ansatz will 'invest' its excitations to describe this trimer state.However, it is not unlikely that in reality this molecule will again induce its own polaron cloud.The double-excitation Ansatz is not capable of describing this behavior, since its excitations are already used up in the description of the bound state itself.While for the ground-state energy this may not have a significant effect, for the wave function and the description of the full excitation spectrum it can certainly be of importance.
Another shortcoming of our implementation of the double-excitation Ansatz is its description of the intermediate-coupling Bose polaron in presence of explicit interboson repulsion, which especially gives problems in the cubic term (see App. C).When the repulsion is described on the Bogoliubov level there is no problem, and the results from the double-excitation Ansatz agree with perturbation theory up to third order [67].However, when explicitly taking into account the repulsion beyond the Bogoliubov approximation, consistently describing the repulsion with the double-excitation Ansatz would require an accurate description of the interactions in the background BEC.In this case it is more natural to describe the polaron cloud on the same footing as the background BEC, as done using a coherent-or Gaussian-state Ansatz.
Moreover, the double-excitation Ansatz brings with it the subtle problem of how to count the number of particles in the polaron cloud.One would expect that the double-excitation Ansatz gives rise to at most two excitations.This is indeed true if one does not allow the doubleexcitation Ansatz to have a component in the mode of the BEC.However, with a coherent-state or Gross-Pitaevskii approach, one does not have this restriction.If indeed the double-excitation Ansatz is also allowed a component in the mode of the BEC, then the cross-terms with the BEC will give rise to equally large particle numbers in the polaron cloud as in the coherent-state case.This could potentially explain the qualitative difference between the large number of particles the coherent-state approaches predict in the polaron cloud [93][94][95] and the apparent success of variational approaches with only few excitations [67,69,90].

Narrow Feshbach resonances
In this work we use a single-channel model for the boson-impurity interactions, which is best suited to describe broad and isolated Feshbach resonances [8].In Refs.[67,69] a two-channel model has been used, which is also applicable to narrow resonances.In this case, the multichannel nature of the interactions can lead to an effective three-body repulsion.This will have a similar effect as the intrinsic interboson repulsion and therefore help to suppress the instability.This may therefore lead to a shift of the boundary of the instability region in Figs. 1  and 7.

Analytical model
The form of Fig. 6b) and c), with its first-order transition ending in a critical point, followed by a smooth crossover, is remarkably similar to well-known diagrams of first-order phase transitions such as the liquid-gas phase transition.To strengthen this connection we attempt to understand the Bose polaron phenomenology in simpler terms.To this end, we develop a analytical model to qualitatively reproduce the key features of Fig. 1 and Fig. 6. Surprisingly, we find that we can achieve this with a much simplified Gaussian-state Ansatz.As we will see, even when restricting the variational Ansatz to only l = 0 and l = 1 angular momentum modes, and only a single Gaussian basis function per angular momentum mode, the model already qualitatively reproduces the important physics.
We thus start the derivation of the analytical model by writing

Coherent states and effective scattering length
To build our understanding, we start by first studying a coherent-state Ansatz and omitting the interboson interactions in the Hamiltonian.We replace the real-space Gaussian bosonimpurity potential by a short-range separable interaction with a Gaussian cutoff function in momentum space, The scattering length a for this potential is related to the coupling constant g as Using only the coherent-state part of our Ansatz (i.e.ξ = 0), one finds the energy as a function of φ and σ φ to be given by where T φ and V φ are the expectation values of the kinetic and interaction energies Eq. ( 19) can be trivially minimized with respect to φ and σ φ .The resulting value of σ φ = 5σ g is independent of any of the other parameters.This leads to where Note that this expression of the energy is remarkably similar to the energy found from mean field theory assuming a weakly repulsive BEC within the Bogoliubov approximation [74].In the case of Eq. ( 12) from Ref. [74], the origin of the shift a 0 of the scattering length in the denominator can be traced to the interboson repulsion limiting the size of the polaron cloud.
In our case it is the exponent of the Gaussian basis function that limits the size of the cloud.Note that the interboson repulsion is not included yet in the analytical model up to this point.
Including the repulsion would prevent the divergence of the energy in Eq. ( 23).
If we want to compare our analytical result with the full model we can define an effective scattering length This effective scattering length diverges when a bound state appears in our model.We can also write g in terms of a eff as With this replacement, the coherent-state energy for the polaron without background repulsion is recovered perfectly.

Gaussian states and the polaronic instability
Having considered the coherent-state case, we now include also the Gaussian part as in Eq. ( 16) into the wave function, still without including the interboson repulsion.The additional variational parameters are ξ 0 , ξ 1 , σ 0 and σ 1 .We keep σ φ =5σ g fixed, since the coherent part is the dominant part in the polaron regime.Expanding the energy functional up to quadratic order in φ 2 , ξ 0 and ξ 1 , we find In this expression, the quantities T 1 and T Lφ are given by The quantities T 0 , V 0 and T L0 are defined similarly to T φ , V φ , and T Lφ by substituting σ φ by σ 0 .The T Lφ and T L0 terms originate from the Lee-Low-Pines term in the Hamiltonian.The terms S 0 and S 1 are overlap integrals, given by Minimizing the energy functional (26) with respect to σ 0 , σ 1 , ξ 0 and ξ 1 now yields σ 0 = σ φ = 5σ g and σ 1 = 15 2 σ g .This implies T φ = T 0 and T Lφ = T L0 = T L .For ξ 0 and ξ 1 one then finds The parameters ξ i characterize the Gaussian part of the wave function.Thus, already from Eqs. ( 31) and ( 32) we can see a sign of the instability, which will occur when the denominator of ξ 1 vanishes and hence ξ 0 and ξ 1 diverge.From Refs.[70,71] we know that in the lowdensity limit this happens at the Efimov scattering length a − , where the three-body Efimov bound state arises.Hence, if we derive the value of a eff where the divergence occurs, we can extract the value of a eff,− in our model with the simplified wave function.This is given by: The scattering length a eff,− is negative.We thus see that even in our simple model a threebody bound state appears from the continuum in a Borromean way, i.e., it arises before a two-body bound state is possible.Furthermore, the linear scaling of a eff,− with the length scale of the potential σ g , behaves according to the expectations of van der Waals universality [82][83][84].Finally, note that for a light impurity M 2 µ 2 r → 1 and for a heavy impurity This implies that a eff,− → −11.2 σ g for light impurities and a eff,− → −∞ for heavy impurities.The limits of the mass-dependence of a eff,− are therefore also physical.Quantitatively, however, one should note that the realistic mass dependence of a − is stronger than found here.Now we can plug in the results for ξ 0 and ξ 1 to obtain an energy functional in terms of only φ: The structure of this equation with a linear, quadratic and quartic term in φ is reminiscent of the paradigmatic Landau model for phase transitions, where φ would correspond to the order parameter.Our scenario, with a positive quadratic term and a negative quartic term, corresponds to the case of a first-order phase transition.Here the linear term, which depends on the density of the background BEC, adds an effective external field to the description (similar to a external magnetic field in the theory of phase transitions in magnetic materials).Note that all the terms depend explicitly on the boson-impurity interaction strength.as a function of the normalized "order parameter" φ for several effective scattering lengths.The magnitude of the scattering length increases from top to bottom.The mass ratio here corresponds to Li-Cs and n 0 σ 3/2 g = 1.25 10 −5 .
The energy functional (35) is plotted in Fig. 9 for different values of a eff .Here S 0 , as defined in Eq. ( 29), serves to normalize the Gaussian basis function.The dashed lines indicate the combined contribution of the linear and quadratic parts, while the solid lines show the full result from Eq. (35).
For small a eff the function has a minimum at small φ corresponding to the polaron state.In this regime, the quartic term plays no role yet.For increasing a eff the value of φ at the polaron minimum increases.Therefore the quartic term becomes more and more important.As evident from Eq. ( 35), the quartic term is also directly dependent on a eff , further enhancing its importance for increasing a eff .
At some point the quartic term overcomes the quadratic term and the polaron minimum disappears: the polaron becomes unstable.Because no interboson repulsion is included to stabilize the energy functional, φ will grow indefinitely beyond this point.Our model breaks down in this limit since ξ 0 and ξ 1 cease to be small parameters, and higher order terms in ξ will be important in Eq. (26).
The point of polaronic instability does not correspond to the point of a phase transition in the Landau model, because a phase transition in the thermodynamic sense happens when the ground-state of the system changes its character.This would be at the point where any cluster becomes lower in energy than the polaron for the first time.The metastability of the polaron can instead be interpreted as a form of hysteresis.The point of the true phase transition is not properly described in the energy functional (35) and requires to improve the model further as we will discuss in the following.
At the point of instability, both the first and second derivative of the energy with respect to φ vanish and we can use this to find the density and scattering length determining the "stability boundary" of the polaron state.This critical density, as a function of the scattering length, is given by Finding the inverse equation for a eff as a function of n 0 can be done by solving a third-order polynomial, yielding an analytical, but lengthy expression.Remarkably, the only dependence on the mass in equation Eq. ( 36) is via a eff,− , see Eq. (33).

Including interboson repulsion
To obtain an even clearer analogy to the theory of phase transitions, we now include the interboson repulsion to stabilize the cluster states.The only interboson repulsion term which qualitatively matters for the description of the behavior at strong coupling is the quartic term in the Hamiltonian of Eq. ( 4).The quadratic and cubic term are mostly important to determine the shape of the polaron cloud at long distances, which we do not attempt to describe in this simplified model.Therefore, we only keep the quartic term.Note that this is the opposite approach compared to what is usually done in the Bogoliubov approximation.
To see all the relevant effects we need to expand the energy functional to the third instead of only second order in φ 2 , ξ 0 , and ξ 1 .At this point we do not optimize the exponents of the Gaussian basis functions again, but simply use the optimized exponents obtained in the previous step.This yields the energy functional Here we have defined The numbers U 00 , U 10 , U ′ 10 and U 11 are defined in the Appendix in Eqs.(D.27-D.30).Minimizing Eq. ( 37) with respect to ξ 0 and ξ 1 gives Taking the limit φ → 0, we can find a new value for a eff,− as a function of U, When we again expand the energy functional up to quartic order in φ we retrieve Eq. ( 35) but with a eff,− replaced by a eff,− (U).To recover the effect of the stabilization of the polaronic collapse, we need to go to higher order in φ.We find Since ξ 0 has φ 2 in multiple arguments of the denominator and the numerator, this does not allow a simple expression in terms of a eff,− (U).

Note further, that the quartic interboson repulsion term
from Eq. ( 37) has been absorbed into the quartic term originating from the Gaussian part of the state.44), where the energy is plotted in units of µ r n −2/3 0 as a function of the normalized φ for several scattering lengths.The magnitude of the scattering length is indicated by the colorbar and it increases for the graphs from top to bottom.The mass ratio chosen here corresponds to Li-Cs and n 0 σ 3/2 g = 1.25 10 −5 .The three subfigures correspond to three different values of the interboson scattering length from left to right given by a B = 0.7, 0.9 and 1.1 L U , respectively.Here L U = 4.5 σ g .In Fig. 10 the value of the energy functional of Eq. ( 44) is plotted as a function of the order parameter φ S 0 for several boson-impurity (indicated via the color of the lines) and interboson scattering lengths.
For weak repulsion [Fig.10a)] we find a double well picture of a shallow well corresponding to the polaron and a deeper well corresponding to the Efimov cluster.Quantitatively, the cluster has a much weaker binding energy than in the original model, mainly because the exponents of the Gaussian basis functions are optimized for the polaron and not the clusters.
Qualitatively, the physics is very similar as for the full model.In Fig. 10a), the first-order transition is clearly apparent.It occurs when the well corresponding to the cluster on the right becomes lower in energy than the well associated with the polaron state.However, the polaron first remains a metastable local minimum, even though the cluster state is lower in energy.Then, at some critical scattering length, the barrier protecting the polaron disappears and there is a sudden transition from the polaron into the cluster state.This can still be interpreted as a form of polaronic instability.
If the interboson repulsion is increased, the double well picture no longer applies.In Fig. 10b) the energy functional is shown close to the critical point.Here we see that there is no first-order transition or polaronic instability any more, but still the ground-state changes rapidly in character for some critical scattering length.Finally, for even stronger repulsion, in Fig. 10c), no sign of a transition remains, and we are deeply in the smooth crossover regime.
In Fig. 11 we capture this behavior in a single figure, by showing the polaron energy as a function of the boson-impurity and boson-boson scattering length.Indeed, the figure is remarkably similar to Fig. 6c), showing how well our analytical model captures this behavior.Note that in the grey area in the bottom left, a eff < a eff,− (U).Here our model breaks down.In particular, in this figure we clearly see the line of first-order transition, where the polaron ceases to be the ground-state, and the line of instability, where the polaron becomes completely unstable.These two points merge in the critical point, where the energy functional takes the form as shown in Fig. 10b).At this critical point, where the line of first-order transition terminates, the phase transition turns into a second-order one.
g .The red line indicates the polaronic instability and the white dashed line indicates the point where the first-order phase transition happens, where the cluster energy crosses below the polaron energy.Both of these lines end at the critical point.In the grey region in the bottom left, our analytical model is not applicable.Fig. 11 as a whole, as well as Fig. 6b) and c), are remarkably reminiscent to the phase diagram corresponding to the gas-liquid phase transition.In this analogy the polaron state corresponds to the gaseous state and the cluster state to the liquid state.The gas-liquid transition is also a first-order phase transition, up to a critical point, after which the state is called a supercritical fluid.In the regime where the polaron state is metastable, this would be analogous to a supercooled gas: a gas cooled below its condensation point.Note, however, that in the polaron case, this phase diagram is realized in the quantum regime at zero temperature.
As a final result, in Fig. 12 we plot the "Bose polaron phase diagram" as a function of the mass ratio and the interboson repulsion, based on our analytical model Eq.(44).We see that the stability diagram is remarkably similar to Figs. 1 and 7, even on a quantitative level. 5This shows that the form of the Bose-polaron phase diagram is remarkably robust.

Conclusion and outlook
We have characterized the behavior of the attractive Bose polaron, i.e., the ground-state of a mobile impurity interacting with a Bose-Einstein condensate, across a large parameter regime of boson-impurity and boson-boson scattering lengths, BEC densities and impurity-to-boson mass ratios.Thus, we have brought together the qualitatively different results from several studies into a single, unified theoretical picture.To this end, we have compared two state-ofthe-art variational methods: the Gaussian-state and double-excitation approaches.We developed a computationally efficient technique by expressing the variational functions in terms of a set of spherical Gaussian basis functions.
We found that the polaron experiences an instability as predicted in Refs.[70,71] for weak repulsion and light impurities.This instability turns into a smooth crossover for larger g .In the red area of the diagram the polaron undergoes an instability as the boson-impurity interaction strength is swept across a Feshbach resonance.The critical scattering length a c of the instability is indicated by the colormap.In the white regime the polaron undergoes a smooth crossover instead.repulsion or heavier impurities.Most of the experiments will naturally be in the crossover regime, but the instability regime should be experimentally accessible with light impurities in BECs with a small interboson scattering length.
We developed a simple analytical model capturing the phenomenology of the instability and crossover.From this model it becomes apparent that the physics of attractive Bose polarons can be understood in the language of the Landau model of first-order phase transitions.It is in fact strongly reminiscent of the gas-liquid phase transition, where the polaron state is analogous to the gaseous phase and the cluster state to the liquid phase.Clearly, in the case of a single impurity we cannot truly speak of a phase transition in our model, but it is likely that this will turn into a proper quantum phase transition when considering a finite density of impurities.We have furthermore shown that the stability diagram of the polaron is reproduced surprisingly well by the simple analytical model.A topic inviting further study would be a more detailed characterization of the critical behavior at the critical point where the phase transition turns second order.
Given that the properties of the ground-state differ strongly across the parameter regimes, an interesting question remains how these effects would manifest themselves in the dynamics and the excitation spectrum of the polaron.Especially interesting would be to see whether some resonant behavior occurs at the onset of the polaronic instability, in analogy to the Efimov resonance [70,71].Dynamical calculations have so far been carried out mostly using (extended) Gross-Pitaevskii equations [106] or equivalently, with a coherent-state approach [74,107,108].To see the effects mentioned above one would need to explicitly allow for the correlations between the bosons induced by the impurity.In principle, the Gaussianstate approach can also be applied for real-time evolution, but this is computationally more challenging due to the more oscillatory nature of the wave function.For this reason, the parameterization of the variational functions used here in terms of Gaussian basis functions, is likely not optimal for real-time evolution.
From dynamical calculations also the polaron spectral function can be extracted, which can then be compared to experimental data.For a good comparison, it will also be important to assess what is the role of recombination.This will lead to spectral broadening especially for many-body bound states, and this might therefore make it difficult to experimentally observe such states.
Furthermore, also interesting effects are predicted to appear in the finite temperature behavior of the Bose polaron [109][110][111][112], although several approaches have again found conflicting results.Therefore, it would be a promising avenue of study to see whether a unified model can also be developed capturing the finite-temperature properties for varying mass ratios and background repulsion.
In the context of Bose-Fermi mixtures, mediated interactions by fermions in Bose-Einstein condensates have been studied already both theoretically [113][114][115] and experimentally [32,33] in a regime of finite fermion density.However, in our work the fermion-mediated interactions are caused by a different mechanism.Namely, they arise from the Efimov effect, instead of the many-body effects originating from the Fermi surface of the fermions.Since we show that the impurity-mediated interactions play a profound role in the Bose polaron model, it can be expected it will also be important in the general description of finite-density mixtures at strong interactions.We hope our work can provide a starting point for fruitful studies in this direction.In particular, we hope that further inspiration can be drawn from quantum chemistry to also tackle the bipolaron problem.
Bose polarons have recently been observed for the first time in two-dimensional semiconductors coupled to a microcavity [116].Even though the standard Efimov effect does not exist in two dimensions, it is still possible that in certain conditions bound states with more than two particles can form.Therefore, similar phenomena as we have predicted here in the context of cold atomic gases might well be found in these solid state platforms.
Note added.During the preparation of the manuscript a related work has appeared [117].Here the authors introduce another variational approach to the Bose polaron problem, which allows to extract excited states in the spectrum.However, the authors include only limited interboson correlations, meaning that the Efimov effect, which is underlying most of our results, is not captured.

B Imaginary time evolution
We optimize the energy of our state using imaginary time evolution.Here we derive the equations of motion.We use McLachlan's variational principle where E is the energy, given as where we assume the norm of the state to be 1.The energy term is needed in Eq. (B.1) to ensure conservation of the norm.Now the derivative ∂ τ can be replaced by: Inserting this into Eq.(B.1) gives: The parameter vector x consists of N , N * , φ, φ * , ξ and ξ * .The derivatives ∂ x i |ψ〉 are given by: Let us now denote x ′ as the vector x without N and N * .If we take x j = N in Eq. (B.4) we arrive at Multiplying this equation by |N | 2 N and adding it to the complex conjugate, we exactly recover the equation for the norm (B.12)This shows that indeed the norm is conserved.For the other equations: Inserting now Eq.B.11 gives: (B.15) where η * j = ∂ φ j E. In practice, taking a simplified equation for φ gives the same result with comparable computational efficiency.
The energy E is most simply expressed in terms of F and G because of Wick's theorem.If we thus define we can express these derivatives in terms of ξ via This yields the total expression Inserting this into Ref.B.17 and multiplying with (SX S) −1 from the left (transposed) and right, yields This derivation can straightforwardly be applied also to the case of real-time evolution.The equations of motion for ξ for real-time evolution are also given for a general case in Appendix G of Ref. [76].In other works using the Gaussian-state variational approach the equations of motion are usually formulated in terms of the covariance matrix Γ [71,76,119].These equations of motion can be extracted directly from the equations for ξ and ξ * .However, we find that for imaginary time evolution, it is numerically more stable to work with ξ instead of Γ .

C Treatment of the interboson repulsion
The interboson repulsion term of the Hamiltonian is given by ĤU The Gaussian-state expectation value with respect to this term given by For the ground-state scenario, φ is real, meaning that φ * = φ.The same holds for F and G.
Note that the F -terms typically dominate over the G-terms in the polaron regime, since F ∼ ξ and G ∼ ξ 2 , and ξ is generally small.Furthermore, F is usually negative, so that the φ 2 F * terms give a negative energy contribution, whereas G is always positive.Directly using this energy functional gives rise to problems, since the n 0 F and n 0 φ F terms give rise to a negative energy contribution even in absence of interactions between the impurity and the bath.This is caused by the fact that the background BEC is described by a coherentstate, which treats the interactions in the BEC on the level of the Born approximation.Because the Gaussian-state approach is now capable of renormalizing the interactions surrounding the impurity, this gives rise to an artificial three-body potential.We resolve this be removing these terms from the energy functional.
Consequently, the interboson interactions between the bosons in the polaron cloud and the bosons in the BEC are not renormalized by the Gaussian-state Ansatz.As we will show below (see Fig. 13), the best results are obtained when we replace the corresponding coupling constants in the energy functional by the coupling constant U B = 8a B m πL U , just like is done for the energy of the background BEC.This gives the correct scattering length on the level of the Born approximation.As a result, all the interactions which are not renormalized are treated on the level of the Born approximation.
This procedure leads to the energy functional Importantly, the interactions between the bosons in the polaron cloud, corresponding to the lower two lines in Eq. (C.2), are still fully renormalized by the Gaussian-state Ansatz.In the strong-coupling regime, these terms are the most important (see e.g.Fig. 2b)), and also the repulsion in cluster states in absence of a background BEC is still described fully beyond the Born approximation.
For the double-excitation Ansatz In this case it is more ambiguous how to take the Born approximation, since this is not a mean-field approach.The best strategy depends on the parameter regime.To be most consistent with the Gaussian-state case we choose to use the following energy functional This energy functional is chosen to reproduce the mean-field result in the weak-coupling limit.The quartic term is still described exactly, which is most important for the main part of our work.
To understand better the role of the Born approximation and the comparison between the Gaussian-state and double-excitation methods, we show in Fig. 13 the corresponding polaron wave functions for an infinitely heavy impurity and an intermediate negative boson-impurity scattering length (ak n ) −1 = −2 (the parameters are such as in Fig. 2e)).This is the polaronic regime, and if the interboson repulsion would be properly treated, the results should be in close agreement with each other.We first observe that the Gaussian-state result with the energy functional of our choice (red solid) agrees well with the coherent-state result including the Born approximation (black dashed).If we do not take the Born approximation in the quadratic and cubic term of the Gaussian-state energy functional (red dash-dotted) we find that the healing length of the BEC is not correctly captured and that the results are more similar to the coherent-state result without the Born approximation (black dotted).
The double-excitation result is close to the Gaussian-state result, but still noticeably different.In this parameter regime this is an artefact of our way of dealing with the interboson repulsion, which gives problems in the cubic term of the double-excitation Ansatz, especially in the polaronic regime.In the weak-coupling limit, where only the quadratic term matters, the Gaussian-state and double-excitation results agree with each other.For the parameters of Fig. 13 the energy difference between the Gaussian-state and double-excitation Ansatz is approximately 3%.Even though this is significant, this difference is still very small compared to the qualitative differences between the two approaches observed in Sec. 5. Furthermore, in the bound-state regime, it is the quartic term in the energy functional which is strongly dominant, and this is described fully for both approaches.Hence, we do not expect our approximations in the treatment of the interboson repulsion to have an impact on our conclusions for the strong-coupling Bose polaron.

D Computation of Gaussian integrals
The expectation values with respect to the kinetic, boson-impurity interaction and LLP-terms yield relatively simple expressions.
For the kinetic energy we find = π(l + 1) [(2l + 3)!!] 2 (2l + 1)M 2 2l+6 (σ l+1,i + σ l, j ) 5/2+l (σ l+1,p + σ l,q ) 5/2+l .(D.4) For the interaction term with the Gaussian boson-impurity interaction potential (10) it is most convenient to carry out the integral in real space, and define the Fourier transform χ of the Gaussian orbitals.Then we find χlm (σ, r ) = 1 (2σ) l+3/2 Y l m (θ , φ) exp −  Here we use the combined indices j 1 = ( j 1 l 1 m 1 ), where the first index labels the exponent α j 1 , for which α i = 1 4σ i .The integral is nontrivial because the interaction potential depends on the distance between the bosons (r 1 − r 2 ).To perform the angular and radial integrals we need to write this in a different form.Indeed, following Ref.[120] we can decompose the term according to   The expression can be simplified when taking our symmetries into account.For the coherent part for example, l 1 = l 2 = l 3 = l 4 = 0.In this case we find 2 , 3 2 ; 3 2 ; Now we have all we need to find the expressions for U 00 , U 10 , U ′ 10 and U 11 from Eq. (37).For this, the integrals need to be multiplied by the prefactor U 2L 2

U
. Furthermore, we need to add the prefactors following the Fourier transform to real space, see Eq. (D.5).Using this, we can fill in the above equations for the case of our analytical model:

Figure 1 :
Figure1: a) Energy of the Bose polaron as a function of the inverse impurityboson scattering length 1/a in the regime where it undergoes an instability (bold) or crossover (dashed) as a function of the interaction strength.In the instability regime, the polaron state becomes unstable at the scattering length indicated by the dot.At this point the polaron decays into a cluster.The energy of the cluster before the instability is drawn with the dash-dotted line (only shown in the regime where it is lower than the polaron energy).The bold and dash-dotted lines are computed using a Gaussian-state Ansatz, and the dashed line using a double-excitation Ansatz (see Sec. 2).b) Stability diagram indicating whether an instability or crossover will occur as a function of the impurity-boson mass ratio M /m and the interboson repulsion scattering length a B .In the lower left part of the diagram, the scattering length of the instability a c is indicated by the colormap.The crosses indicate the parameters corresponding to the lines in subfigure a).For both a) and b) scales characterizing the interaction ranges of the boson-impurity and boson-boson potential, L g = L U = 2L vd w are chosen to be equal (see Sec. 2.3).For the density of the BEC we have chosen a typical value of n 0 = 10 −5 L −3 g (approximately 10 14 cm −3 ).
Figure1: a) Energy of the Bose polaron as a function of the inverse impurityboson scattering length 1/a in the regime where it undergoes an instability (bold) or crossover (dashed) as a function of the interaction strength.In the instability regime, the polaron state becomes unstable at the scattering length indicated by the dot.At this point the polaron decays into a cluster.The energy of the cluster before the instability is drawn with the dash-dotted line (only shown in the regime where it is lower than the polaron energy).The bold and dash-dotted lines are computed using a Gaussian-state Ansatz, and the dashed line using a double-excitation Ansatz (see Sec. 2).b) Stability diagram indicating whether an instability or crossover will occur as a function of the impurity-boson mass ratio M /m and the interboson repulsion scattering length a B .In the lower left part of the diagram, the scattering length of the instability a c is indicated by the colormap.The crosses indicate the parameters corresponding to the lines in subfigure a).For both a) and b) scales characterizing the interaction ranges of the boson-impurity and boson-boson potential, L g = L U = 2L vd w are chosen to be equal (see Sec. 2.3).For the density of the BEC we have chosen a typical value of n 0 = 10 −5 L −3 g (approximately 10 14 cm −3 ).

Ĥ0 1 (2π) 3 d
(r − R) b † r br + 1 2 r ′ r V BB (r ′ − r ) b † r ′ b † r br ′ br .(1) Here E b g is the energy of the background BEC without the impurity.We denote r = d 3 r and k = 3 k.The chemical potential is set to the mean field value µ B = n 0 U B , where n 0 is the density of the BEC and U B the coupling constant on the level of the Born approximation.To approximate the ground-state of the Hamiltonian we consider a variational Ansatz of the type |ψ〉 = Ûn 0 ÛLLP Â(x )|0〉 .

Figure 2 :
Figure 2: Properties of an infinite mass impurity immersed in an interacting BEC of density n 0 = 10 −5 L −3 g and for a B = L g = L U .a) Polaron energy as a function of the inverse boson-impurity scattering length.b) Bosonic density for three different scattering lengths from Gaussian states as a function of the distance to the impurity.c-e) Number of additional particles at distance R from the impurity [see Eq. (13)] for c)(ak n ) −1 = 2, d) (ak n ) −1 = 0, and e) (ak n ) −1 = −2.For the characterization of the computational approaches, see Tab. 1 and the main text.The legend for the line style and color in a) and b) also apply to c),d) and e).

Figure 4 :
Figure 4: The wave function of the deepest bound clusters at unitarity (i.e., a −1 = 0) found from Gaussian states for the mass ratio M /m = 6/133 and for n 0 = 0.The colors of the lines indicate the value of the interboson repulsion.

Figure 5 :
Figure 5: Polaron energies (colored lines) as a function of the inverse bosonimpurity scattering length for a light impurity (M /m = 6/133), calculated using a,b) the Gaussian-state Ansatz, c,d) the double-excitation Ansatz.In panels a) and c) a B = 1.5L g , in panels b) and d) a B = 2.7L g .The color of the lines indicates their corresponding background density.In all panels the black dashed lines indicate the dimer and trimer (with triangles) energies.The black solid lines indicate the energies from the two types of cluster from the Gaussian-state approach at n 0 = 0.In figure a) the point of polaronic instability is indicated with the pairs of dots and dotted lines.

Figure 6 :
Figure 6: a) Critical scattering length a c of the polaronic instability as a function of the density n 0 , for various values of the interboson scattering length a B .At the square endpoints of the line the gap between the polaron and cluster states has closed, and the instability turns into a crossover.b,c) Polaron energy as a function of scattering length for b) varying density and fixed a B = 1.2L g and c) varing a B and fixed density n 0 L 3 g = 10 −4.5 .The red lines indicate the critical scattering length of instability, such as in figure a).The white dotted line indicates the scattering length where the energy of the cluster crosses the polaron energy.Between the dotted and dashed lines the polaron is thus metastable.The mass ratio is chosen to be M /m = 6/133 and L U = 2.3L g .

Figure 7 :
Figure 7: Stability diagrams of the polaron for densities a) n 0 L 3 g = 10 −6 and b) n 0 L 3 g = 10 −6 as a function of the impurity-boson mass ratio M /m and the interboson repulsion scattering length a B .In the lower left part of the diagram, when the boson impurity scattering length is swept across the Feshbach resonance, the polaron will experience an instability at scattering length a c indicated by the colormap.In the upper right part the polaron will smoothly cross over into a small cluster.For this plot L g = L U = 1.

Figure 8 :
Figure 8: Colormaps of the polaron energy as a function of the mass ratio M /m and interboson scattering length a B from the Gaussian-state Ansatz (first column)and double-excitation Ansatz (second column).In the third column the ratio of the energies from these approaches is shown.The different rows correspond to different inverse boson-impurity scattering lengths, given by a-c) −0.05L −1 g , d-f) 0L −1 g , f-i) 0.05L −1 g .Here L g = L U and the density is given by n 0 = 10 −5 L −3 g .

Figure 10 :
Figure 10: Energy functional of Eq. (44), where the energy is plotted in units of µ r n

Figure 11 :
Figure 11: Colormap of the energy of the polaron as a function of the effective impurity-boson scattering length a eff and the interboson repulsion scattering length a B as obtained from our analytical model.The density n 0 σ 3/2 g = 2.5 10 −5 and

Figure 12 :
Figure 12: Stability diagram of the polaron as a function of the impurity-boson mass ratio M /m and the interboson repulsive scattering length a B , as obtained from our analytical model.The density n 0 σ 3/2 g = 1.25 10 −5 and L U = 2σ 1/2g .In the red area of the diagram the polaron undergoes an instability as the boson-impurity interaction strength is swept across a Feshbach resonance.The critical scattering length a c of the instability is indicated by the colormap.In the white regime the polaron undergoes a smooth crossover instead.

Figure 13 :
Figure 13: Polaron wave functions for an infinite mass impurity, n 0 L 3 g = 10 −5 , a B = L U = L g and (ak n ) −1 = −2, such as in Fig. 2e).The lines labeled with GS, CS1, and CS3 correspond to the parameters listed in Tab. 1.The GSU-line corresponds to the Gaussian-state Ansatz, but with U 2 = U 3 = U instead of U B , and the DE-line corresponds to the double-excitation Ansatz with the interboson repulsion energy functional as in Eq. (C.6).