From Efimov Physics to the Bose Polaron using Gaussian States

Since the Efimov effect was introduced in 1970, a detailed theoretical understanding of Efimov physics has been developed in the few-body context. However, it has proven to be challenging to describe the role Efimov-type correlations play in many-body systems such as quenched or collapsing Bose-Einstein condensates (BECs). To study the impact the Efimov effect can have in such scenarios, we consider a light impurity immersed in a weakly interacting BEC, forming a Bose polaron. In this case, the higher-order correlations are localized around the impurity, making it more feasible to develop a theoretical description. Specifically, we employ a Gaussian state variational Ansatz in the reference frame of the impurity, capable of both capturing the Efimov effect and the formation of the polaron cloud. We find that the Efimov effect leads to a cooperative binding of bosons to the impurity and the formation of a many-body bound state. As a result, the polaron is not the ground state, but rendered a metastable excited state which can decay into these Efimov clusters. While this decay is slow for small interaction strengths, it becomes more prominent as the attractive scattering length increases, up to the point where the polaron becomes completely unstable. This critical scattering length can be interpreted as a many-body shifted Efimov resonance, where the scattering of two excitations of the bath with the polaron can lead to bound state formation. Compared to the few-body case, the resonance is shifted to smaller attractive scattering lengths due to the participation of the polaron cloud in the cooperative binding process. This corresponds to an intriguing scenario of polaron-assisted chemistry, where many-body effects lead to enhanced signal of the chemical recombination process, which can be directly probed in state-of-the-art experiments.


