Self-stabilized Bose polarons

The mobile impurity in a Bose-Einstein condensate (BEC) is a paradigmatic many-body problem. For weak interaction between the impurity and the BEC, the impurity deforms the BEC only slightly and it is well described within the Fr\"ohlich model and the Bogoliubov approximation. For strong local attraction this standard approach, however, fails to balance the local attraction with the weak repulsion between the BEC particles and predicts an instability where an infinite number of bosons is attracted toward the impurity. Here we present a solution of the Bose polaron problem beyond the Bogoliubov approximation which includes the local repulsion between bosons and thereby stabilizes the Bose polaron even near and beyond the scattering resonance. We show that the Bose polaron energy remains bounded from below across the resonance and the size of the polaron dressing cloud stays finite. Our results demonstrate how the dressing cloud replaces the attractive impurity potential with an effective many-body potential that excludes binding. We find that at resonance, including the effects of boson repulsion, the polaron energy depends universally on the effective range. Moreover, while the impurity contact is strongly peaked at positive scattering length, it remains always finite. Our solution highlights how Bose polarons are self-stabilized by repulsion, providing a mechanism to understand quench dynamics and nonequilibrium time evolution at strong coupling.


I. INTRODUCTION
Impurities in a Bose-Einstein condensate (BEC) exhibit a multitude of fundamental physical phenomena: the formation of quasiparticles [1], Efimov bound states [2,3], synthetic Lamb shifts [4,5], Casimir interactions induced by a fluctuating medium [6,7], and quantum criticality [8]. Current experiments with ultracold atomic gases are investigating several of these effects reaching far into the strong-coupling regime [4,[8][9][10][11][12][13]. For understanding experimental observations it is thus vital to develop a theoretical model that applies to impurity systems at strong coupling and that can address both ground state and non-equilibrium phenomena. At weak interactions between the impurity and the BEC, the impurity deforms the BEC only slightly and physics is well described within the Fröhlich model in terms of long-wavelength phonon excitations [1,[14][15][16][17][18][19][20][21][22][23][24]. For strong local attraction, instead, the Fröhlich model is incomplete and needs to be amended by quadratic terms that absorb and re-emit phonons [25]. Importantly, these term are also required to correctly capture the formation of bound states between the impurity and bath atoms which is a crucial ingredient for the physics of polorans at strong coupling [26]. Variational wave functions based on a single phonon excitation [27,28] are able to describe the single occupation of such a bound state. However, the bosons that make up the BEC tend to bunch, and at strong coupling it is energetically favorable to occupy the bound state multiple times leading to a gain of several times the binding energy. This process, recently observed with Rydberg atoms immersed in a BEC [29], is described by a coherent state variational ansatz [25,26,30,31] that allows for an arbitrary number of excitations and strong local deformations of the condensate [32][33][34][35][36][37][38].
Generally, the application of the variational principle for the determination of the ground state is only viable if the Hamiltonian is bounded from below. Under this condition a stable solution can be found, and theoretical approaches should rely only on such approximations that preserve stability of the underlying problem 1 . For the strongly interacting Bose polaron, the resulting theoretical challenge can be understood in the simple toy model illustrated in Fig. 1. Here a static attractive impurity is represented by a local potential well of finite range around the impurity. This well acts as a microtrap within the BEC. At strong coupling beyond a scattering resonance [39], the potential well (red dashed line) is deep enough to admit a bound state with energy ε i = −ε B < 0 (red bar). In absence of boson repulsion, in the manybody ground state all bosons would occupy the bound state and the ground-state energy E 0 = −N ε B is indeed unbounded from below in the thermodynamic limit: the whole BEC is collapsed onto the impurity.
Naturally, a local boson repulsion counteracts this process by balancing the impurity attraction [40] and thus providing a lower bound to the ground-state energy. This mechanism can be understood in terms of a singlesite Bose-Hubbard model [41,42] with local repulsion U 2 n i (n i − 1) on the impurity site r i , that competes with a local attractive potential energy −ε B n i for occupation number n i . For U > 0 the ground state has a finite occupation n i 2ε B /U , which is nonperturbative in the strength of the interaction: the size of this polaron dressing cloud grows for weaker repulsion. It is crucial to adequately capture this repulsive effect in the theoretical description of Bose polarons.
Previous variational approaches have included the boson repulsion only at the level of the Bogoliubov approximation. Here the interaction between the Bogoliubov quasiparticles is neglected, and thus bosons in the bound state fail to generate the compensating pressure required to ensure the stability of the ground state. As a consequence this approach falsely predicts instead a dynamical instability in presence of boson repulsion toward infinite occupation of the bound state in the strong coupling regime [30,31,36]. This shows the crucial importance of including the boson repulsion beyond the Bogoliubov approximation.
In this work, we present a stable variational approach to the Bose polaron problem at strong coupling. Our approach applies to arbitrary dimension and impurityboson scattering lengths and it provides a basis for the study of dynamical properties of Bose polarons. In this way our work complements and extends previous approaches using quantum Monte Carlo to determine ground-state properties [43], studies of the role of boson repulsion in one-dimensional systems [44][45][46][47][48], or a recently developed nonlocal Gross-Pitaevskii (GP) theory for nonequilibrium dynamics [40]. Moreover, previous works raised the question to which extent polaron properties are universal in the short-range limit [49,50]. In the following we study this question across the Feshbach resonance, and in particular in the regime where a bound state is supported by the impurity-bath potential, extending previous work that considered how the properties of the Bose polaron depend on the range of the interactions, both for the finite-range boson repulsion in the nonlocal GP theory [40] and for finite-range impurity potentials [42,51,52].
Specifically, we study the effect of local boson repul-sion and finite-range attractive impurity potentials employing an inhomogeneous variational state that allows for large dressing clouds and strong local deformations of the BEC. In Sec. II we introduce the stable Bose polaron model and discuss its solution within Gross-Pitaevskii theory. In Section III we minimize the resulting GP energy functional and obtain the condensate wave function around the impurity. We find that even if the bare impurity potential admits a bound state, the emerging effective potential does not, thus providing a simple mechanism for the self-stabilization of Bose polarons. Section IV presents our results for the polaron energy, the size of the polaron dressing cloud and the Tan contact across the resonance. The universality of the polaron energy is discussed in Sec. V, and we show that the energy at unitarity depends universally on the effective range, as long as Efimov states can be neglected. Finally, in Sec. VI we compare variational approaches to the Bose polaron and discuss which Hamiltonians and energy functionals can provide rigorous bounds on the ground-state energy.

II. MODEL
We consider a single impurity particle immersed in an interacting Bose gas. The combined system is described by the Hamiltonian Here,p andx denote the momentum and position of the impurity of mass m I , andp B,i andx B,i characterize the bosons i = 1, . . . , N of mass m B . The boson interaction V BB (x) is assumed to be repulsive and the impurity-boson interaction V IB (x) is attractive.
The coupled system of impurity and bosons is conveniently analyzed in the reference frame comoving with the impurity. This is achieved by a canonical transformation introduced by Lee, Low, and Pines [53] and elaborated on by Girardeau [15]. The transformation S = exp(ix · ip B,i ) leads to the LLP Hamiltonian (2) in which the impurity operators have been eliminated, and p 0 denotes the conserved total momentum of the system. In the comoving frame, the impurity potential V IB (x) acts as a static external potential centred at the origin, while the kinetic term proportional to ∼ 1/m I accounts for the recoil of the impurity. For bosons in the vicinity of the impurity this term leads to induced interactions between bosonic particles in addition to their inter-boson repulsion V BB (x). Homogeneous Bogoliubov Approximation.-In order to appreciate the importance of the adequate inclusion of boson repulsion we briefly review approximations to the model (1) that are frequently applied to the study of the Bose polaron problem. In the formalism of second quantization Eq. (1) readŝ Here V is the system volume,d † p andâ † p are the creation operators of impurity and bosons, respectively, and V BB (q) and V IB (q) are the Fourier transforms of the respective interactions in Eq. (1). Next, the ladder operators of bosons are shifted using a simple canonical coherent state transformation leading toâ † p →â † p + δ p,0 √ N 0 . This is then followed by the crucial Bogoliubov approximation: All terms beyond quadratic order bosonic operators are neglected, leading to the truncated Hamiltonian where, for simplicity, we have chosen the example of a boson contact interaction of strength g BB . As a result of the Bogoliubov approximation the bosonic part of the model can be diagonalized using the standard Bogoliubov rotation While the Bogoliubov approximation thus allows to obtain a simple dispersion relation for bosonic quasiparticles, it, however, captures the effect of repulsion only within the Bogoliubov coefficients u p and v p . They give rise, e.g., to the modified quasiparticle dispersion ω p = p ( p + 2n 0 g BB ) where, notably, the non-deformed, homogenous boson density n 0 appears. It turns out that this approximate account of boson repulsion is insufficient to self-stabilize the Bose polaron. Indeed, the neglect of terms beyond quadratic order is responsible for the apparent instability of the truncated Hamiltonian (4) 2 . Instead, it is crucial to keep all terms beyond quadratic order which allows the repulsive interactions to act as a stabilizing counter term to the impurity attraction. Including these terms allows then to expand the theory not simply around the homogenous BEC but around a BEC that is already deformed due to the presence of the impurity (for a discussion of the one-dimensional case see Ref. [47]). As discussed in the following, this in turn allows one to effectively map the strong coupling Bose polaron problem onto a weakly interacting one.

A. Gross-Pitaevskii functional
Following this strategy we focus on a variational approach to the ground state of the full Hamiltonian (2). We use a product state [54] Ψ where the condensate wave function φ(x) is normalized to the condensate particle number d d x |φ(x)| 2 = N 0 3 . Both φ(x) and the ground-state energy are found by minimizing the resulting variational Gross-Pitaevskii (GP) energy functional Here µ is the chemical potential and we assume weak boson repulsion represented in three dimensions by a contact interaction of strength g = 4πa BB /m B with boson scattering length a BB > 0. Note that normal ordering of the impurity kinetic term in Eq. (2) contributes to the boson kinetic term with reduced mass m −1 The energy functional Eq. (7) exhibits two important limiting cases: (i) for an infinitely heavy impurity m I → ∞, the normal-ordered kinetic recoil term in the first line of (7) vanishes and the standard GP energy functional for bosons in a static external potential is recovered; (ii) for a Bose polaron at rest (p 0 = 0), and a radially symmetric impurity potential V IB (|x|), the wave function φ(x) is spherically symmetric and the recoil term again vanishes -it only re-appears beyond the product ansatz when boson correlations are included [40,55,56].
We find the condensate wave function φ(x) by minimizing the GP functional (7) in the thermodynamic limit subject to the boundary conditions |φ(x → 0)| < ∞, |φ(|x| → ∞)| = √ n 0 in terms of the condensate density n 0 far away from the impurity. For a radially symmetric impurity potential V IB (|x|) at rest (p 0 = 0), the ground-state wave function is spherically symmetric and real. In the solution of the GP functional the energy is universally expressed in units of the BEC bulk chemical potential µ = gn 0 and the distance from the impurity r = |x| is measured in units of the modified healing length ξ = 1/ √ 2m red µ = 1/ 8π(m red /m B )a BB n 0 , which involves the reduced mass and is therefore larger than the usual bulk healing length ξ 0 = 1/ √ 2m B µ of the BEC without impurity.
We define the polaron energy functional relative to the energy of the unperturbed BEC with wave function φ 0 = √ n 0 . For the threedimensional case it is conveniently expressed in terms of the scaled radial function u(r) = rφ(r)/ √ n 0 as Here the boundary conditions for φ(x) translate to u(0) = 0 and u(r → ∞) = r. Numerically, E[u] is minimized by global optimization on r ∈ [0, L] with L = 10 ξ and an r-grid spacing ∆r = 0.05 ξ much smaller than the potential range.
We thus obtain the scaling solution of the GP equation (GPE) in units of µ and ξ for a given dimensionless potential shapeṼ (x = r/ξ) = V IB (r)/µ. This scaling solution is universal for arbitrary values of the condensate density n 0 , boson scattering length a BB and mass ratio m red /m B as long as GP theory applies [40,42], with these parameters entering only indirectly via µ and ξ. For comparison with experiment, and in order to visualize the effect of boson repulsion, the universal solution can be rescaled to obtain the specific solution for desired values of the boson scattering length a BB and the mass ratio m red /m B in density units of energy E n = 2 n 2/3 0 /2m red and length n −1/3 0 .

B. Impurity potential
For the impurity potential V IB (r) we consider two different functional forms. This allows us to study the universality of the Bose polaron by analyzing how quasiparticle properties depend on the potential shape and range. Specifically, we consider an attractive Gaussian potential, and an exponentially decaying potential, both of depth V 0 > 0 and range R. The low-energy scattering properties of the impurity and boson are characterized by the impurity-boson scattering length a IB and the effective range r eff , which determine the leading terms of the effective range expansion of the momentum-dependent scattering phase shift in three dimensions. They are found by numerically solving the Schrödinger equation for the scattering of a boson with the impurity in the center-of-mass system via the potential V IB (r). Equivalently, one may solve the first-order nonlinear variable phase equation [57], which yields Here a(r) and r eff (r) obey the boundary conditions a(0) = 0, r eff (0) = 0, and account for the phase shift accumulated by the scattering wave function (for the generalization of the variable phase equation to singular potentials see [58]). Correspondingly, the differential equations (11) are integrated from r = 0 . . . ∞ which yields the scattering length a IB = a(r → ∞) and effective range r eff = r eff (r → ∞).

III. EFFECTIVE POTENTIAL
First, we present results for the condensate profile around the impurity that we find by minimizing the Gross-Pitaevskii energy functional (7)- (8).
For attractive impurity-boson interaction, the wave function φ(r) -which, in the comoving frame, directly yields the impurity-boson density-density correlation function g (2) IB (r) = |φ(r)| 2 /n 0 -is enhanced near the impurity, as shown in Fig. 2(a). In this figure we have chosen a potential V IB (r) (dashed line in Fig. 2(b)) that is sufficiently deep to support a two-body bound state at energy −ε B < 0. Correspondingly, the potential is characterized by a positive impurity-boson scattering length, in Fig. 2(b) chosen as a IB ≈ 4ξ > 0.
We thus realize a scenario as described in the introduction, which is mimicked by a microtrap or a single-site Bose-Hubbard model. While at the two-body level it suggests a dynamical instability where the occupation of the bound state would grow without bounds, many-body effects come to the rescue. Indeed, due to the repulsion between bosons, each additional boson trying to participate in the formation of the Bose polaron is subject to the effective potential that results from by the competition of the bare attractive potential V IB (r) (dashed line in Fig. 2(b)) and the repulsion created by the already existing polaron cloud.
IB (r) as determined from the condensate wave function φ(r) as function of the distance r from the impurity (blue solid). For an attractive impurity potential the wave function is enhanced near the impurity as compared to that of an unperturbed BEC φ0(r) = √ n0 (orange dashed). (b) The bare impurity potential (orange dashed) of Gaussian shape (V0/µ = 5.1651, R/ξ = 0.81892) has an effective range r eff = ξ and positive scattering length aIB = 4ξ, and correspondingly admits a bound state. An extra bosonic test particle is, however, subject to the effective potential (blue solid) that is weakened by the repulsion from the polaron cloud; while still attractive the effective potential is characterized by an effective, negative scattering length a IB,eff = −0.1ξ (and renormalizedr eff = −207ξ) which thus no longer supports a bound state. Hence, additional bosonic quantum fluctuations lead only to a weak, additional dressing of the impurity particle.
In essence, this effect can be understood as arising in an effective density-functional theory (DFT) for the BEC particles in the presence of the impurity: the bosons already attracted to the impurity screen the attractive potential and make it shallower, as shown by the blue line in Fig. 2(b). Moreover, we find that they create a small repulsive barrier at intermediate distance. As a consequence, the effective potential no longer admits a bound state and it is correspondingly characterized by a negative effective impurity-boson scattering length a IB,eff ≈ −0.1ξ. Hence, the single-particle excitation spectrum for each additional boson is bounded from below: the dynamical instability is replaced by Bose polarons self-stabilized by their dressing cloud.

IV. BOSE POLARON ENERGY AND CONTACT
The value of the Gross-Pitaevskii functional (7) evaluated at the ground-state wave function determines the polaron energy E = E[φ] relative to the homogeneous BEC. The energy is shown in Fig. 3(a) as a function of the dimensionless impurity-boson interaction ξ/a IB . It is always negative for an attractive impurity potential. In Fig. 3(b) we present the energy for a mobile impurity of arbitrary mass in units of E n which depends on the BEC density n 0 . Specifically, we show results for two BEC gas parameters, which for equal mass of boson and impurity, i.e., m red /m B = 1/2, correspond to values n 0 a 3 BB = (4π) −3 = 5.0 × 10 −4 (n 0 ξ 3 = 1) at stronger boson repulsion and n 0 a 3 BB = (8π) −3 = 0.63 × 10 −4 (n 0 ξ 3 = 2.8) at weaker boson repulsion. We find that, in absolute terms, the polaron binding energy is larger for weaker boson repulsion where the BEC can be more strongly deformed and thus acquire more attractive potential energy.
At weak impurity-bath attraction 1/a IB n 1/3 0 −1 the polaron energy approaches the mean-field value irrespective of a BB . Within mean-field theory the polaron energy diverges to −∞ at unitarity 1/a IB = 0. Variational approaches based on the trunctated Hamiltonian Eq. (4) predict that the inclusion of Bogoliubov corrections is not sufficient to heal this instability, but rather results in a shift of the instability to the repulsive side of the Feshbach resonance [30]. In contrast, we find that going beyond the Bogoliubov approximation by working with the full Hamiltonian (1) allows the boson repulsion to stabilize the polaron at a finite ground-state energy that smoothly crosses over from the attractive to the repulsive side of the Feshbach resonance. The deformation of the BEC is also reflected in the number of bosons participating in the formation of the polaron dressing cloud As shown in Fig. 3(c), for a smaller gas parameter (orange dashed line) the impurity attracts a larger polaron cloud because the bosons are less repulsive. Naturally, this larger dressing cloud corresponds to the larger polaron binding energy found in Fig. 3(b).
Our results in Fig. 3 are shown for a constant range R across the Feshbach resonance as applicable to experiments where the microscopic range of interactions can typically not be tuned synchronously with the scattering length. Thus the effective range r eff varies in dependence on a IB : The inset of Fig. 3(c) shows the effective range that is obtained from the scattering phase shift using Eq. (11). At constant potential range Rn 1/3 0 = 1, the effective range r eff reaches a minimum on the repulsive side of the resonance (a IB > 0) and grows towards both weak-coupling limits a IB → 0 ± .
The variation of the polaron energy with a IB defines the impurity contact parameter [38,40,42,51,59] which characterizes the impurity-boson correlations g (2) IB (r) at short distances outside the impurity potential: . (15) The contact is shown in Fig. 4: we find that it reaches a maximal value on the repulsive side of the resonance. Similar to the energy, we find that the contact is larger for smaller boson repulsion a BB (blue solid line) where the BEC is more strongly deformed and thus g IB (r) is enhanced (see Fig. 2). At weak attractive interaction the contact approaches the ground-state value of an impurity in an ideal BEC (a BB = 0) [38] C mf = 16π 2 n 0 a 2 IB , which agrees with the mean-field result prediction.

V. UNIVERSALITY
Finally, we test the notion of universality of the Bose polaron by studying different shapes and ranges R of the impurity potential. Generally, we work in the regime where the potential range is much larger than the boson scattering length, R a BB , where local Gross-Pitaevskii theory has been shown to be applicable [35]. To test universality in this regime, we specifically compare the predictions following from the Gaussian potential V gauss (r) in Eq. (9), and the exponential potential V expon (r) in Eq. (10) for various ranges R.
We find that when tuning the depth of the potentials to obtain equal scattering length a IB = −∞ at equal range R, the polaron energies are very different. Instead, if R is tuned to yield the same effective range r eff for both potentials, as shown in Fig. 5(a), remarkable agreement between the polaron energies is found. Indeed, as shown in Fig. 5(b), for both potentials, the polaron energy approximately follows a linear scaling law for r eff 0.2ξ, Remarkably, this universality is found not only for the purely attractive Gaussian and exponential potentials but extends also to interactions featuring a repulsive contribution such as a Lennard-Jones potential V LJ (r) ∼ λR 10 /r 12 − R 4 /r 6 , as illustrated in Fig. 5(b). This shows that in the Bose polaron problem the momentumdependent scattering phase shift is probed in the regime where the effective range expansion is valid. Our result complements a recent GPE study which found a power-law scaling of the unitary polaron energy at ranges shorter than the healing length, E/E n ∼ (r eff /ξ) 1/3 for a BB r eff ξ in the case of a square-well potential [42]. The GPE approach can be extended to even shorter ranges r eff a BB by using a nonlocal generalization of Gross-Pitaevskii theory [40].
Also on the attractive side of the resonance, see Fig. 5(c) for 1/a IB n 1/3 0 = −1, we find a universal effective range dependence of the polaron energy, albeit not a linear one. Note that to cover the domain of effective ranges shown in Fig. 5(d), for our parameters the Gaussian potential, unlike the exponential, must be tuned over such depths that it supports more than one bound state. Consequently, there is no longer a unique mapping from r eff to E. This is illustrated in Fig. 5(d) for the repulsive side, at 1/a IB n 1/3 0 = 1: universality can at best hold in the vicinity of the first bound state, and only for a limited interval of r eff values. Furthermore, since the Efimov effect modifies the polaron energy spectrum [2,3,49], universality can also depend on the three-body parameter [50,[60][61][62].

VI. DISCUSSION
The variational principle provides a powerful tool to compute both ground-state and dynamical properties of quantum many-body systems. However, in order to make full use of its predictive power it is essential to understand the limitations of approximations applied to the Hamiltonian to be analyzed. In this regard the strongcoupling Bose polaron is a case in point. The full Bose polaron Hamiltonian (1) -and equivalently Eq. (3)-is bounded from below for repulsive boson interaction, and hence variational wave functions give rigorous bounds on the ground-state energy. For finite Bose repulsion, the polaron energy remains finite for any Bose-impurity scattering length, including resonant interactions, and the ground state represents a strong-coupling Bose polaron that is self-stabilized by its own dressing cloud.
In contrast, when applying the Bogoliubov approximation to the full model (1) by truncating terms of higherthan quadratic order in the boson operators, the resulting, truncated HamiltonianĤ in Eq. (4) becomes unbound from below. Crucially, this results in an instability of the Bose polaron problem that is solely an artefact of this approximation.
Quite remarkably, the instability of the truncated HamiltonianĤ becomes, however, only evident when considering wavefunctions that account for more than two phonon excitations from the homogenous BEC. For instance, for a Chevy-type wavefunction [25,27,28] that itself is truncated at the single excitation level, the terms beyond quadratic order in the repulsive interactions have vanishing expectation value. Thus, incidentally, the Chevy ansatz yields the same prediction when applied to both the full and the truncated model. Thus due to its tremendous simplicity the Chevy ansatz becomes immune to the instability of approximate HamiltonianĤ .
However, while being a well-defined approach, the simple Chevy ansatz misses the fact that at the weak Boson repulsion (as typically present in cold gases) the polaron cloud-even within the full model Eq. (1)-can contain an exceedingly large number of bosons ( Fig. 3(b)). Such a large local deformation of the BEC is naturally captured by the inhomogeneous wave function (6). Crucially, while accounting for an arbitrary number of boson excitations, when applied to the full Hamiltonian (1) it still leads to a bounded energy functional (7). Its solution shows that the smaller the boson repulsion and the wider the impurity potential, the larger the polaron cloud becomes. The product state approach is complementary to the Chevy ansatz including its extensions to multi-boson excitations [3,49,63], and it becomes particularly accurate for soft impurity potentials where it is justified to ignore bosonic correlations. Remarkably, the case of extremely soft potentials is realized with Rydberg excitation immersed in BECs. In this case it was predicted that up to hundreds of atoms can be bound to the single impurity leading to the creation of Rydberg polarons [64]. Since for Rydberg impurities the range of interactions R dramatically exceeds the interparticle distance, our local GP theory applies and provides a so far missing explanation as to why the experimental observation of Rydberg polarons [29] is described exceptionally well by a coherent state approach [26].
Recently also first steps to the understanding of the complementary, intermediate regime of short-range impurity potentials -with yet large dressing clouds-has been achieved by using an extension to nonlocal Gross-Pitaevskii theory [40]. In conjunction with our present result, these combined new approaches resolve a fundamental shortcoming of the Bogoliubov approximation: while the formulation of the interacting Bose gas in terms of Bogoliubov quasiparticles is exact, the additional approximation to neglect the residual interaction between phonons is not. In particular, the quadratic Bogoliubov mean-field Hamiltonian is unbounded from below for the strong-coupling Bose polaron, most obviously in the regime where a two-body bound state appears on the repulsive side of the resonance [30,31].
As discussed above, the Chevy-type ansatz applied to the truncated the Bogoliubov Hamiltonian [28] still yields a finite polaron energy since it is of such low order in boson excitations that it is not sensitive to the truncated part of the Hamiltonian. However, when the coherent state ansatz or higher-order excitation extensions of the Chevy ansatz are applied to the truncated Hamiltonian (4) they can lead to divergencies in the ground-state energy. In the case of the coherent state ansatz [30,36] the divergence is due to the large occupation of excitations in the vicinity of the impurity not counteracted by boson repulsion. Instead the local extension to the truncated Bogoliubov approach studied in the present work (see also [42,50,51]) as well as its nonlocal counterpart [40] provides a stable starting point for strong-coupling Bose polaron dynamics, and we showed how it can find an effective description in terms of a renormalization of the impurity-boson potential (Fig. 2).
Beyond our treatment of the two-particle impurityboson correlations, three-body and higher-order correlations give rise to the Efimov effect and three-body recombination. The Efimov effect can occur either between one impurity and two bosons [3,49,63] or between two impurities and one boson [2,6]. These few-body effects can be captured by Gaussian variational wave functions [55,56]. Alternatively, extensions of the Chevy ansatz to two or more independent bosonic excitations [3] can be employed. The latter approach was applied in the analysis of the truncated model (4) and universal scaling depending on the three-body parameter was found [49,63]. It remains an interesting open question how universal three-body physics carries over to the many-body case in a dense bosonic medium when the full Hamiltonian Eq. (1) is considered [65].
Finally, in ultracold atomic gases the boson repulsion originates from an attractive van-der-Waals potential between atoms. This results in the existence of deeply bound molecular states into which atoms can decay in three-body recombination. These deeply bound states are neither accounted for in variational approaches nor in quantum Monte Carlo [43]. Up to now most experiments probe the metastable Bose polaron state of matter on transient time scales where the deeply bound statesarising both from the fundamentally attractive Bose-Bose and Bose-impurity potentials-can be ignored, and thus the stable variational approach discussed in our work is well applicable. However, as one starts to explore longer time scales or the build-up of more complex correlated states of impurities, understanding the impact of the fundamental dissipative nature arising from deeply bound states becomes essential and requires the development of new theoretical approaches to quantum impurity problems.