I. INTRODUCTION
To describe the properties of a many-body system, an accurate understanding of the relevant interactions and correlations between its microscopic constituents is crucial. However, even when the few-body physics is understood, extracting the emergent properties of the system as a whole is not a simple task. A good example is the Efimov effect [2,3]. This fascinating few-body effect was predicted by Efimov in 1970 in the context of nuclear physics. Efimov showed that there exists an infinite series of three-body bound states appearing close to the unitarity point of a Feshbach resonance with binding energies obeying a geometric scaling law. It took more than 30 years to demonstrate this effect experimentally, first in cold atomic gases [4] and later in Helium [5]. After that there were more experimental observations and a good theoretical understanding of the few-body physics has been developed [3]. However, the task to understand the effect that Efimov-like correlations have on quantum many-body systems has proven to be tremendously challenging and has attracted interest in the context of quenched or collapsing Bose-Einstein condensates (BECs) [6][7][8][9][10][11][12][13][14], especially at unitary interactions. Here the rapid buildup of three-body and higher-order correlations renders the description a challenge whose solution remains elusive.
A different approach to try to understand the effect of Efimov-like correlations in a many-body setting is to study the problem of a mobile impurity immersed in a weakly interacting BEC, leading to the formation of a Bose polaron. Importantly, this problem exhibits the heteronuclear instead of the homonuclear Efimov effect and the relevant three-body processes always involve the impurity. As a result, the most important three-body correlations are localized around the impurity, making the theoretical description more feasible. In particular, this allows the use of variational approaches [15][16][17], which is much more difficult when higher-order correlations throughout the BEC have to be accounted for. Also experimentally, probing the Bose polaron can give new insight, due to the available technique of radiofrequency (rf) spectroscopy for its observation [18][19][20].
The concept of the Bose polaron as an impurity immersed in a bath of bosonic excitations was introduced by Landau [21] to describe an electron dressed by phonons in a solid. A model that is commonly used in condensed matter to describe this scenario is the Fröhlich model [22]. However, the Fröhlich Hamiltonian does not contain all necessary terms to describe the cold atom Bose polaron [23], since it does not allow for bound state formation. Therefore, to describe processes for which bound state physics is important, such as the Efimov effect, the Fröhlich model needs to be extended. Furthermore, the description of the wavefunction needs to include the relevant interboson correlations. For this reason, the Bose polaron problem in the cold-atom context has proven to be challenging. Where the Fermi polaron problem can 0 1/(a )  of the polaron found using Gaussian States (blue) and mean-field theory (red dash-dotted) as a function of the inverse scattering length (1/aΛ). In grey the wave number of the Efimov trimer (solid line), and larger Efimov clusters (dashed lines) is shown. Note that the dashed lines are only schematic and not all possible Efimov clusters are displayed. For small scattering lengths (the green area), the polaron is the ground state of the system and is thus a stable quasiparticle. For |a| > |a (min) − |, the smallest scattering length at which an Efimov cluster can be formed, the polaron turns into a metastable excited state, which can decay into the large Efimov clusters via many-body scattering processes. These scattering processes become more likely as the absolute value of the scattering length increases, up to the critical scattering length a * , at which the polaron is no longer well-defined. At this point scattering of the polaron with one or two more bath particles can lead to decay into an Efimov cluster. In contrast, in the simpler mean-field picture, the polaron remains stable up to the scattering length a0, at which an infinite number of particles pile up in two-body bound states formed with the impurity.
be understood with relatively simple theoretical models [24][25][26], various descriptions of the strong-coupling Bose polaron have yielded widely varying results [15][16][17][27][28][29][30]. In Refs. [15,17] variational Ansätze with a small number of excitations were used, and a smooth crossover between the polaron and an Efimov cluster was predicted, corroborated by Quantum Monte Carlo results in Ref. [27]. However, when using a coherent state variational Ansatz which is capable of truly describing a deformation of the BEC instead of only excitations on top of the homogeneous BEC, one finds that the polaron becomes unstable due to an infinite number of excitations piling up on the impurity [16]. This instability is prevented by explicitly taking into account interboson repulsion [29][30][31][32][33], but still the polaron cloud will contain a large number of particles. What the coherent state approach clearly lacks, however, are interboson correlations, meaning the Efimov effect cannot be described. Also, while renormalization group theory [28] predicts the polaron to be unstable already for attractive scattering lengths, no clear connection to Efimov physics is evident.
In this work we systematically extend the coherent state Ansatz by fully including two-body interboson correlations, and three-body correlations including the impurity. To achieve this, we use Gaussian States [34] in the reference frame of the impurity. As a result, we can both describe the many-body physics of the creation of a macroscopic polaron cloud and the few-body physics of the Efimov effect. In particular, we study how these two phenomena affect each other. We choose to focus on light impurities because in this case the Efimov effect is much more prominent than in the case of equal mass [3,[35][36][37][38][39], since a light impurity can more efficiently mediate interactions between the bosons compared to a heavy impurity.
We find that the polaron is no longer the ground state of the extended Fröhlich Hamiltonian when three-body correlations are taken into account, due to the presence of large Efimov clusters much lower in the energy spectrum. However, the polaron remains a metastable excited state up to some critical scattering length a * , at which the energy barrier protecting the polaron from decaying into an Efimov cluster disappears. This critical scattering length can be interpreted as a many-body shifted Efimov resonance.
To explain our results, we first need to understand how the heteronuclear Efimov effect behaves for more than three particles. For the homonuclear case this is well understood [3,40] and four-and five particle Efimov clusters have been experimentally observed [41,42]. In this case larger and larger clusters are more and more tightly bound. This can be understood from simple dimensional arguments [3]. When the atoms have pairwise interactions the binding energy grows with N 2 , whereas the kinetic energy only grows with N . However, the properties of larger Efimov clusters are not as well explored for the heteronuclear case. When the bosons are not interacting or interacting repulsively, and only interact attractively with the impurity, this means that now also the binding energy scales with N . This leads to a non-trivial competition between the kinetic and potential energy, and the breakdown of the simple dimensional argument. As a result, providing a general picture of the binding energy of clusters with increasing size, is more difficult for the heteronuclear case. We show here that this effect of larger clusters being more tightly bound still persists for noninteracting bosons when the boson-impurity interaction is modeled with a single-channel model. We will refer to this effect as a "cooperative binding" effect. This term is imported from chemistry, where it refers to the effect that the particles cooperate with each other to lead to ever stronger binding. This occurs for example in the binding of oxygen molecules to hemoglobin.
Strikingly, we find that the existence of this cooperative binding effect has a profound impact on the stability of the polaron. We have illustrated this in Fig. 1. In this figure we schematically plot the wavenumber (which is proportional to the square root of the energy) of the polaronic state found from Gaussian (blue solid line) and coherent states (red dash-dotted line). The result from coherent states we refer to as mean field theory. The grey solid line indicates the position of the lowest threebody Efimov state in vacuum, appearing at a − . The grey dashed lines indicate the lowest Efimov clusters of increasing particle number, which shift to the left as a function of particle number. We prove using simple arguments that there is a minimum scattering length a (min) − , needed to form any Efimov cluster which is also indicated in the figure.
The background color of the figure indicates the stability of the polaron. In the green area, no Efimov cluster can be formed yet, meaning that the polaron is the ground state of the extended Fröhlich Hamiltonian. When crossing a (min) − into the blue area, the polaron will no longer be the absolute ground state, but a metastable excited state. For scattering lengths |a| < |a * |, the decay of the polaron into Efimov clusters, such as indicated with the wiggly vertical arrows, is extremely slow and needs simultaneous scattering of a large number of bosons on the polaron. However, as the scattering length increases, both the number of particles in the polaron cloud increases and the number of particles needed to form a bound state decreases. For both of these reasons the decay processes will become more and more important. When a * is crossed into the red area, scattering of one or two additional bosons is sufficient to cause the breakdown of the polaron into an Efimov cluster. Since these decay processes are included in our Gaussian state Ansatz, the local minimum of the polaron on our variational manifold disappears, meaning it is no longer a metastable quasiparticle. Since the BEC provides a coupling between the different particle number clusters, we find a cascade of the wavefunction into ever larger Efimov clusters.
Remarkably, |a * | is smaller than |a − |, because the polaron can decay immediately into larger clusters than the Efimov trimer, needing smaller scattering lengths due to the cooperative binding. In experiments, the formation of these large Efimov clusters leads to the rapid chemical recombination into more deeply bound states not included into our model. Furthermore, the shift of the Efimov resonance, means that the experimentally observed resonant recombination signal will also be shifted. This provides a fascinating example of chemistry in a quantum medium and polaron could-enhanced reactivity [1].
The energy of the polaron calculated with Gaussian states closely follows the energy of the mean field polaron, up to the scattering length a * . If no three-body correlations are included, the polaron remains stable even across unitarity up to the scattering length a 0 . Here the mean field polaron breaks down due to the piling up of an infinite number of bosons in a two-body bound state with the impurity. The Gaussian state polaron breaks down due to decay into an Efimov cluster.
The structure of the main body of our paper is the following: after showing our theoretical methods in section II, we prove in section III that our Gaussian state Ansatz incorporates Efimov physics. We present evidence for the cooperative binding effect, extending to large particle numbers, and we elucidate its mechanism. We demonstrate that many-body Efimov clusters form the ground state of the Hamiltonian. Then, in section IV, we use the Gaussian sate Ansatz to calculate the energy of the polaron as a function of scattering length and density and demonstrate its abrupt instability. Finally, we will discuss in detail how our results can be probed experimentally.

A. The Hamiltonian
To describe the impurity immersed in a three dimensional, infinite, weakly-interacting BEC at zero temperature, we use the extended Fröhlich model introduced in Ref. [23].
Here we treat the impurity with mass M in first quantization, with position and momentum operatorsR andP . The bosons of mass m are described using second quantization, where operatorsâ † k andâ k create/annihilate bosons of momentum k that interact via a contact interaction with scattering length a B determined by a coupling strength g B . The boson-impurity interaction is modeled using a regularized contact interaction and coupling strength g, with corresponding scattering length a. In Eq. (1), dk is shorthand for 3 . We have used a single-channel model to describe the interactions. This is supposed to be a good description of the interactions in case in the vicinity of a broad Feshbach resonance.
We model the bosons within the Bogoliubov approximation and introduce the quasiparticle operatorsb † k and b k . The validity of this approximation will be discussed in detail in the following sections. After the Bogoliubov rotation, we apply the unitary Lee-Low-Pines transformation [43]Û LLP to move to the frame of the impurity, Crucially, this transformation reduces the number of degrees of freedom of the problem, by leading to a transformed HamiltonianĤ =Û † LLPĤ 0ÛLLP , namely: where n 0 is the density of the BEC. The definitions of the quasiparticle dispersion ω k , and the variables, k,k and V (2) k,k are given in the appendix A. Since the impurity momentum operatorP now commutes with the Hamiltonian, it has been replaced by a vector P , which is the total momentum of the system. We set P = 0. Normal ordering the remaining quartic term gives rise to a modified dispersion of the quasiparticles: ω k + k 2 2M . To regularize the boson-impurity contact interaction, we introduce a UV cutoff Λ on all momentum integrals, leading to: where µ r = mM m+M denotes the reduced mass. The value of Λ is related to the effective range of the potential [44]: Λ −1 is proportional to the Van der Waals length l vdw . It therefore also determines the so-called three-body parameter, which is crucial for Efimov physics, since it sets the position of the first Efimov resonance a − . The value of Λ in our model should be chosen so that the obtained value of a − matches the experimental value of the system of interest.
To reduce computational complexity, we make use of the spherical symmetry of the problem and transform from plane wave modes to spherical waves: The Hamiltonian then takes the form: where we have defined: k = Λ 0 dk. The second term, which we will denote byĤ QLLP , is written without the transformation to spherical waves. The full expression after the transformation is very lengthy, and given in Appendix B.
This quartic term will later be shown to be extremely important to describe the Efimov effect. Indeed it is the only quartic term of the Hamiltonian, thus it is fully responsible for describing the effective interboson interactions mediated by the impurity. Because this term contains the vector product k 1 · k 2 , it contains the information about the relative direction of the bosonic momentum, allowing the bosons to distribute kinetic energy more efficiently. This manifests itself in a coupling between consecutive angular momentum terms.
Note that even though we have set P = 0, this does not mean that the impurity is stationary. Even though P = 0, P 2 is not zero in the lab frame, and the kinetic energy of the impurity is crucially important for the results.

B. Gaussian states
Throughout this work we use a variational method based on a Gaussian state variational Ansatz. Using Gaussian states to describe the Bose polaron was already suggested in Ref. [45], but in their description the beyond-Fröhlich parts of the Hamiltonian were missing. These terms are necessary for describing bound state formation and Efimov physics. Furthermore, the Gaussian part of the wave function was only treated perturbatively.
Here we use Gaussian states in the formalism discussed in the review paper by Shi et al. [34]. Note that since we use the Lee-Low-Pines transformation, this entangles the impurity degrees of freedom with bosons, effectively giving a beyond-Gaussian Ansatz in the lab frame. The notation we use is most similar to the presentation in Ref. [46] The variational Ansatz in the reference frame of the impurity has the form: In this equationΨ klm = (b klm ,b † klm ) T is the Nambu vector containing the bosonic creation and annihilation operators, Σ z = σ z δ(k − k )δ l,l δ m,m with Pauli matrix σ z , Φ klm = Ψ = (φ klm , φ * klm ) T is the coherent displacement vector, and ξ is the Hermitian correlation matrix. We do not denote the state of the impurity, because in the frame of the impurity this is trivial: the impurity is stationary at the center of the frame. Furthermore, our Ansatz acts on the state |BEC , which is the state of the background BEC.
The variational parameters we use here are the coherent displacement Φ and the covariance matrix Γ klm,k l m = {δΨ klm , δΨ † k lm } of the fluctuation field δΨ klm =Ψ klm − Φ klm . The covariance matrix Γ is related to ξ via the symplectic matrix S = exp(iΣ z ξ) as Γ = SS † . We use two complementary algorithms to optimize the variational parameters and to find the minima on our variational manifold: imaginary time evolution and iterated Bogoliubov theory. In principle these two algorithms are equivalent [47], but dependent on the situation either one is computationally more efficient. The equations of motion (EOM) for the imaginary time evolution are given by [34]: where we define The expectation value of the Hamiltonian can be calculated using Wick's theorem. With iterated Bogoliubov theory, one sequentially solves η = 0 for φ to update the displacement and diagonalizes H using a symplectic matrix S, to update Γ = SS † . Performing a single iteration of this algorithm is equivalent to solving the Gross-Pitaevskii equation (solving for φ) and calculating the equivalent of the Lee-Huang-Yang correction (diagonalizing H ) [46]. This is especially efficient compared to imaginary time evolution in case of very slow dynamics and when η is linear in φ.
As we will discuss in the following, for some of our calculations we will fix the expectation value of the particle number during the optimization. This works only for the imaginary time evolution, where it is achieved by adding a chemical potential term µ NN to the Hamiltonian, with particle number operatorN . Then the chemical potential µ N is dynamically adjusted to keep the particle number fixed. This can be done using [48]:

III. EFIMOV CLUSTERS AND COOPERATIVE BINDING
Cooperative binding is an important aspect of the Efimov effect, since a three-body Efimov state is bound when no two-body bound state can be formed yet. In this section we demonstrate that this cooperative binding can in fact persist for an arbitrary number of particles. We consider the scenario of an impurity interacting with bosons without interactions among themselves via a single-channel model. The cooperative character becomes evident from our finding that adding an extra particle to a bound state also lowers the energy of the other particles, or in other words, that a N + 1-body bound state is easier to form than an N -body bound state. In terms of the illustration in Fig. 1, this means that the grey dashed lines move further and further to the left for increasing N . These strongly bound Efimov clusters are lower in energy than the polaron and therefore the polaron is not the ground state but an excited state of the Hamiltonian.
A. The two-and three-body problem

Variational Ansätze
Before moving on to the description of large Efimov clusters, it is instructive to first consider the two-and three-body problem (i.e. with one, or two bosons plus the impurity). To achieve this we set the background density n 0 in Eq. (6) to zero, and we consider non-interacting bosons (a B = 0). This yields the Hamiltonian For the two-and three-body problem, we use the exact variational Ansätze In these expressions angular momentum conservation has been taken into account and the total angular momentum has been set to zero. For the two-body problem (Eq. (15)), this means that the single boson always has angular momentum zero, whereas in the three-body problem (Eq. (16)), the two bosons always have opposite angular momenta in the frame of the impurity. Moreover the prefactor only depends on m with the sign (−1) m , because of the spherical symmetry.
Despite the complicated form of H QLLP , the EOM for real-time evolution are relatively simple when these symmetries are implemented, and one finds

The non-interacting problem
First, we consider the non-interacting problem (g = 0). In the non-interacting problem the momenta of each particle are good quantum numbers. As a result, the two-body problem is trivial, since the total momentum is fixed to zero, and the impurity and the boson will thus always have opposite momenta, giving an energy of k 2 2µr . TheĤ QLLP -term does not play any role here.
In contrast, for the three-body problem, the EOM in the impurity frame (Eq. (18)) are non-trivial, since thê H QLLP term couples the different angular momentum modes. The lowest energy for fixed momenta of the two bosons k 1 and k 2 and total momentum P = 0 occurs when the two particles move in opposite directions. To describe this, the directions of movement of the bosons need to be correlated, which is achieved by the coupling between the angular momentum modes.
The non-interacting EOM of the three-body problem can be mapped to the problem of two coupled harmonic oscillators. This can exactly be solved using a Bogoliubov rotation. In this case an energy is obtained of

2M
, corresponding to an impurity momentum of k 1 − k 2 . This means that the bosons indeed move in opposite directions. In the corresponding wave function, every angular momentum mode is equally populated. Remarkably, this implies that an infinite number of modes is needed to solve the non-interacting problem in this framework. When g < 0, as considered below, lower angular momentum modes are favored, since only the zeroth angular momentum mode has an (attractive) coupling with the impurity. In case of bound-state formation, this means that convergence can be achieved with a smaller number of modes.

Including interactions
Now we consider the case of a finite, negative g. In the two-body problem, a bound state is formed for positive scattering lengths, which vanishes into the scattering continuum at infinite scattering length. However, for the three particles, due to the Efimov effect, a bound state can already be formed for finite negative scattering lengths. Indeed, already at this point cooperative binding is demonstrated, which is driven by theĤ QLLPterm of the Hamiltonian. WhenĤ QLLP is not included (leading to the second line of Eq. (18)), the EOM become separable, and the EOM of the two-body problem are retrieved. This shows that theĤ QLLP -term must be crucial for the Efimov effect. Since this term originates from the kinetic energy of the impurity, this shows that reducing the kinetic energy of the impurity is the driving force of the Efimov effect.
The full three-body problem can be solved numerically using sparse matrix diagonalization employing linear or, more efficiently, logarithmic grids in k. This way, the ground state and low-lying excited Efimov states can be found efficiently. In principle this method can be generalized to arbitrary interaction potentials between the bosons and the impurity. In that case the matrix elements of the potential may need to be precomputed numerically, since they typically do not have a straightforward representation in spherical wave momentum space. Including boson-boson interactions in this framework leads to complicated expressions for the relevant matrix elements and removes the sparsity of the matrices that represent the equations of motion, greatly increasing numerical complexity.
In Table I we show the Efimov scattering length a − following from our model for 6 Li-atoms interacting with several types of bosons. We see that a − is smallest for light impurities, meaning that the cooperative binding is the strongest in that case. This can easily be understood from its mechanism: reduction of the kinetic energy of the impurity. Because the kinetic energy of the impurity for a fixed momentum is larger for light impurities this implies that reducing this energy can have a larger impact. This can also be understood directly from thê H QLLP -term of the Hamiltonian, which scales inversely with the mass of the impurity. Since a − quickly grows as the mass ratio decreases, observing Efimov physics experimentally can most easily be achieved by immersing a light impurity in a bath of heavy bosons. In this paper, we will use the mass ratio of the experimentally available system of a 6 Li impurity immersed in a BEC of 133 Cs [37,38] as an example.

B. Many-body Efimov problem
To show that the cooperative binding effect of the Efimov effect also persists for larger particle numbers, one needs to demonstrate that the binding energy per boson E/N monotonically increases with N .
In order to study N -body bound states, the methods used in the previous subsection need to be extended to higher particle numbers. However, this route of exactly solving the N -body problem rapidly becomes computationally intractable. Instead, we adopt a different approach, based on a variational wave function that is not a particle number eigenstate, but rather a superposition of states with different particle numbers: a Gaussian state. In this Ansatz pairwise correlations between the bosons are fully taken into account in the frame of the impurity, which translates to three-body correlations between two bosons and the impurity in the lab frame. Instead of fixing the number of particles, we now fix the expectation value of the particle number N . We can show that E/ N grows monotonically with N , which is a strong indication that also E/N grows with N.
Note that the functional form of the correlations in a Gaussian state, is different from the correlation functions in the Jastrow-formalism [49] often used in condensedmatter physics, which was for example also used to describe the polaron in Ref. [27,50]. In our case, the covariance matrix, which characterizes the correlations, is optimized variationally without restrictions on the functional form of the correlations.
In the following we will now highlight our various findings.

Gaussian states incorporate 2-and 3-body physics
First, we demonstrate that a Gaussian state after the Lee-Low-Pines transformation is capable of describing Efimov physics. To this end, we show that the EOM of the two-and three-body problem can be retrieved exactly by linearizing the EOM for the Gaussian states around the vacuum. The procedure of linearizing the EOM is described in detail in Ref. [46] and we retrace here the main steps. Without loss of generality, ξ can be written as The equations of motion (EOM) for real time evolution are given by [46] The linearized EOM of Γ can be rewritten in terms of δξ In this equation, D is obtained by symplectic diagonalization of H (φ 0 , Γ 0 ) by a matrix S 0 with Γ 0 = S 0 S † 0 . The subscript 12 indicates the off-diagonal block in the Nambu basis. Writing we can diagonalize our equations around the vacuum to find Expressions for η, E and ∆ are given in the appendix, which are functions of the quantities: .
Note that despite the fact that the right-hand side of these expressions contains the angular momentum label m, the matrix elements of G and F formulated this way do not depend on m.
Finally, expanding to first order in δφ and δξ we find: Since δF = iδξ when linearizing around the vacuum, it is evident that these equations exactly correspond to the EOM for the two-and three-body problem: Eq. (17) and 18.

Cooperative binding with Gaussian states
To find evidence for the cooperative binding effect, we optimize the ground state energy of the system by imaginary time evolution. In this procedure we include a dynamically changing chemical potential to fix the average number of particles (not the density), given by: During the optimization with fixed particle number, the coherent part of the wave function goes to zero, leaving a purely Gaussian State. This is not surprising, since this allows for the maximum amount of correlation between the particles.
In Fig. 2 we show a contour plot of the energy per particle as a function of N and a. The bold line indicates the scattering length where the energy per particle becomes negative: the critical scattering length a c for bound state formation. Note that since the Hamiltonian conserves particle number, a c is a smooth line only due to a classical average over different particle number sectors. In the limit of N going to zero, a c → a − , the Efimov scattering length at which a 3-body bound state appears. This is expected, since we showed previously that in the vacuum limit the EOM for the 3-body problem are exactly recovered. As N is increased, |a c | monotonically decreases. Importantly, as soon as N is sufficient to lead to bound state formation, E/ N monotonically increases with N at fixed a, clearly demonstrating cooperative binding. This result seems in good correspondence with recent work by Blume [51], who also shows a strongly increasing energy as a function of particle number.
It is important to realize that our method is variational not for fixed N but for fixed N . This means that the energy we find for a given N = N may be lower than the lowest particle number eigenstate of N particles. This is not because of a coupling between the particle number eigenstates in the Hamiltonian, but because E( N ) is a concave function. This means that having a larger spread in particle number gives rise to a lower energy. Furthermore, the particle number fluctuations do not become relatively small for large N , such as often the case in statistical physics or for coherent states. Indeed, the standard deviation of the particle number for the Gaussian states scales with N . This standard deviation is actually maximized, since including a larger spread lowers the energy. The Gaussian form of the wave function does restrict the particle number statistics, meaning not arbitrary superpositions of particle number eigenstates are allowed. These features of our variational results with Gaussian states may be surprising, but do not invalidate our strong evidence for the cooperative binding effect.
One obvious question is whether a c will go to zero if N goes to infinity. A hint of how to approach this question is provided by the insight that the cooperative binding effect is driven by the quartic LLP-term originating from the kinetic energy of the impurity. A lower bound on |a c | can therefore be found by setting the kinetic energy of the impurity to zero. In this case the Hamiltonian becomes: (30) Formally, this corresponds to the Hamiltonian of bosons interacting with a infinitely heavy impurity. There is however an important subtlety: the relation between the scattering length and the coupling constant is determined by the mass of the impurity as described in Eq. (4). In our effective model, we therefore get a different scattering length a ef f corresponding to an impurity with infinite mass, which is related to g via In the spectrum of the Hamiltonian of Eq. 30 a bound state crosses the continuum when a ef f = −∞. Inserting this into Eq. 31 immediately yields the effective scattering length For Li-Cs this value is given by a c,lim = −0.07Λ −1 . As can be seen in Fig. 2 ≤ a c,lim . We stress that this limit is fundamental and independent of the variational method. With more advanced Ansätze the limit can only be approached faster or more closely.
In any case, from comparison with Table I it is evident that a c,lim is approximately two orders of magnitude smaller than a − . Note that such a minimal required scattering length to form an Efimov cluster does not exist for the homonuclear case, since there the pairwise interactions between all of the participating particles drive the cluster formation.
The increase in binding energy as a function of particle number and scattering length shown in Fig. 2   a decrease in the spatial extent of the Efimov clusters, which we define as As shown in Fig. 3, the size of the Efimov clusters monotonically decreases as a function of a and N . Note that R parametrizes the average distance of the bosons to the impurity, and not the distance of the particle furthest away. This means that the total extent of the Efimov cluster wave function can be much larger than shown here, especially close to a c . It is also obvious that these cooperative Efimov clusters are fairly tightly bound, with RΛ being of order 1.
Finally we corroborate our qualitative explanations of the cooperative binding mechanism with quantitative numerical results. In Fig. 4 we plot the relative values of the different energetic contributions: the kinetic energy of the bosons, the kinetic energy of the impurity, and the interaction energy, as a function of the scattering length and the particle number N . In this figure we clearly see that the reduction of the kinetic energy of the impurity compared to the interaction energy, is the driving force of the cooperative binding mechanism. Where the ratio between the interaction energy and the bosonic kinetic energy only changes marginally when the number of particles is increased, the kinetic energy of the impurity relative to the interaction energy can decrease with more than a factor two. Note that the kinetic energy of the impurity does not decrease in an absolute sense when more bosons are added. Adding more bosons leads to tighter binding and thus larger kinetic zero-point energy.

Further aspects
We have not discussed the structure of the manyparticle Efimov clusters in detail. One of the reasons is that these clusters are inherently unstable due to recombination into deeply bound molecular states, not included in our model. Because of this rapid recombination, the specific structure of these clusters is not relevant for experiments. Our main interest in this work is to show that the polaron, as discussed in the next section, can decay into these clusters, which is independent of their microscopic structure. For that physics, the cooperative binding effect is essential, making this our main focus. The results we obtain in Figs. 2, 3 and 4 are not universal since for N 1 the bound states that are formed are very tightly bound (as shown in Fig. 3). Instead, we predict a breakdown of universality due to the piling up of bosons around the impurity. The cooperative binding mechanism leading to this breakdown is, in turn, universal though, since the drivingĤ QLLP -term is independent of the interaction potential.
The accumulation of bosons on the impurity can be prevented or limited by including an interboson repulsion. Due to the tightly-bound nature of the bound states not only the interboson scattering length but also the range of the interboson interactions will play an important role, meaning a simple contact interaction will not be sufficient to describe this effect and more realistic potentials will need to be employed. Whether ultimately Van der Waals universality remains even for deeply bound states and states with more particles, is an interesting open question.
Furthermore, also using a two-channel model, where the interaction between the impurity and the bosons is modeled by a coupling to a closed molecular channel, may modify the results. Especially in the closed-channel dominated limit only a single boson can interact with the impurity at a time, resulting in an effective three-body repulsion [17]. Finally note that the coherent part of the Gaussian state may also contribute to cooperative binding, but more weakly than the Gaussian part. This is due to the mixing between the coherent term and the Gaussian term in the expectation value ofĤ QLLP . For the next section this is important, because only the coherent state term couples to the linear Fröhlich term in the Hamiltonian, responsible for the polaron formation.

IV. THE EFIMOV EFFECT AND THE BOSE POLARON
We have shown that large and strongly bound Efimov clusters already form at small negative scattering lengths. The ground state in our system is therefore not the Bose polaron, but a state where all particles are bound to the impurity. Questions concerning the fate of the polaron immediately arise. Does a polaron state still exist and if yes, will it still be stable? Can it decay into the cooperative Efimov clusters? Finding answers to these questions is the topic of this section.

A. Stable and metastable polaron
We are now in the position to discuss in more detail the spectrum shown in Fig. 1. Starting at the left of this figure, as discussed below Eq. (32) in the previous section, there is a minimum scattering length a can be very small for light impurities, yielding only a small region where the polaron is truly stable.
The first clusters that can form at a (min) − contain an exceedingly large number of particles. Since smaller clusters are not stable yet, this implies that all particles required for the bound state formation need to come together without the possibility to cascade through longlived intermediate states. Since cold atomic gases are dilute, the rate of the required N -body scattering processes is extremely small. Therefore the polaron state, which only leads to a relatively small and long-ranged deformation of the condensate at small scattering lengths, will not experience rapid decay into these clusters and remain metastable slightly beyond a (min) − . One can also explain mathematically from our Hamiltonian why the polaron should remain a metastable state. For small scattering lengths and a small number of excitations forming the polaron cloud, the contribution of theĤ QLLP term of the Hamiltonian is negligible. Since we are looking at attractive scattering lengths, no twobody bound states can be formed yet. This implies that in Eq. (3) the quadratic kinetic term dominates over the quadratic interaction terms. The model that therefore remains resembles closely the original Fröhlich model with a linear term driving polaron formation and a quadratic term counteracting this. This model must have a minimum of the energy as a function of particle number. This is illustrated graphically in Fig. 5, where the energy landscape of the Hamiltonian is plotted as a function of the number of excitations around the impurity for fixed density. The three subfigures correspond from left to right to the stable, metastable and unstable regime.
For 0 > a > a (min) − , the polaron is the global minimum of our energy landscape and we are in the regime illustrated in Fig. 5a). When a (min) − is crossed this does not immediately affect the local minimum since the barrier protecting it is still very large. Only for a very large number of excitations a bound state can be found. Math-ematically, as long as the particle number at the local minimum is small enough that the quartic term does not play a big role yet, the assumptions discussed before remain valid. This changes, however, for large scattering lengths or densities. In this case the particle number at the polaron minimum found from this simple model, is so large that the quartic term can no longer be neglected. When this occurs, the barrier protecting the polaron from decaying into Efimov clusters starts to disappear, as depicted in 5b). Eventually, at critical scattering length a * , this barrier completely disappears leading to a breakdown of the polaron, entering the regime of 5c).

B. The breakdown of the polaron
Now we will discuss in more detail the process leading to the breakdown of the polaron. At the point where the polaronic minimum disappears, the polaron state is no longer stable against a specific type of perturbations. These perturbations can be found from our model by linearizing the EOM such as explained in section III B 1. The resulting processes then represent single and double excitations on top of the Gaussian State. Thus, the polaron becomes unstable when a single or double excitation (or a linear combination thereof) can result in the formation of a bound state. Once a bound state is formed, the cooperative binding effect directly implies that at that scattering length also larger bound states can be formed with an ever increasing binding energy. This means that more and more bosons will accumulate on the impurity, leading to the eventual collapse of the system onto the impurity.
The value of the scattering length a * where the polaronic instability occurs, can be calculated numerically as a function of the density. To this end, we use the following procedure: we start at small scattering lengths where the local minimum corresponding to the polaron can be easily found. We then incrementally increase the coupling strength. In every step the polaronic properties are calculated with iterated Bogoliubov theory, using the polaron wave function at the previous scattering length as a starting point. As the critical scattering length a * is reached, the barrier preventing the polaron from decaying into many-body bound states disappears, leading to an abrupt divergence of the energy and particle number.

Polaronic instability and shift of the Efimov resonance
In Fig. 6 we use solid lines to plot a * as a function of the interparticle distance n −1/3 0 Λ for two values of the interboson scattering length a B Λ = 0.01 and a B Λ = 0.1. As a colormap, we additionally show the spectral weight Z of the polaron. The weight Z is given by the overlap of the wave function with the vacuum state, and it can experimentally be measured using injection rf-spectroscopy [18,19]. For Gaussian states, Z is given by: where I is the identity matrix. The red region close to the solid line in Fig. 6a), defines a region of dynamical instability, to be discussed further below. First we consider the low-density limit, where we find that a * → a − , i.e., the polaron ceases to exist exactly at the three-body Efimov scattering length. This finding can be understood as follows. For small densities, the polaron cloud is extremely dilute, meaning the impurity is practically free. A bound state for a free impurity plus two excitations from the background can be formed when the first three-body bound state crosses the continuum: at a − . This transition from a free impurity to a three-body Efimov state gives rise to a sharp drop of the quasiparticle weight from 1 to 0.
As the density is increased, a polaron cloud surrounding the impurity forms which contains excitations from the BEC. There are multiple consequences of this polaron cloud formation. First, as can be seen from 6, it leads to a reduction of the quasiparticle weight for increasing density. Second, due to the increased density around the impurity, scattering on the polaron can immediately lead to the formation of clusters of more than three particles. Importantly, due to the cooperative binding effect these larger clusters can already be formed at smaller scattering lengths than a − . This means that, as shown in Fig.  6, |a * | shifts to smaller scattering lengths as a function of density. This gives rise to the interpretation of the polaronic instability as a many-body shifted Efimov resonance, where the polaron now takes over the role of the free impurity as a collision partner of two additional bosons. This reasoning also implies that the time-scale associated with this instability is on the same order as the time-scale associated with three-body recombination.
One fascinating aspect of this finding, is that the shift of a * is continuous. While the average over particle number sectors in our discussion of cooperative binding in section III and Fig. 2 was purely classical, here the Hamiltonian coherently couples the different particle number sectors. This means that instead of obtaining a classical average over three-, four-and five-body Efimov states, the clusters that are formed in a BEC contain a quantum mechanical superposition of these different particle states. This originates from the quantum mechanical nature of the BEC, giving rise to the linear Fröhlich term in the Hamiltonian. This highlights an intriguing aspect of chemistry in a medium of a quantum nature.

Mechanism of the polaronic instability
We now study the mechanism leading to the polaronic instability in more detail. To interpret our results it is important to introduce the relevant length scales of the problem. These are: • the size of the Efimov trimer, which is on the order of a − .
• the average interparticle distance, parametrized by n • the size of the polaron cloud, which is determined by the modified healing length of the BEC: ξ B = (8πn 0 a B µ r /m) −1/2 [29].
In Fig. 7 we plot the number of excitations in the polaron cloud computed with a Gaussian (coherent) state Ansatz as a function of the scattering length as dashed (dotted) lines. Along with this, we show as solid lines the critical particle number needed to form a bound state with an energy lower than the polaron energy at the given scattering length and background density. In other words, we compare the number of particles contained in the polaronic excitation cloud with the number of particles that is needed to form a bound state. This figure is closely related to Fig. 5: the dashed and solid lines indicate the position of respectively the black circle and the black square in Fig. 5b).
The crosses at the end of the dashed lines indicate the scattering length a * of the polaronic instability. This instability is not captured by the coherent state Ansatz. In between the crosses and the diamonds that appear on some lines, we find regions of dynamical instability, discussed in more detail below. The black line in this figure shows the critical scattering length for bound state formation as a function of particle number in absence of a background BEC, such as shown in Fig. 2. The color indicates the background density. The difference between figure a) and b) is the interboson scattering length a B . First we compare the dotted lines and the dashed lines. We see that the dotted and dashed lines for the same density coincide for small a and N ex . However, for larger a, the particle number of the polaron calculated using Gaussian states grows much faster. Close to the critical scattering length this difference grows to 1-2 orders of magnitude. Now we also consider the solid lines. Naively one would expect that the onset of bound state formation would occur when the dashed lines cross the solid lines, but this appears not to be the case. In the low-density regime, the dashed lines in fact cross the solid lines, meaning that the polaron is stable, even when containing more than enough particles to form a bound state. The crucial insight to understand this behaviour, is that ξ B a − . Thus, even though the particle number is very large, the polaron cloud is so extended that the number of particles within the range a − is still too small to facilitate bound state formation. In contrast, in the high-density regime, the instability exactly occurs when the dashed lines hit the solid lines. In this regime the healing length that determines the size of the polaron cloud is smaller than, or of similar size as, the Efimov state. Therefore a bound state can be formed immediately when the polaron cloud contains the required number of particles. In the intermediate density regime, there is a crossover between these two regimes. Consistent with this interpretation, we find that the larger a B , and thus the smaller ξ B , the lower the density at which this transition occurs.
We see that |a * |, indicated with the crosses in Fig. 2, shifts towards smaller values for larger densities. From these graphs we find that the main mechanism of this shift is an increase of the number of excitations in the polaron cloud and their density as a function of background density (horizontal shift of dashed lines). Another contribution comes from the vertical shift of the solid lines. This is due to linear Fröhlich interaction term in the Hamiltonian, which leads to stabilization of the coherent polaron cloud due to the background density. This coherent part of the wavefunction also participates in the cooperative binding, leading to a downward shift of the solid lines. This is slightly counteracted by the interbosonic repulsion and the modified quasiparticle dispersion, which plays the strongest role for large densities and a B .
In our framework the most important role of a B is to determine the healing length. If we compare Figs. a) and b) we see that the healing length mainly determines the number of particles in the polaron. This is also reflected in the quasiparticle weight Z in Fig. 6 a) and b).

Properties of the metastable polaron
Above we have focused on the breakdown of the polaron described by Gaussian states. Another interesting question is how much the properties of the polaron are altered compared to the coherent state approach [16] by including interboson correlations in the regime where the polaron is metastable. In Fig. 8 we plot the energy difference ∆E = E − E M F with the mean field result: where a 0 is defined in appendix A. Surprisingly, the effect on the energy appears to be very small. We see that there is only a significant correction to the polaron energy close to the instability. Even there, the correction is smaller than 10 %. This leads to the conclusion that introducing correlations in the variational Ansatz leads to a decrease in spectral weight and destabilization at some point, but that the energy of the polaron is still very well described using mean-field theory. We have already seen by comparing the dashed and dotted lines in Fig. 7 that the number of excitations in the polaron cloud found from coherent or Gaussian states are very different. For small particle numbers and scattering lengths, they still coincide, since the quarticĤ QLLP -term has only a small contribution. When this term starts to be important, for scattering lengths approaching a * , the difference rapidly grows however. This also leads to a much smaller Z-factor for the Gaussian State result, as can be understood from Fig. 6, by comparing the dashed contours of Z M F compared to the colormap. This reflects the general notion that while variational energies may work well, the same may not apply for the wave functions.

Dynamical Instability
Our numerical results show a region of dynamical instability. A dynamical instability occurs when the variational parameters correspond to a minimum on the variational manifold, but to a saddle point of its tangent space. We identify the presence of the dynamical instability by linearizing the real time EOM around the minimum found by the iterated Bogoliubov theory. The system is stable when the symplectic diagonalization of the linearized time-evolution operator yields only positive real eigenvalues. However, in case of a dynamical instability, one finds imaginary eigenvalues, which corresponds to a negative direction in the Hessian. This manifests itself as an instability in the real-time evolution. The dynamical instability indicates that the variational manifold is no longer suitable to describe this metastable state of the system in this region. To solve this problem, one would need to incorporate even higher order correlations in the calculation. In our Figs. 6 the dynamical instability occurs in the small red region attached to the solid line at intermediate densities for a B Λ = 0.01. This is exactly the region where the polaron is very large in number of particles and extent, but where the density at the impurity is too small to lead to bound state formation. A dynamical instability can also be found when extending our plots to higher density. Here our results are no longer valid, because when n −1/3 0 Λ > 1, one can no longer speak of universal long range physics, and short range physics will dominate the behaviour of the system. This corresponds to interparticle distances comparable in scale to the length scale of the interaction potentials. In gases of cold atoms the typical densities are orders of magnitude away from this limit. A dynamic instability was also found for the Bose polaron with coherent states in Ref. [52]. Although the character and position of this instability is different from the one we find here, the origin may be related.

The dependence on the mass of the impurity
Finally, we study the mass-dependence of our results. To this end, we investigate a significantly different mass ratio and choose the Li-Na system as an example. In this case we fix the interboson scattering length to a B Λ = 0.1 and we plot a * as a function of the density. The result is shown in figure 9. We plot a * both alongside the quasiparticle weight and the energy correction to the mean-field result.
The most obvious difference is the change in the scale of the scattering length. From the few-body problem we know that a − becomes much larger when mass ratio m/M decreases (see Table I) and we find that the same holds in the many-body case. This also leads to a change in the scale of the density, which is not a surprise, since the number which matters is the dimensionless ratio between the interparticle distance and the size of the Efimov state, which scales with a − . Except for these changes in scale there are no striking physical differences between the Li-Na and Li-Cs systems studied here. One small quantitative difference is the difference in the polaron energy between the Gaussian state and the mean field theory results. This difference is slightly larger for the Li-Na mass ratio than for Li-Cs.

A. Limitations of our approach
Our current approach has some limitations. First, we have only studied very small interboson scattering lengths. In realistic experiments a B Λ will be of order unity. We have limited ourselves to small a B to preserve the validity of the Bogoliubov approximation. To ensure this, we tested that going beyond the Bogoliubov approximation with these interboson scattering lengths only leads to small deviations from our results (< 10% for a * ).
Going beyond the Bogoliubov approximation comes with the challenge of strongly increasing numerical complexity and the necessity and difficulty to properly regularize the interboson interactions. Studies dealing with these issues are left for future work. Nevertheless, we expect that for experimentally realistic a B and for light impurities our results still hold qualitatively. For heavy impurities, however, we expect that the cooperative binding will be strongly suppressed due to the interboson repulsion.
Second, since we have used imaginary time evolution, we could only study the properties of local energy minima on our variational manifold. For that reason, we could in particular not go beyond a * to study, e.g., higher-order Efimov states or the repulsive polaron. Furthermore, calculating spectral functions with a Gaussian state variatonal manifold is challenging for various reasons. Therefore the fate of the polaron branch in our model beyond the instability remains an important open direction of study. We consider it unlikely that the branch will completely disappear and we hypothesize that it will be broadened due to the decay into the Efimov clusters. We expect that the resulting width of the spectral line will be inversely proportional to the timescale of the decay. At this point this timescale would be similar to the time scale of three-body recombination.
Third, our approach considers up to three-body correlations between two bosons and the impurity. It is thus a natural question to ask how our results would generalize when including even higher order correlations. We expect that when up to N -body correlations are included, the value of a * would asymptotically connect to the resonance of the N -body Efimov cluster. Due to the coop-erative binding effect, its absolute value would in turn be smaller than the |a * | we find here. As discussed before, however, we do not expect the polaron branch in the spectrum to suddenly end at the point of instability, but we expect it to be broadened by the time scale associated to its decay. Since cold atomic gases are typically very dilute, N -body scattering in experiments is typically not very prominent for N > 3. As a result, we expect that the inclusion of beyond three-body correlations would in practice not have a strong effect on the observed results.
Finally our results may not fully apply to closedchannel dominated Feshbach resonances, because we use a single channel model. In the case of closed-channel dominated resonances there are effective repulsive threebody interactions counteracting the cooperative binding, which are not included in our model.

B. Comparison to other theoretical approaches
We have already compared our results quantitatively with the mean field result obtained from a coherent state approach [16,29]. The most drastic difference is that the coherent state approach does not include three-body correlations, which are needed to describe the Efimov effect, and that it therefore completely misses the cooperative binding effect and the presence of the Efimov clusters. In this model, the polaron would be the ground state for all negative scattering lengths and even beyond unitarity. The polaron energy will then diverge at the positive scattering length a 0 , due to a large number of excitations piling up on the impurity. The scattering length a 0 can be interpreted as a shifted version of the unitarity point, so the point where a two-body bound state crosses the continuum. In contrast to the Gaussian state model, where a * is shifted to smaller coupling strengths compared to a − , the shift from unitarity to a 0 is a shift to stronger coupling strengths. This highlights also the different mechanism of the shifts. For Gaussian states the shift is predominantly caused by the cooperative binding mechanism, whereas the interboson repulsion on the Bogoliubov level causes the shift in case of coherent states.
Interboson interactions play an important role in our model by setting the healing length of the BEC, and hence the extent of the polaron cloud. Guenther et al. [29] show that explicitly including interboson interactions can prevent the collapse of an infinite number of bosons onto the impurity. However, this only plays a role for large numbers of bosons and high densities in the polaron cloud. In principle, explicit interboson repulsion will also limit the number of particles that an Efimov cluster can host. However, for the light impurities and small scattering lengths we consider, the size of the Efimov cluster where this effect will play an important role is much larger than the Efimov clusters into which the polaron state initially decays. Therefore this interboson repulsion would not qualitatively change the critical scattering length at which the polaron instability occurs.
Using a renormalization group approach [28], it was found that the polaron breaks down for a finite negative scattering length. This was attributed to phase and particle number fluctuations and no connection to the Efimov effect was made. Since the general picture of a polaronic instability at negative scattering lengths matches with our results, it remains an interesting open question how the two pictures are related.
In Refs. [15,17,39] the interplay between Efimov physics and Bose polaron formation was studied for a limited number of excitations. For a mass ratio m/M = 1 a smooth crossover between the Bose polaron and the lowest energy Efimov cluster was predicted [15,17] and this was supported by quantum Monte Carlo results [27]. In this case the most deeply bound Efimov cluster contained a limited number of particles. In our case the lowest energy Efimov cluster contains an infinite particle number, and there can hence be no smooth crossover from the polaron into this state. These qualitative differences can be explained by the different mass ratios that were used and the effective three-body repulsion which is included in the two-channel model used in Refs. [15,17]. In Ref. [39], radiofrequency injection spectra were predicted for the same mass ratio we considered. However, since the number of possible excitations from the BEC in the employed diagrammatic approach was limited, the possibility of forming many-particle Efimov clusters was not included.

C. Experimental implementation
Our predictions have not yet been tested experimentally [18][19][20]. In the Bose polaron experiments that were performed so far, the mass ratios that were used led to a strong suppression of the Efimov effect. As a result, any Efimov features would only have been observable at very large negative or positive scattering lengths. In this work we focus on the example of 6 Li impurities in a BEC of 133 Cs, but for Li in Rb or K, as also available in experiments [53,54], the results would be very similar.
The typical densities of above BECs are in the range of 10 13 -10 15 cm −3 [55]. When we assume that the threebody parameter Λ ≈ l −1 vdw and we take the Van der Waals-length of Li-Cs of 45a 0 this gives a regime of n −1/3 0 Λ ∼ 40 − 200. This means that the low and intermediate density regime of our results can be probed, and a clear shift of the Efimov resonance should be observable. Tuning the density in the experiment for a given BEC may in practice be difficult. However, in any case the resonance position in a BEC can be compared with the resonance position in a thermal gas, which should give comparable results to our low density limit.
In the particular case of the Efimov resonances observed for Li-Cs at positive Cs-Cs scattering length [56,57], there is a subtlety which plays a role. In this case the resonance corresponding to the lowest energy Efimov is suppressed due to coupling to a shallow Cs 2 bound state. Therefore the first observed Efimov resonance corresponds to the second Efimov state and appears at a scattering length of around −2000a 0 . If we use this Efimov resonance as the lowest resonance in our model, we find a larger three-body parameter Λ −1 = 8 l vdw . Due to the larger size of the Efimov state, this means that now the high density regime of our results could be probed: n −1/3 0 Λ ∼ 5 − 25. Furthermore, when the mass ratio m/M is smaller such for 6 Li in a BEC of 23 Na, the value of a − is naturally larger. This implies that n −1/3 0 a − is larger and relatively higher densities are reached as well.
We propose combining two experimental approaches to test our theoretical predictions. The first approach would be to perform loss measurements such as typically used to observe Efimov resonances [4,37,38,41,42,56,57]. In this case the magnetic field should adiabatically be ramped to form a polaron. Then at a given final scattering length a, the magnetic field should be kept fixed to measure the loss arising from recombination. When now this final scattering length is varied, one should observe a clear enhancement of the loss as a * is reached. Whether this appears as a resonant feature or as the onset of a regime where the three-body recombination is enhanced, is still and open question and a subject of further study. This may also depend on the density. The loss measurements may be more efficiently performed by using the recently introduced photoassociative ionization technique [58].
The second approach would be to perform rf-injection spectroscopy such as in Refs. [18,19] for lighter impurities. We expect that a clear drop in quasiparticle weight and broadening of the spectral line of the polaron will be observed when the scattering length of the polaronic breakdown is approached and crossed. When doing ejection spectroscopy as in Ref. [20] and the ground state polaron is actually prepared, enhanced three-body loss due to the Efimov cluster formation will be the most important observable. The formation of tightly bound Efimov clusters may also give rise to higher frequency contributions in the rf-spectrum because of their large amount of kinetic energy.
In conclusion, observation of our theoretically predicted phenomena are in reach with current state-ofthe-art experimental techniques.

VI. CONCLUSION AND OUTLOOK
We use a Gaussian state variational Ansatz to describe the Bose polaron problem and the Efimov effect. We find that the cooperative binding caused by the Efimov effect leads to the formation of many-particle Efimov clusters. The cooperative binding is driven by the reduction of the kinetic energy of the impurity. Since these Efimov clusters are lower in energy than the Bose polaron, the polaron is not the ground state of the extended Fröhlich Hamiltonian but rather exists as a metastable excited state. This excited state loses its stability at a critical scattering length a * that can be interpreted as a manybody shifted Efimov resonance. We predict that while the mean-field energy of the polaron is reliable up to the point where the polaron becomes unstable, the inclusion of interboson correlations leads to a strong decrease in the spectral weight. Our results can be experimentally probed by a combination of rf-spectroscopy and threebody loss measurements of light impurities immersed in BECs. The parameter regimes discussed in our work are experimentally feasible, if suitable magnetic fields are found that feature both very small interboson scattering length and a large boson-impurity scattering length. We expect our results to also hold, up to a quantitative shift, for larger interboson scattering lengths.
Future interesting directions include the study of the real-time dynamics of the polaron [50,52,59] using Gaussian states, to understand how the Efimov cluster formation occurs in real-time and whether indeed resonant behaviour can be observed at the scattering length a * . Furthermore, the scope of our results can be extended to include finite total momentum, explicit interboson repulsion and by going beyond the single-channel model. This will respectively allow the study of the dispersion relation of the polaron, the repulsive side of Feshbach resonances and closed-channel dominated resonances. Our methods can also be extended to rotating impurities [60][61][62] or bipolarons [63][64][65] to study the effect of interboson correlations in those systems.
Finally, it would be fascinating to explore further the connection to quenched BECs [6][7][8][9][10][11][12][13][14]. In this work, we have shown that the formation of a polaron cloud around an impurity can lead to a modification of Efimov physics. One natural question to ask, is whether a similar effect occurs in BECs quenched to attractive scattering lengths. While one can certainly not employ the language of polarons or polaron clouds in this case, it is still pairwise correlations between the bosons that lead them to cluster closer together, which is also the essence of polaron physics. Hence, it will be interesting to explore whether a shift of Efimov resonances can also be observed in said scenarios. It may be possible to explore such effects by applying the cumulant expansion method as described in Ref. [13,14] to attractive scattering lengths.
Another question is how similar the polaronic instability is to the collapse of a BEC. In our model with bosonic interactions on the Bogoliubov level, the single impurity induces a first-order phase transition of the entire system at a (min) − . Even though this is prevented in practice by explicit interboson repulsion, it still indicates that the impurity can have a profound effect on the medium. Furthermore, the study of how large an Efimov cluster in a medium can become before it will lead to recombination and loss, remains another interesting open question. These questions highlight how the study of the impurity scenario can give new insights into the dynamics of quenched or collapsing BECs.