Strong coupling Bose polarons in a BEC

TL;DR

This paper employs a non-perturbative renormalization group method to analyze the Bose polaron problem in a BEC, revealing phase transitions, effective mass behavior, and quantum fluctuation effects across interaction regimes.

Contribution

It provides a unified theoretical framework for Bose polarons, including a detailed phase diagram and insights into the polaron-to-molecule transition near Feshbach resonances.

Findings

Prediction of a sharp increase in effective mass near Feshbach resonance

Identification of a transition to states below polaron energy at weaker interactions

Logarithmic corrections to polaron energy due to quantum fluctuations

Abstract

We use a non-perturbative renormalization group approach to develop a unified picture of the Bose polaron problem, where a mobile impurity is strongly interacting with a surrounding Bose-Einstein condensate (BEC). A detailed theoretical analysis of the phase diagram is presented and the polaron-to-molecule transition is discussed. For attractive polarons we argue that a description in terms of an effective Fr\"ohlich Hamiltonian with renormalized parameters is possible. Its strong coupling regime is realized close to a Feshbach resonance, where we predict a sharp increase of the effective mass. Already for weaker interactions, before the polaron mass diverges, we predict a transition to a regime where states exist below the polaron energy and the attractive polaron is no longer the ground state. On the repulsive side of the Feshbach resonance we recover the repulsive polaron, which has a finite lifetime because it can decay into low-lying molecular states. We show for the entire range of couplings that the polaron energy has logarithmic corrections in comparison with predictions by the mean-field approach. We demonstrate that they are a consequence of the polaronic mass renormalization which is due to quantum fluctuations of correlated phonons in the polaron cloud.

Tables

Table 1. Table 1: Dimensional analysis of quantum fluctuations around the MF polaron state at low energies, k ≪ 1 / ξ much-less-than 𝑘 1 𝜉 k\ll 1/\xi . Engineering dimensions were used for deriving this table. Relevant and marginal terms are highlighted.

scaling for $k\ll 1/\xi$	$d=1$	$d=2$	$d=3$
${\hat{a}}_{𝒌}$	${\Lambda }^{-1}$	${\Lambda }^{-3/2}$	${\Lambda }^{-2}$
$\frac{1}{M}{\int }^{\Lambda }{d}^{d}k{d}^d{k}^{\prime }𝒌\cdot {𝒌}^{\prime }:{\hat{\Gamma }}_{𝒌}{\hat{\Gamma }}_{{𝒌}^{\prime }}:$	${𝚲}^{\mathrm{𝟎}}$	${𝚲}^{\mathrm{𝟎}}$	${𝚲}^{\mathrm{𝟎}}$
$G_{+}{\int }^{\Lambda }{d}^{d}k{d}^d{k}^{\prime }:\delta {\hat{n}}_{𝒌}\delta {\hat{n}}_{{𝒌}^{\prime }}:$	${\Lambda }^1$	${\Lambda }^2$	${\Lambda }^3$
$G_{-}{\int }^{\Lambda }{d}^{d}k{d}^d{k}^{\prime }:{\hat{\vartheta }}_{𝒌}{\hat{\vartheta }}_{{𝒌}^{\prime }}:$	${𝚲}^{-\mathrm{𝟏}}$	${𝚲}^{\mathrm{𝟎}}$	${\Lambda }^1$

Table 2. Table 2: The same analysis as in Tab. 3 is shown, but for the high energy regime where k ≫ 1 / ξ much-greater-than 𝑘 1 𝜉 k\gg 1/\xi .

scaling for $k\gg 1/\xi$	$d=1$	$d=2$	$d=3$
${\hat{a}}_{𝒌}$	${\Lambda }^{-3/2}$	${\Lambda }^{-2}$	${\Lambda }^{-5/2}$
$\frac{1}{M}{\int }^{\Lambda }{d}^{d}k{d}^d{k}^{\prime }𝒌\cdot {𝒌}^{\prime }:{\hat{\Gamma }}_{𝒌}{\hat{\Gamma }}_{{𝒌}^{\prime }}:$	${𝚲}^{-\mathrm{𝟑}}$	${𝚲}^{-\mathrm{𝟐}}$	${𝚲}^{-\mathrm{𝟐}}$
$G_{+}{\int }^{\Lambda }{d}^{d}k{d}^d{k}^{\prime }:\delta {\hat{n}}_{𝒌}\delta {\hat{n}}_{{𝒌}^{\prime }}:$	${𝚲}^{-\mathrm{𝟏}}$	${𝚲}^{\mathrm{𝟎}}$	${\Lambda }^1$
$G_{-}{\int }^{\Lambda }{d}^{d}k{d}^d{k}^{\prime }:{\hat{\vartheta }}_{𝒌}{\hat{\vartheta }}_{{𝒌}^{\prime }}:$	${𝚲}^{-\mathrm{𝟏}}$	${𝚲}^{\mathrm{𝟎}}$	${\Lambda }^1$

Table 3. Table 3: Dimensional analysis of quantum fluctuations around the MF polaron state at low energies for k ≪ 1 / ξ much-less-than 𝑘 1 𝜉 k\ll 1/\xi . The RG flow of the coupling constants is included and the regime where neither of the coupling constants G ± subscript 𝐺 plus-or-minus G_{\pm} diverges is considered. Relevant and marginal terms are highlighted.

scaling for $k\ll 1/\xi$	$d=1$	$d=2$	$d=3$
${\hat{a}}_{𝒌}$	${\Lambda }^{-1}$	${\Lambda }^{-3/2}$	${\Lambda }^{-2}$
$\frac{1}{ℳ}{\int }^{\Lambda }{d}^{d}k{d}^d{k}^{\prime }𝒌\cdot {𝒌}^{\prime }:{\hat{\Gamma }}_{𝒌}{\hat{\Gamma }}_{{𝒌}^{\prime }}:$	${𝚲}^{\mathrm{𝟎}}$	${𝚲}^{\mathrm{𝟎}}$	${𝚲}^{\mathrm{𝟎}}$
$G_{+}{\int }^{\Lambda }{d}^{d}k{d}^d{k}^{\prime }:\delta {\hat{n}}_{𝒌}\delta {\hat{n}}_{{𝒌}^{\prime }}:$	${\Lambda }^1$	${\Lambda }^2$	${\Lambda }^3$
$G_{-}{\int }^{\Lambda }{d}^{d}k{d}^d{k}^{\prime }:{\hat{\vartheta }}_{𝒌}{\hat{\vartheta }}_{{𝒌}^{\prime }}:$	${𝚲}^{\mathrm{𝟎}}$	${𝚲}^{\mathrm{𝟎}}$	${\Lambda }^1$

Table 4. Table 4: The same as Tab. 3 but in the high energy regime where k ≫ 1 / ξ much-greater-than 𝑘 1 𝜉 k\gg 1/\xi .

scaling for $k\gg 1/\xi$	$d=1$	$d=2$	$d=3$
${\hat{a}}_{𝒌}$	${\Lambda }^{-3/2}$	${\Lambda }^{-2}$	${\Lambda }^{-5/2}$
$\frac{1}{ℳ}{\int }^{\Lambda }{d}^{d}k{d}^d{k}^{\prime }𝒌\cdot {𝒌}^{\prime }:{\hat{\Gamma }}_{𝒌}{\hat{\Gamma }}_{{𝒌}^{\prime }}:$	${𝚲}^{\mathrm{𝟎}}$	${𝚲}^{\mathrm{𝟎}}$	${𝚲}^{-\mathrm{𝟏}}$
$G_{+}{\int }^{\Lambda }{d}^{d}k{d}^d{k}^{\prime }:\delta {\hat{n}}_{𝒌}\delta {\hat{n}}_{{𝒌}^{\prime }}:$	${𝚲}^{\mathrm{𝟎}}$	${𝚲}^{\mathrm{𝟎}}$	${\Lambda }^1$
$G_{-}{\int }^{\Lambda }{d}^{d}k{d}^d{k}^{\prime }:{\hat{\vartheta }}_{𝒌}{\hat{\vartheta }}_{{𝒌}^{\prime }}:$	${𝚲}^{\mathrm{𝟎}}$	${𝚲}^{\mathrm{𝟎}}$	${\Lambda }^1$

Equations

$$\hat{\mathcal{H}}_{\rm B}=g_{\text{IB}}^{(0)}n_{0}+\int^{\Lambda_{0}}d^{3}k~{}\Biggl{[}\omega_{k}\hat{a}^{\dagger}_{\bm{k}}\hat{a}_{\bm{k}}+V_{k}e^{i\bm{k}\cdot\hat{\bm{x}}}\left(\hat{a}_{\bm{k}}+\hat{a}^{\dagger}_{-\bm{k}}\right)\Biggr{]}+\frac{\hat{\bm{p}}^{2}}{2M}+\\ +\frac{g_{\text{IB}}^{(0)}}{(2\pi)^{3}}\int^{\Lambda_{0}}d^{3}k~{}d^{3}k^{\prime}~{}e^{i(\bm{k}-\bm{k}^{\prime})\cdot\hat{\bm{x}}}~{}\left(\cosh\theta_{k}\hat{a}^{\dagger}_{\bm{k}}-\sinh\theta_{k}\hat{a}_{-\bm{k}}\right)\left(\cosh\theta_{k^{\prime}}\hat{a}_{\bm{k}^{\prime}}-\sinh\theta_{k^{\prime}}\hat{a}^{\dagger}_{-\bm{k}^{\prime}}\right).$$

$${\hat{d}}_{\mathbf{k}}=\cosh {\theta }_k{\hat{a}}_{\mathbf{k}}-\sinh {\theta }_k{\hat{a}}_{-\mathbf{k}}^{†},$$

$$W_k=\cosh {\theta }_k-\sinh {\theta }_k={\left(\frac{k^2{\xi }^2}{2+k^2{\xi }^2}\right)}^{1/4}.$$

$$\frac{1}{g_{\text{IB}}^{(0)}}=\frac{m_{\mathrm{red}}}{2\pi a_{\text{IB}}}-{\int }^{{\Lambda }_0}\frac{d^{3}k}{(2\pi {)}^3}\frac{2m_{\mathrm{red}}}{k^2},$$

$$g_{\text{IB}}^{(0)}=-\frac{{\pi }^2}{m_{\mathrm{red}}}{\Lambda }_0^{-1}-\frac{{\pi }^3}{2m_{\mathrm{red}}{a}_{\text{IB}}}{\Lambda }_0^{-2}.$$

$$\hat{\mathcal{H}}_{\rm LLP}=g_{\text{IB}}^{(0)}n_{0}+\int^{\Lambda_{0}}d^{3}k~{}\Biggl{[}\omega_{k}\hat{a}^{\dagger}_{\bm{k}}\hat{a}_{\bm{k}}+V_{k}\left(\hat{a}_{\bm{k}}+\hat{a}^{\dagger}_{\bm{k}}\right)\Biggr{]}+\frac{1}{2M}\left(\bm{p}-\int^{\Lambda_{0}}d^{3}k~{}\bm{k}\hat{a}^{\dagger}_{\bm{k}}\hat{a}_{\bm{k}}\right)^{2}+\\ +\frac{g_{\text{IB}}^{(0)}}{(2\pi)^{3}}\int^{\Lambda_{0}}d^{3}k~{}d^{3}k^{\prime}~{}\ \left(\cosh\theta_{k}\hat{a}^{\dagger}_{\bm{k}}-\sinh\theta_{k}\hat{a}_{-\bm{k}}\right)\left(\cosh\theta_{k^{\prime}}\hat{a}_{\bm{k}^{\prime}}-\sinh\theta_{k^{\prime}}\hat{a}^{\dagger}_{-\bm{k}^{\prime}}\right).$$

$$\begin{gathered} {\hat{\mathcal{H}}}_{\mathrm{toy}}=g_{\text{IB}}{n}_0+\omega {\hat{a}}^{†}\hat{a}+g_{\text{IB}}\sqrt{n_0}W\left({\hat{a}}^{†}+\hat{a}\right)+ \\ +\frac{g_{\text{IB}}}{4}\left[(W^2+W^{-2}){\hat{a}}^{†}\hat{a}+(W^2-W^{-2})\hat{a}\hat{a}+\text{h.c.}\right]. \end{gathered}$$

$${\hat{U}}_{\text{MF}}=\exp \left({\alpha }^{\text{MF}}{\hat{a}}^{†}-\text{h.c.}\right),$$

$$\begin{gathered} {\tilde{\mathcal{H}}}_{\mathrm{toy}}={\hat{U}}_{\text{MF}}^{†}{\hat{\mathcal{H}}}_{\mathrm{toy}}{\hat{U}}_{\text{MF}}={\beta }_{\text{MF}}{g}_{\text{IB}}{n}_0+\omega {\hat{a}}^{†}\hat{a}+ \\ +\frac{g_{\text{IB}}}{4}\left[(W^2+W^{-2}){\hat{a}}^{†}\hat{a}+(W^2-W^{-2})\hat{a}\hat{a}+\text{h.c.}\right]. \end{gathered}$$

$$\delta \hat{n}$$

$$\hat{\vartheta }$$

$${\tilde{\mathcal{H}}}_{\mathrm{toy}}=n_0{W}^2{A}_{-}{\hat{\vartheta }}^2+\frac{A_{+}}{4n_0{W}^2}\delta {\hat{n}}^2$$

$$A_{\pm }=\omega +g_{\text{IB}}{W}^{\pm 2}.$$

$$g_{\text{IB},\pm }=-\omega /W^{\pm 2},$$

$$|{\psi }_{\mathrm{MF}}⟩={\hat{U}}_{\text{MF}}|0⟩,$$

$${\hat{U}}_{\text{MF}}[{\alpha }_{\mathbf{k}}^{\text{MF}}]=\exp \left({\int }^{{\Lambda }_0}{d}^{3}k{\alpha }_{\mathbf{k}}^{\text{MF}}{\hat{a}}_{\mathbf{k}}^{†}-\text{h.c.}\right).$$

$$\frac{\delta \mathcal{H}[{\alpha }_{\mathbf{k}}^{\text{MF}}]}{\delta {\alpha }_{\mathbf{k}}^{\text{MF}}}=0$$

$${\alpha }_k=-{\beta }_{\text{MF}}\frac{V_k}{{\Omega }_k^{\text{MF}}},{\Omega }_k^{\text{MF}}={\omega }_k+\frac{k^2}{2M},$$

$${\beta }_{\text{MF}}={\left[1+\frac{g_{\text{IB}}^{(0)}}{(2\pi {)}^3}{\int }^{{\Lambda }_0}{d}^{3}k\frac{W_k^2}{{\Omega }_k^{\text{MF}}}\right]}^{-1},$$

$$E_0^{\text{MF}}=n_0{\beta }_{\text{MF}}{g}_{\text{IB}}^{(0)}+\frac{g_{\text{IB}}^{(0)}}{(2\pi {)}^3}{\int }^{{\Lambda }_0}{d}^{3}k{\sinh }^2{\theta }_k.$$

$$E_0^{\text{MF}}=n_0{\beta }_{\text{MF}}{g}_{\text{IB}}^{(0)}=n_0\frac{2\pi }{m_{\mathrm{red}}}{\left(1/a_{\text{IB}}+1/a_{\text{IB},+}^{\text{MF}}\right)}^{-1},$$

$$a_{\text{IB},+}^{\text{MF}}=\pi {\left[2{\Lambda }_0-m_{\mathrm{red}}^{-1}{\int }_0^{{\Lambda }_0}dk{k}^2\frac{W_k^2}{{\Omega }_k^{\text{MF}}}\right]}^{-1}\ge 0.$$

$$\mathcal{H}[{\alpha }_{\mathbf{k}}]=E_0^{\text{MF}}+\frac{1}{2}\int d^{3}k{d}^3{k}^{\prime }\delta \text{uline}{\alpha }_{\mathbf{k}}^{†}\text{uuline}H(\mathbf{k},{\mathbf{k}}^{\prime })\delta \text{uline}{\alpha }_{{\mathbf{k}}^{\prime }},$$

$$\delta {\hat{n}}_{\mathbf{k}}$$

$${\hat{\vartheta }}_{\mathbf{k}}$$

$$\begin{gathered} \tilde{\mathcal{H}}=E_0^{\text{MF}}+{\int }^{{\Lambda }_0}{d}^{3}k{d}^3{k}^{\prime }\frac{\mathbf{k}\cdot {\mathbf{k}}^{\prime }}{2M}:{\hat{\Gamma }}_{\mathbf{k}}{\hat{\Gamma }}_{{\mathbf{k}}^{\prime }}:+ \\ +{\int }^{{\Lambda }_0}{d}^{3}k{\Omega }_k^{\text{MF}}\left[\frac{:\delta {\hat{n}}_{-\mathbf{k}}\delta {\hat{n}}_{\mathbf{k}}:}{4n_0{W}_k^2}+n_0{W}_k^2:\delta {\hat{\vartheta }}_{\mathbf{k}}\delta {\hat{\vartheta }}_{{\mathbf{k}}^{\prime }}:\right]+ \\ +{\int }^{{\Lambda }_0}{d}^{3}k{d}^3{k}^{\prime }\left[\frac{G_{+}}{2n_0}:\delta {\hat{n}}_{\mathbf{k}}\delta {\hat{n}}_{{\mathbf{k}}^{\prime }}:+2n_0{G}_{-}:{\hat{\vartheta }}_{\mathbf{k}}{\hat{\vartheta }}_{{\mathbf{k}}^{\prime }}:\right]. \end{gathered}$$

$${\hat{\Gamma }}_{\mathbf{k}}={\hat{a}}_{\mathbf{k}}^{†}{\hat{a}}_{\mathbf{k}}+{\alpha }_k^{\text{MF}}\left({\hat{a}}_{\mathbf{k}}^{†}+{\hat{a}}_{\mathbf{k}}\right).$$

$$\hat{\psi }(\mathbf{r})\approx \sqrt{n_0+\delta \hat{n}(\mathbf{r})}{e}^{i\hat{\vartheta }(\mathbf{r})},$$

$$[\delta {\hat{n}}_{\mathbf{k}},\delta {\hat{\vartheta }}_{-{\mathbf{k}}^{\prime }}]=i\delta (\mathbf{k}-{\mathbf{k}}^{\prime }).$$

$$\int d^{3}k|\delta {\alpha }_k{|}^2\stackrel{!}{=}1.$$

$$\frac{\delta }{\delta {\vartheta }_k}\left(\mathcal{H}[{\vartheta }_k]-{\mu }_{-}{n}_0{W}_k^2\int d^{3}k{\vartheta }_k^2\right)=0.$$

Full text

I Introduction

The dressing of mobile impurities with collective excitations of a surrounding many-body system – i.e. the formation of polarons Landau1946 ; Landau1948 – is a ubiquitous phenomenon in physics. Its consequences, including for example the enhancement of the effective mass PolaronsAdvMat2007 ; alexandrov2009advances , have been observed in a large class of systems, ranging from electrons in a semiconductor coupled to phonon excitations in the host lattice Froehlich1954 ; Mahan2000 and exciton polaritons interacting with a surrounding Fermi sea Sidler2016 to mixtures of ultracold fermions Schirotzek2009 ; Koschorreck2012 and bosons Catani2012 ; Fukuhara2013 ; Scelle2013 ; Jorgensen2016PRL ; Hu2016PRL with contact interactions.

Understanding the physics of polarons at intermediate couplings, where the phonon excitations in the polaron cloud strongly correlated, has posed a challenge for decades Feynman1955 ; Devreese2013 ; Grusdt2015Varenna . Recently, experiments with mixtures of ultracold bosons Catani2012 ; Hu2016PRL ; Jorgensen2016PRL have realized this regime and they have initiated a detailed investigation of Bose polarons for a large range of couplings. Despite intense theoretical effort in understanding the physics of mobile impurities in ultracold quantum gases Mathey2004 ; Bruderer2007 ; BeiBing2009 ; Tempere2009 ; Casteels2011 ; Casteels2012 ; Rath2013 ; Shashi2014RF ; Kain2014 ; Li2014 ; Grusdt2014BO ; Yin2015PRA ; Grusdt2015RG ; Grusdt2016RG ; Vlietinck2015 ; Levinsen2015 ; Christensen2015 ; Ardila2015 ; Grusdt2016 ; Shchadilova2016 ; Shchadilova2016PRL ; Parisi2016 ; Kain2016 , a number of questions remains unresolved so far.

Recent calculations of the spectral function of the Bose polaron based on a time-dependent coherent state approach Shchadilova2016PRL , closely related to the mean-field (MF) description of the polaron Lee1953 ; Shashi2014RF , suggested the existence of at least two different parameter regimes. While the attractive polaron corresponds to a single sharp feature in the spectrum clearly visible at weak couplings, an entire set of lines associated with molecular bound states was predicted on the repulsive side of the Feshbach resonance at energies below the repulsive polaron peak. It remained unclear, however, whether these regimes are connected by a smooth crossover or a sharp phase transition corresponding to a non-analytical ground state energy. On the one hand the energy of the MF polaron state diverges at a critical interaction strength Shchadilova2016PRL , suggesting a sharp transition. On the other hand, variational wavefunctions including only a few excitations Li2014 ; Levinsen2015 ; Jorgensen2016PRL as well as closely related T-matrix calculations Rath2013 predict a smooth crossover of the ground states between the two regimes.

Another question concerns the dependence of the polaron energy on the effective range associated with the impurity-boson scattering potential, which represents the first correction to the scattering amplitude beyond pure contact interactions. For example, the ground state energy of an interacting Bose gas has a sub-leading contribution due to quantum fluctuations which scales logarithmically with the effective range of the boson-boson scattering Wu1959 ; Sawada1959 ; Hugenholtz1959 . Similar corrections were predicted for impurities by a renormalization group (RG) analysis of the effective Fröhlich polaron model Grusdt2015RG and by analytical calculations using perturbation theory Christensen2015 .

Here we develop a unified description of Bose polarons at strong couplings, and in the limit of a weakly interacting BEC. Our analysis is based on a non-perturbative RG approach which we use to derive the polaron phase diagram. In particular, we investigate one of the theoretically most challenging regimes, where the polaron cloud contains many phonons which become correlated due to the interactions induced by the mobile impurity. An accurate description of these correlations requires either numerical quantum Monte Carlo calculations Prokofev1998 ; Ardila2015 ; Vlietinck2015 ; Parisi2016 , variational wavefunctions with many degrees of freedom Shchadilova2016 ; Kain2016 or an RG approach Grusdt2015RG ; Grusdt2015Varenna ; Grusdt2016RG .

As a main result we find evidence for two sharp transitions where the nature of the polaron changes abruptly. Our work extends earlier calculations, which were either based on the effective Fröhlich Hamiltonian Mathey2004 ; Bruderer2007 ; BeiBing2009 ; Tempere2009 ; Casteels2011 ; Casteels2012 ; Shashi2014RF ; Kain2014 ; Li2014 ; Grusdt2014BO ; Yin2015PRA ; Grusdt2015RG ; Grusdt2016RG ; Vlietinck2015 ; Grusdt2016 or relied on approximate trial wavefunctions Li2014 ; Levinsen2015 ; Jorgensen2016PRL . We include two-phonon terms going beyond the Fröhlich model Rath2013 and treat them accurately up to high energies Rath2013 ; Shchadilova2016PRL . Our RG analysis moreover reveals that the Fröhlich terms in the Hamiltonian are the most relevant ones for describing the Bose polaron away from the resonance where multiple bound states appear, although corrections are important for making quantitative predictions.

The paper is organized as follows. We start by summarizing our results in Sec. II. We proceed by introducing the model describing an impurity in a BEC in Sec. III. A simplified toy model is analyzed in Sec. IV to gain intuitive understanding for the polaron phase diagram. In Sec. V we summarize the mean-field theory on which the RG method is based, and discuss molecular bound states. In Sec. VI we present our RG method and the RG flow equations. Regularization schemes for divergencies of the coupling constants are also introduced. In Sec. VII we present calculations of the polaron energies which we compare to predictions by the simpler Fröhlich model. Corrections to the phase diagram by the RG are discussed in Sec. VIII. We close with a summary and by giving an outlook in Sec. IX.

II Summary of results

The focus of this paper is to understand the properties of the Bose polaron, which corresponds to a sharp peak in the spectral function $I(\omega)$ at weak couplings. Changes in the nature of the Bose polaron at strong couplings are accompanied by strong modifications of the spectrum, which we discuss next. In Fig. 1 we summarize our main results and sketch the shape of $I(\omega)$ for different impurity-boson scattering lengths $a_{\rm IB}$, which can be tuned by a Feshbach resonance in experiments with cold atoms.

(I) Attractive polarons.– On the attractive side sufficiently far from the resonance, i.e. for $1/a_{\text{IB}}\ll-1/\xi$ less than the inverse of the healing length $\xi$ of the BEC, the spectrum consists of a delta function peak $Z\delta(\hbar\omega-E_{0})$ located at the energy $E_{0}$ of the attractive polaron ($①$ in Fig. 1). The corresponding quasiparticle residue $Z\approx 1$ is close to one far from the Feshbach resonance (for finite $\xi$). All the remaining spectral weight, $1-Z$, is located at energies $E>E_{0}$ above the polaron energy and corresponds to shake-off processes of phonons from the polaron cloud as discussed in Refs. Rath2013 ; Shashi2014RF . This regime can be described within the Fröhlich model.

For stronger interactions, but still on the attractive side, an increasing amount of spectral weight moves into the incoherent part of the spectrum. Otherwise the shape of the spectral function corresponding to the attractive polaron remains unchanged ($②$ in Fig. 1). With increasing interaction strength the quantitative deviations from predictions by the Fröhlich model begin to grow.

(II) Attractive polarons at strong coupling.– At a critical value $a_{\rm IB}=a_{\text{IB},-}^{\rm RG}<0$ on the attractive side of the resonance, the shape of the spectrum changes qualitatively ($③$ in Fig. 1). Here a phase transition takes place and the ground state is no longer defined. We emphasize that this conclusion relies on the Bogoliubov approximation for describing the host BEC atoms, which does not include interactions between Bogoliubov quasiparticles. Within this approximation we predict states at arbitrarily negative energies, giving rise to a finite spectral weight $I(\omega)>0$ even for $\omega<E_{0}$ below the polaron energy. Note however that this spectral weight quickly diminishes for $\omega$ below $E_{0}$. The value $a_{\text{IB},-}^{\rm RG}\sim\xi$ is of the order of the BEC healing length $\xi$ and depends on the impurity mass $M$; In particular $|a_{\text{IB},-}^{\rm RG}|=\infty$ for a localized impurity, $M=\infty$.

In this regime beyond $a_{\text{IB},-}^{\rm RG}$, the peak in the spectrum at energy $E_{0}$ corresponding to the attractive polaron acquires a finite linewidth. It can be associated with the inverse lifetime of the polaron, which can decay into the states appearing at energies below $E_{0}$. In the spectral function shown in Fig. 1, the broadened attractive polaron peak can not be distinguished from the even broader incoherent part of the spectrum associated with excitations of the polaron cloud at energies above $E_{0}$. Note, however, that these shake-off processes at high energies give rise to an asymmetric shape of $I(\omega)$ around $\omega=E_{0}$ which is visible in Fig. 1.

For small but non-vanishing phonon-phonon interactions, beyond the Bogoliubov approximation considered here, we expect a continuous transition into this regime around $a_{\text{IB},-}^{\rm RG}$. When such phonon non-linearities are sufficiently weak, however, we think that the features in the spectral function discussed above remain valid.

When the Feshbach resonance is approached further from the attractive side, the energy $E_{0}$ of the attractive polaron starts to decrease dramatically and diverges at a second critical value $a_{\rm IB}=a_{\text{IB},+}^{\rm RG}$. At the same time the spectral weight in the polaron peak around $E_{0}$, broadened by the finite lifetime, becomes strongly suppressed as a consequence of the dressing with a diverging number of phonons, see Fig. 2. Meanwhile the incoherent part of the spectrum gains weight. The polaron mass also diverges at $a_{\text{IB},+}^{\rm RG}$, see Fig. 3.

These features suggest another transition located at $a_{\text{IB},+}^{\rm RG}$, where the nature of eigenstates at low energies changes abruptly. This transition is also predicted by MF theory Shchadilova2016PRL , albeit at a different interaction strength $a_{\rm IB}=a_{\rm IB,+}^{\text{MF}}$, where the polaron changes from attractive to repulsive. The shift of the transition to $1/a_{\rm IB,+}^{\rm RG}<1/a_{\rm IB,+}^{\text{MF}}$ predicted by the RG is due to additional polaronic mass renormalization which is not taken into account by the MF theory.

(III) Instability of the RG.– Between the two critical values predicted by the RG approach and MF theory, i.e. for $1/a_{\rm IB,+}^{\rm RG}<1/a_{\rm IB}<1/a_{\rm IB,+}^{\text{MF}}$, the RG becomes unstable and diverges (region III in Fig. 1). It predicts an infinite number of phonons in the polaron cloud, caused by quantum fluctuations of the mobile impurity which are not included in MF theory. The effective polaron mass is also predicted to be infinite here. The range of scattering lengths within this regime depends on the impurity mass and the instability disappears for a localized impurity. To describe correctly the physics for these scattering lengths, phonon-phonon interactions neglected in our approach should be included. They are required to establish a balance with the attractive interactions between phonons induced by the mobile impurity.

The coherent state analysis of Ref. Shchadilova2016PRL revealed a regime in the spectrum where the eigenstates correspond to mesoscopic bound states, on the repulsive side close to the Feshbach resonance, i.e. for $1/a_{\rm IB}\gtrsim 1/a_{\rm IB,+}$. It was predicted that the weakly bound molecule state existing in this regime can be populated by multiple Bogoliubov phonons at the same time, an effect related to the formation of superpolarons observed for Rydberg atoms in a BEC Schmidt2016 ; Schlagmuller2016 . The instability of the RG is located precisely in this regime, $1/a_{\rm IB}\gtrsim 1/a_{\rm IB,+}^{\rm RG}$, which also indicates the presence of mesoscopic bound states. In this case ($④$ in Fig. 1) we expect a broad and featureless spectrum from the coherent state analysis Shchadilova2016PRL , difficult to distinguish qualitatively from the spectrum in the region $③$ of Fig. 1.

(IV) Repulsive polarons.– On the repulsive side of the Feshbach resonance, for $1/a_{\rm IB}>1/a_{\rm IB,+}^{\text{MF}}>0$, we predict a repulsive polaron peak at positive energies $E_{0}>0$ ($⑤$ and $⑥$ in Fig. 1). Because molecular bound states exist at lower energies (not shown in the spectrum), the polaron peak broadens and acquires a finite lifetime Rath2013 ; Li2014 ; Shchadilova2016PRL . This effect is particularly pronounced close to the resonance ($⑤$ in Fig. 1).

For weaker interactions, $1/a_{\rm IB}\gg 1/\xi$, the quasiparticle weight $Z$ of the repulsive polaron approaches unity, associated with a corresponding loss of weight in the incoherent part of the spectrum. For sufficiently weak couplings ($⑥$ in Fig. 1) the Fröhlich model can be used to describe the repulsive polaron.

(V) Efimov states.– In Refs. Levinsen2015 ; Sun2017 the possibility of Efimov bound states existing in the presence of the BEC has been discussed. Such states can not be captured by our description of the Bose polaron. In the spectral function shown in Fig. 1, they were predicted to appear first on the attractive side Levinsen2015 ; Sun2017 and extend all the way to the repulsive side. Although we do not expect dramatic changes of the polaron branch in the spectrum due to Efimov states, we note that they could give rise to additional broadening of the polaron peak. They, too, should become dressed by phonons from the surrounding BEC, but understanding the effect of strong dressing close to the Feshbach resonance remains an open problem.

Relation to experiments.– In experiments with ultracold atoms Catani2012 ; Spethmann2012 ; Fukuhara2013 ; Scelle2013 ; Hohmann2015 ; Jorgensen2016PRL ; Hu2016PRL ; Rentrop2016 both polaron branches (attractive and repulsive) can be explored. Recently the first observations of strong coupling Bose polarons in radio-frequency absorption spectra have been reported in two independent experiments Jorgensen2016PRL ; Hu2016PRL . In Appendix A we compare the results of our RG calculations to the measured absorption spectra Jorgensen2016PRL ; Hu2016PRL . We find that our results are consistent with the experimental observations.

To distinguish different theoretical predictions for the polaron properties measurements are required which go beyond the spectral function. For strong couplings we expect the quasiparticle weight $Z$ of the polaron to be small, see Fig. 2, and therefore most of the spectral weight corresponds to excited states. Ultracold atoms also offer the possibility to study dynamics of polarons and explore states which are inaccessible in the spectrum because of their negligibly small quasiparticle residue by preparing them adiabatically.

For example, by studying the dynamics of polarons oscillating in a harmonic trap, their effective mass can be accurately measured Catani2012 . In Fig. 3 we show our prediction of the effective polaron mass for an impurity in a BEC, for parameters as in the experiment by Jørgensen et al. Jorgensen2016PRL . A similar behavior is obtained for parameters as in the experiment of Hu et al. Hu2016PRL . We expect that the large increase of the polaron mass close to the resonance can be observed experimentally, enabling a more detailed comparison of different theoretical descriptions of strongly coupled Bose polarons in a BEC. The effective mass can also be used to study polaronic self-trapping, for a discussion see Appendix B.

Relation to other theoretical approaches.– We emphasize that the goal of our theory is to provide an accurate description of the polaron branch in the spectrum. However the RG formulated for the polaron also provides indirect indications for the existence of bound states, in the form of divergencies of the coupling constants and imaginary parts of the polaron energy Schmidt2013 . Similar RG flows in the complex plane were introduced for the description of Efimov four-body bound states Schmidt2010 ; Moroz2010 . This connects our work to the T-matrix theory of Ref. Rath2013 and to variational treatments of the impurity-BEC problem Li2014 ; Levinsen2015 where molecular bound states were explicitly included in the analysis. In contrast to these approaches we start from a MF theory of the polaron Lee1953 ; BeiBing2009 ; Shashi2014RF ; Shchadilova2016PRL which allows for an infinite number of phonon excitations (see also Ref. Rath2013 ). In addition we include correlations between phonons non-perturbatively and solve the resulting RG flow equations fully self-consistently. As a consequence the predictions of our RG method in the strong coupling regime are vastly different from those of truncated basis methods including only a few quasiparticle excitations at a time.

Like in Refs. Rath2013 ; Li2014 ; Levinsen2015 ; Jorgensen2016PRL ; Shchadilova2016PRL , our discussion is based on the use of Bogoliubov theory, assuming non-interacting phonons to describe the weakly interacting BEC. We expect that the Bogoliubov approximation is justified for the description of the polaron branch, as long as only virtual phonons are loosely attached to the impurity. On the other hand, the instability of the RG at strong couplings is expected to disappear when phonon-phonon interactions are included as in Refs. Ardila2015 ; blinova2013single . The properties of the bound states, which are not the main subject of this paper however, may depend more sensitively on residual phonon-phonon interactions.

An important difference between various theoretical approaches to polaron problems is the predicted dependence of the polaron energy on the large-momentum (UV) cut-off. In Ref. Grusdt2015RG a logarithmic divergence of the polaron energy with the UV cut-off was identified in the Fröhlich model, for which indications were found in diagrammatic quantum Monte Carlo calculations Vlietinck2015 . Using Gaussian variational wavefunctions, this divergence of the Fröhlich polaron energy was also confirmed Shchadilova2016 . Here we show analytically that the RG predicts the same logarithmic UV divergence as in the Fröhlich model for the full Bose polaron problem. It can be regularized by introducing a finite range interaction which effectively introduces a UV cut-off in the microscopic Hamiltonian, see also Ref. Tempere2009 . This log-divergence is a direct manifestation of quantum fluctuations in the polaron cloud which are not included by other approaches Rath2013 ; Li2014 ; Levinsen2015 ; Shchadilova2016PRL . Note that a closely related log-divergence was confirmed in a rigorous higher-order perturbative analysis Christensen2015 .

III Model

We consider an ultracold impurity atom of mass $M$ inside a weakly interacting, homogeneous BEC in three-dimensions. The condensate density will be denoted by $n_{0}$. The interaction of the impurity with the bosons (mass $m_{\rm B}$) in the BEC is modeled by a contact interaction characterized by the scattering length $a_{\text{IB}}$. The interaction between the bosons, characterized by the scattering length $a_{\rm BB}$, is assumed to be sufficiently weak so that Bogoliubov theory is valid for the description of the BEC and its excitations Pitaevskii2003 ; Pethick2008 .

Following Refs. Tempere2009 ; Grusdt2015Varenna ; Shchadilova2016PRL we obtain the following Hamiltonian ($\hbar=1$):

[TABLE]

While the first line correspond to the Bogoliubov-Fröhlich Hamiltonian Froehlich1954 , the second line describes processes involving two phonons, where bosons of finite momentum scatter off the impurity. Here we introduced Bogoliubov phonons $\hat{a}_{\bm{k}}$ by the relation

[TABLE]

where $\hat{d}_{\bm{k}}$ annihilates a boson at momentum $\bm{k}$. Their dispersion relation is given by $\omega_{k}=ck\sqrt{1+k^{2}\xi^{2}/2}$, with $\xi=1/\sqrt{2m_{\rm B}g_{\rm BB}n_{0}}$ and $c=\sqrt{g_{\rm BB}n_{0}/m_{\rm B}}$ being the healing length and speed of sound in the BEC, respectively. $g_{\rm BB}$ denotes the strength of the delta interactions between the bosons. In this paper we will mostly use $c/\xi$ and $\xi$ as units of energy and length, making the following relation useful: $m_{\rm B}=1/(\sqrt{2}c\xi)$. Moreover, the momentum and position operators of the impurity are defined by $\hat{\bm{p}}$ and $\hat{\bm{x}}$, and the Fröhlich type impurity-phonon interaction is described by $V_{k}=g_{\text{IB}}^{(0)}\sqrt{n_{0}}(2\pi)^{-3/2}W_{k}$, where

[TABLE]

The definition of $\theta_{k}$ can be found for example in Pitaevskii2003 ; Pethick2008 . Note that $W_{k}^{-1}=\cosh\theta_{k}+\sinh\theta_{k}$.

The microscopic interaction strength $g_{\text{IB}}^{(0)}$ is related to the scattering length $a_{\text{IB}}$ through the Lippmann-Schwinger equation. After introducing the UV cut-off $\Lambda_{0}$ for the regularization of the contact interaction one obtains

[TABLE]

where the reduced mass is defined by $m_{\rm red}^{-1}=M^{-1}+m_{\rm B}^{-1}$. Note that the right-hand side depends explicitly on the UV cut-off, and, as will become clear later, for $\Lambda_{0}\to\infty$ we need to keep the following terms in the expansion,

[TABLE]

The Hamiltonian (1) can be simplified by eliminating the impurity coordinate $\hat{\bm{x}}$. This is achieved by the Lee-Low-Pines (LLP) canonical transformation Lee1953 , which leads to the following Hamiltonian in the polaron frame depending explicitly on the conserved momentum $\bm{p}$,

[TABLE]

For more details on the LLP transformation see e.g. the reviews Devreese2013 ; Grusdt2015Varenna . In the following we will mostly be concerned with the spherically symmetric case $\bm{p}=0$, but a generalization of our RG to $\bm{p}\neq 0$ can be found in Appendix C.

IV Toy model: 0d polaron

Before we discuss the full model Eq. (1), let us gain some intuitive understanding of the polaron phase diagram. To this end we assume that the impurity is localized at $\bm{x}=0$ and has an infinite mass $M\to\infty$. Moreover we simplify the Hamiltonian by considering only a single $\bm{k}$ mode of Bogoliubov phonons, effectively making the system zero-dimensional. This leads us to the following toy model,

[TABLE]

We omitted factors of $2\pi$ and the $\bm{k}$ indices, set $V=g_{\text{IB}}\sqrt{n_{0}}W$ and simplified the two-phonon terms from the second line of Eq. (6). Moreover we dropped a constant energy shift $g_{\text{IB}}(W-W^{-1})^{2}/4$. We are not interested in a physical realization of this model, but rather it serves as a mathematical analog for understanding the different regimes in the spectral function, see Fig. 1.

MF theory.– We start by eliminating linear terms in $\hat{a}$ from the toy model. This is equivalent to a MF treatment of the problem and can be achieved by finding the saddle point of the classical Hamiltonian corresponding to Eq. (7). Formally we perform a unitary transformation

[TABLE]

and make the choice $\alpha^{\text{MF}}=-\beta_{\text{MF}}g_{\text{IB}}\sqrt{n_{0}}W/\omega$ with $\beta_{\text{MF}}=1/(1+g_{\text{IB}}W^{2}/\omega)$. As a result we obtain for the Hamiltonian in the new frame:

[TABLE]

Physically the saddle point corresponds to a MF polaron state, with the energy $E_{0}^{\text{MF}}=\beta_{\text{MF}}g_{\text{IB}}n_{0}$. This is not necessarily the lowest energy state however. The properties of the polaron, and in particular the spectral function $I(\omega)$, also depend on the states in the vicinity of the MF saddle point, which we study next.

Spectrum.– To calculate the excitation spectrum of $\tilde{\mathcal{H}}_{\rm toy}$, we introduce the following conjugate variables,

[TABLE]

with $[\delta\hat{n},\hat{\vartheta}]=i$, which are the analogues of the position and momentum operator for a harmonic oscillator. In our case of a BEC, they correspond to particle number and phase fluctuations. In terms of these new variables, the Hamiltonian becomes

[TABLE]

up to an overall energy offset, where

[TABLE]

Note that particle-number and phase fluctuations are decoupled in Eq. (12).

Now we would like to understand how the eigenenergies of our zero-dimensional toy model Hamiltonian change as the interaction strength $g_{\text{IB}}$ is varied. As in the original model we assume that $\omega>0$ is positive and $0<W\leq 1$. From Eq. (12) we recognize three different regimes with different signs of $A_{\pm}$, which we label by (I), (II) and (IV) in analogy with Fig. 1:

• (I)

When $A_{\pm}>0$ are both positive, the polaron corresponds to the true ground state in our toy model. The spectrum is bounded from below and there exists a discrete set of phonon excitations above the polaron energy.

• (II)

When $A_{+}>0$ and $A_{-}<0$, the polaron corresponds to a meta-stable saddle point. The spectrum of the toy model is dense and unbounded both from below and from above. The polaron is dynamically unstable in this case, with exponentially many phonons accumulating over time.

• (IV)

When $A_{\pm}<0$ are both negative, the polaron is the highest-energy state in our toy model. The spectrum is bounded from above and there exists a discrete set of excited states below the polaron energy. They can be understood as analogs of the molecular bound states found in Refs. Schmidt2016 ; Shchadilova2016PRL . Without additional terms in the Hamiltonian, the polaron is stable because it can not decay into the low-lying eigenstates.

The transitions between the three different regimes are determined by the conditions $A_{\pm}=0$. Therefore the critical interaction strengths are given by

[TABLE]

as sketched in Fig. 4 (a). While the polaron is dynamically stable for strong attraction (IV) where discrete bound states can form at lower energies, it is unstable for moderate attraction (II). In that case we expect phase separation because the condensate is depleted in time. This situation is analogous to a quantum particle with negative mass in a harmonic trapping potential, which only has a bounded trajectory if also the curvature of the potential is negative.

In Fig. 4 (b) we calculate the spectral function of the impurity for the toy model. It shows the same characteristic features as expected for the full three-dimensional polaron model, see Fig. 1 and Ref. Shchadilova2016PRL . In regions (I) and (IV) we observe a series of equidistant spectral lines, which correspond to multi-particle bound states with different numbers $n$ of phonons Schmidt2016 ; Schlagmuller2016 ; Shchadilova2016PRL . Their energies $E_{n}=n\Omega+\beta_{\text{MF}}g_{\text{IB}}n_{0}$ can be determined by diagonalization of the Hamiltonian (9) using a Bogoliubov transformation; $\tilde{\mathcal{H}}_{\rm toy}=\Omega~{}\hat{b}^{\dagger}\hat{b}$ up to an energy shift, and the binding energy is given by $\Omega=\sqrt{A_{+}A_{-}}$. Note that $\Omega$ is imaginary in the dynamically unstable regime (II).

In the following section we repeat our analysis of the MF polaron saddle point for the full three-dimensional problem in Eq. (6). We will show that the same three regimes as identified for the toy model exist in that case. Although we include an infinite number of phonon modes in the higher dimensional case, we will argue that still only a single mode becomes dynamically unstable and acquires negative energies for sufficiently strong attraction. It corresponds to a bound state of phonons to the impurity.

V MF theory of impurities in a BEC

In the polaron frame, i.e. after applying the LLP transformation, a variational MF theory of the polaron can be constructed by assuming a product wavefunction of coherent states

[TABLE]

where the coherent displacement operator is defined by

[TABLE]

This ansatz wavefunction was originally suggested by Lee, Low and Pines Lee1953 and has also been applied to the Fröhlich model in a BEC BeiBing2009 ; Shashi2014RF .

The MF theory for the Bose polaron in a BEC, as described by the Hamiltonian (1), has been formulated in Ref. Shchadilova2016PRL . Because it represents the starting point for our RG approach as well as for our analysis of the MF saddle point, we now provide a brief review of the key results relevant to our discussion. Afterwards we study the polaron phase diagram by extending the MF analysis to include molecular bound states around the saddle point.

V.1 Review of the MF polaron state

Starting from the MF variational energy $\mathscr{H}[\alpha_{\bm{k}}]=\langle 0|\hat{U}_{\text{MF}}^{\dagger}[\alpha_{\bm{k}}]\hat{\mathcal{H}}_{\rm LLP}\hat{U}_{\text{MF}}[\alpha_{\bm{k}}]|0\rangle$, the saddle point equations

[TABLE]

can be used to find a local minimum of the polaron energy. Their unique solution yields the coherent amplitudes corresponding to the MF Bose polaron Shchadilova2016PRL ; Shashi2014RF ,

[TABLE]

where we used spherical symmetry, $\alpha_{k}\equiv\alpha_{\bm{k}}$. The scattering amplitude $V_{k}$ is renormalized by a factor

[TABLE]

which is absent in the Fröhlich model Shashi2014RF where $\beta_{\text{MF}}=1$. The MF variational energy $E_{0}^{\text{MF}}=\mathscr{H}[\alpha_{\bm{k}}^{\text{MF}}]$ is given by

[TABLE]

Note that $g_{\text{IB}}^{(0)}$ needs to be expressed in terms of the scattering length, see Eq. (4), and the second term in Eq. (20) vanishes in the limit $\Lambda_{0}\to\infty$. Therefore we will discard it in our calculations. While $\beta_{\text{MF}}\to\infty$ is UV divergent when $\Lambda_{0}\to\infty$, the product $\beta_{\text{MF}}g_{\text{IB}}^{(0)}$ becomes UV convergent after introducing the scattering length $a_{\text{IB}}$,

[TABLE]

where

[TABLE]

This result is remarkable because it predicts a mere shift of the resonance, $1/a_{\text{IB}}\to 1/a_{\text{IB}}+1/a_{\text{IB},+}^{\text{MF}}$, relative to the weak coupling result $n_{0}2\pi a_{\text{IB}}/m_{\rm red}$.

Sometimes in the calculations below we find it convenient to retain the combination $\beta_{\text{MF}}g_{\text{IB}}^{(0)}$ in formal expressions. It should be kept in mind that $\beta_{\text{MF}}g_{\text{IB}}^{(0)}$ is UV convergent and can be expressed solely in terms of the scattering length $a_{\text{IB}}$, see Eq. (21).

V.2 Polaron phase diagram and MF theory of molecule formation

Now we extend our MF theory and analyze the nature of the saddle point (17) corresponding to the Bose polaron. To this end we study low-energy excited states around the MF polaron state. In the vicinity of the MF saddle point, $\alpha_{\bm{k}}=\alpha_{\bm{k}}^{\text{MF}}+\delta\alpha_{\bm{k}}$, the variational energy can be written as a quadratic form

[TABLE]

where $\delta\uline{\alpha}_{\bm{k}}=(\delta\alpha_{\bm{k}},\delta\alpha^{*}_{\bm{k}})^{T}$ and the Hessian matrix $\uuline{H}$ has entries $\uuline{H}_{1,1}(\bm{k},\bm{k}^{\prime})=\delta^{2}\mathscr{H}[\alpha_{\bm{k}}]/\delta\alpha_{\bm{k}}\delta\alpha_{\bm{k}^{\prime}}$, etc.; note that the dimension of $H$ is determined by the system size.

When the Hessian is positive definite (I in Fig. 5), the unique saddle point corresponds to the global minimum of $\mathscr{H}[\alpha_{\bm{k}}]$ and the MF polaron state is a good candidate for the ground state. If only one eigenvalue of $H$ becomes negative (II in Fig. 5), there exists a dynamically unstable mode Arnold1989 . In this case the MF polaron state becomes metastable and it is no longer a candidate for the ground state. Finally, if a pair of conjugate variables with negative eigenvalues exists (IV in Fig. 5), they describe dynamically stable modes at energies below the MF polaron energy $E_{0}^{\text{MF}}$. In this case, too, the MF polaron state is metastable and not a ground state candidate. These features are in direct analogy with our toy model from Sec. IV.

V.2.1 Quadratic Hamiltonian

We start by expanding the Hamiltonian (6) around the MF saddle point solution. To this end we express the two-phonon scattering processes from the second line of Eq. (6) in terms of

[TABLE]

and obtain an exact expression for $\tilde{\mathcal{H}}=\hat{U}^{\dagger}_{\text{MF}}\hat{\mathcal{H}}_{\rm LLP}\hat{U}_{\text{MF}}$,

[TABLE]

Here $:...:$ denotes normal ordering of the Bogoliubov phonons, and to simplify notations we introduced the operator

[TABLE]

The coupling constants $G_{\pm}=\frac{1}{2}\frac{g_{\text{IB}}^{(0)}}{(2\pi)^{3}}$ are equal, although later we find that they acquire different RG flows.

To understand the physical meaning of the new operators $\delta\hat{n}_{\bm{k}}$ and $\hat{\vartheta}_{\bm{k}}$ we note that, to leading order, the boson field operator $\hat{\psi}(\bm{r})$ can be expressed as

[TABLE]

where $\delta\hat{n}(\bm{r})$ ($\hat{\vartheta}(\bm{r})$) is the Fourier transform of $\delta\hat{n}_{\bm{k}}$ ($\hat{\vartheta}_{\bm{k}}$). Therefore they correspond to particle number and phase fluctuations of the BEC, respectively. Note that $\delta\hat{n}(\bm{r})$ and $\hat{\vartheta}(\bm{r})$ represent a pair of conjugate fields,

[TABLE]

When we consider only spherically symmetric perturbations $\delta\alpha_{k}$ (we assume $\bm{p}=0$), the non-linear term $\sim 1/M$ in the first line of Eq. (26) vanishes in the expression (23) for the variational energy. In this case we end up with a purely quadratic Hamiltonian, from which we can directly read off the Hessian matrix required for our analysis of the MF saddle point. Note that particle-number and phase fluctuations decouple, allowing us to treat these two conjugate variables separately now.

V.2.2 Phase fluctuations

First we consider perturbations $\delta\alpha_{k}=i\beta_{k}^{-}\in i\mathbb{R}$ describing only phase fluctuations, $\langle\hat{\vartheta}_{\bm{k}}\rangle=\beta_{k}^{-}W_{k}^{-1}/\sqrt{n_{0}}\equiv\vartheta_{k}$ and $\langle\delta\hat{n}_{\bm{k}}\rangle=0$. To find out whether the resulting quadratic energy functional $\mathscr{H}[\vartheta_{k}]=E_{0}^{\text{MF}}+\frac{1}{2}\int d^{3}kd^{3}k^{\prime}~{}\vartheta_{k}H_{\vartheta}(k,k^{\prime})\vartheta_{k^{\prime}}$ is positive definite, we search for a minimum $E_{0}^{\vartheta}$ of $\mathscr{H}[\vartheta_{k}]$ on a hyper-sphere around the MF saddle point defined by the condition

[TABLE]

$H_{\vartheta}$ is positive definite only if the minimum energy $E_{0}^{\vartheta}>E_{0}^{\text{MF}}$ is larger than the MF polaron energy; Otherwise the MF polaron state is metastable. In this case the same analysis needs to be repeated for the conjugate number fluctuations to check of which type the instability is.

Condition (30) can be taken into account in the minimization of $\mathscr{H}[\vartheta_{k}]$ by introducing a Lagrange multiplier $\mu_{-}$. Because the hyper-sphere is bounded, its global minimum is also a local minimum and we have to solve

[TABLE]

This equation has the solution

[TABLE]

where $I_{-}=\left[\int d^{3}k~{}4G_{-}^{2}W_{k}^{-2}(\Omega_{k}^{\text{MF}}-\mu_{-})^{-2}\right]^{-1/2}$ is a real-valued normalization constant; The Lagrange multiplier is determined by the following equation,

[TABLE]

and the variational energy at the local extremum is

[TABLE]

Eq. (33) only has solutions for $\mu_{-}<0$, because otherwise the integral on the right hand side diverges. According to Eq. (34) all such solutions describe modes with energies below the polaron saddle point energy. Using Eq. (33) again we find that solutions $\mu_{-}<0$ only exist if

[TABLE]

which is equivalent to the condition that

[TABLE]

Therefore at $a_{\text{IB},-}^{\text{MF}}$ the MF polaron state becomes unstable, giving rise to broadening of the polaron peak in the spectrum as discussed in connection with Fig. 1. In the regime where $1/a_{\text{IB}}>1/a_{\text{IB},-}^{\text{MF}}$ there always exist states with energies arbitrarily deep below that of the MF polaron state, $E_{0}^{\text{MF}}$. Because $0<W_{k}<1$, and by comparison with Eq. (22), we find that

[TABLE]

i.e. the instability is located on the attractive side of the Feshbach resonance.

V.2.3 Particle-number fluctuations

Now we repeat our analysis and study the stability of the MF polaron state with respect to particle-number fluctuations. We consider real-valued variations $\delta\alpha_{k}=\beta_{k}^{+}\in\mathbb{R}$ for which $\langle\delta\hat{n}_{\bm{k}}\rangle=2\beta_{k}^{+}W_{k}\sqrt{n_{0}}\equiv\delta n_{k}$ and $\langle\hat{\vartheta}_{\bm{k}}\rangle=0$. Similar to the case of phase fluctuations we minimize the variational energy $\mathscr{H}[\delta n_{k}]$ on the hyper-sphere defined by Eq. (30).

After introducing the Lagrange multiplier $\mu_{+}$ we obtain the solution

[TABLE]

where $I_{+}=\left[\int d^{3}k~{}4G_{+}^{2}W_{k}^{2}(\Omega_{k}^{\text{MF}}-\mu_{+})^{-2}\right]^{-1/2}$ is a normalization constant. The Lagrange multiplier $\mu_{+}$ is determined by

[TABLE]

which only has a solution for $\mu_{+}<0$ and requires

[TABLE]

As in the case of phase fluctuations the variational energy at the local extremum is given by

[TABLE]

Thus we conclude that the variational energy functional $\mathscr{H}[\delta n_{k}]$ is positive definite for $1/a_{\text{IB}}<1/a_{\text{IB},+}^{\text{MF}}$, see Fig. 5. The MF polaron state becomes unstable with respect to particle-number fluctuations only on the repulsive side of the Feshbach resonance.

V.2.4 Molecule formation in a BEC

Combining our results for particle-number and phase fluctuations around the MF polaron state, we identify three characteristic regimes in the polaron phase diagram, see Fig. 5. On the attractive side of the Feshbach resonance, I in Fig. 5, the MF polaron state is stable and represents a good ground state candidate. At $a_{\text{IB},-}^{\text{MF}}$ the polaron becomes unstable with regard to phase fluctuations, while the conjugate particle-number fluctuations remain stable. In this regime, II in Fig. 5, there exists a dynamically unstable mode at energies below the MF polaron energy.

At $a_{\text{IB},+}^{\text{MF}}$, in addition, the polaron becomes unstable with respect to particle-number fluctuations. In this regime, IV in Fig. 5, the polaron is unstable but there exists a mode at an energy below that of the MF polaron state. This state corresponds to a molecular bound state, as can be seen by considering the two-particle case where $n_{0}\to 0$. The healing length $\xi\to\infty$ diverges in this limit and we find that

[TABLE]

corresponding to unitarity. Indeed, here a weakly bound Feshbach molecule appears at the continuum threshold Chin2010 whose wavefunction in momentum space is correctly described by Eqs.(32), (38).

So far we only discussed the properties of the classical variational energy landscape $\mathscr{H}[\alpha_{\bm{k}}]$ but we ignored that $\hat{\vartheta}_{\bm{k}}$ and $\delta\hat{n}_{\bm{k}}$ are actually non-commuting conjugate variables. In Appendix D we take into account quantum fluctuations and construct a quantum theory of the negative energy modes which we have identified above.

VI Renormalization group theory

Within the MF ansatz (15) the phonons in the polaron cloud are treated as uncorrelated excitations. The goal of this paper is to remedy this shortcoming. To this end we will go beyond the MF description summarized in the last section and include correlations between phonons in the polaron cloud. Our starting point is the Hamiltonian

[TABLE]

which describes quantum fluctuations around the MF state. Our method is a generalization of the RG approach developed for the Fröhlich Hamiltonian in Refs. Grusdt2015RG ; Grusdt2015Varenna ; Grusdt2016RG . Unlike more traditional formulations of the RG Altland2010 , we use only wavefunctions and unitary transformations and there is no need to employ path integrals.

Our approach is very similar in spirit to Wegner’s flow equation method Wegner1994 . In both approaches the basic idea is to subsequently diagonalize the Hamiltonian by applying a series of infinitesimal unitary transformations. We make use of the separation of energy scales of phonons at different momenta, i.e. $\omega_{k_{1}}\gg\omega_{k_{2}}$ for $k_{1}\gg k_{2}$, to formulate a momentum-shell RG for an effective Hamiltonian $\tilde{\mathcal{H}}(\Lambda)$. Here $\Lambda<\Lambda_{0}$ defines a high momentum cut-off beyond which fluctuations have already been accounted for. While Wegner’s approach is more general because it does not rely on the existence of a separation of energy scales, it produces significantly more complicated flow equations (at least for the polaron problem at hand).

The derivation of the RG flow equations is closely related to the treatments presented in Refs. Grusdt2015RG ; Grusdt2015Varenna ; Grusdt2016RG . Therefore we will not provide a description of our method in the main text but we delegated an overview as well as detailed calculations to Appendix E. The remainder of this section is organized as follows. In subsection VI.1 we introduce the Hamiltonian used in the RG and explain the meaning of the coupling constants. Their flow equations are presented in part VI.2. In subsection VI.3 we summarize how the flow of the coupling constants gives rise to the polaron phase diagram. We proceed by a dimensional analysis in subsection VI.4. When resonant three-body collisions are possible the RG has divergencies in the polaron energy which we regularize in subsection VI.5. We also discuss physical implications of these divergencies and draw analogies with the well-understood case of two-particle scattering.

VI.1 Universal Hamiltonian and coupling constants

The universal form of the Hamiltonian in the spherically symmetric case (i.e. for zero total polaron momentum, $\bm{p}=0$) is given by

[TABLE]

The first line corresponds to the Fröhlich polaron case Grusdt2015RG . Here $\Omega_{k}$ can be interpreted as the renormalized dispersion of the phonons in the polaron frame. The second term describes phonon-phonon interactions induced by the impurity, where $\mathcal{M}(\Lambda)\geq M$ is the renormalized impurity mass flowing in the RG. We defined the operator $\hat{\Gamma}_{\bm{k}}$ as in Eq. (27),

[TABLE]

but replaced the coherent MF amplitude $\alpha_{\bm{k}}^{\text{MF}}$ by $\alpha_{k}(\Lambda)$ which is also flowing in the RG Grusdt2016RG ,

[TABLE]

The second line of Eq. (44) describes two-phonon terms beyond the Fröhlich model, see also Eq. (26).

Summarizing, the coupling constants of the RG which depend on the cut-off $\Lambda$ are the renormalized impurity mass $\mathcal{M}$, the two interaction strengths $G_{\pm}$ of particle number and phase fluctuations, as well as the renormalized MF amplitude $\beta$. The energy $E_{0}(\Lambda)$ also flows in the RG procedure and approaches the polaron energy when $\Lambda\to 0$. The associated initial conditions are

[TABLE]

Note that the Hamiltonian $\tilde{\mathcal{H}}(\Lambda)$ in Eq. (44) is exact for $\Lambda=\Lambda_{0}$, while for smaller cut-offs $\Lambda<\Lambda_{0}$ it keeps its form within our approximate RG procedure, see Appendix E. In particular, the coupling constants $G_{\pm}(\Lambda)$ only depend on $\Lambda$, related to an energy scale, but not on the momenta $\bm{k}$, $\bm{k}^{\prime}$ of particle-number and phase fluctuations. A more accurate description including also the momentum dependence was provided for the Fröhlich terms of the Hamiltonian in Ref. Enss2009 .

VI.2 RG flow equations

The RG flow equations for the coupling constants are derived in Appendix E by applying infinitesimal unitary transformations $\hat{U}_{\Lambda}$ to decouple fast from slow phonons. This is achieved by explicitly constructing $\hat{U}_{\Lambda}^{-1}$ perturbatively in $1/\Omega_{k}$. In the spherically symmetric case considered here (i.e. assuming $\bm{p}=0$) we obtain

[TABLE]

Here $\int_{\text{F}}d^{2}k=4\pi\Lambda^{2}$ denotes an integral over a fast phonon shell (F) evaluated at the momentum cut-off $k=\Lambda$.

The RG flow equation (49) of the mass $\mathcal{M}$ is identical to its counterpart in the Fröhlich model Grusdt2015RG . Note however that the MF amplitude $\alpha_{k}$ is renormalized as a consequence of the two-phonon terms, which modifies the resulting mass as compared to the Fröhlich case.

The flow of the groundstate energy is given by

[TABLE]

where $\int_{\text{S}}d^{3}p=\int_{0}^{\Lambda}dp~{}4\pi p^{2}$ denotes an integral over the slow phonons. The first term on the right hand side is identical to the expression derived for the Fröhlich model and describes the effect of quantum fluctuations of the impurity, see Ref. Grusdt2015RG . As in the case of the effective mass $\mathcal{M}$ the result is renormalized because the MF amplitude $\alpha_{k}$ is modified by two-phonon terms. The last term in Eq. (52) is new and entirely due to two-phonon scattering. Note that the RG flow of $E_{0}(\Lambda)$ depends on the MF amplitude $\alpha_{p}$ at lower energies, which is determined fully self-consistently in every RG step as described in Ref. Grusdt2016RG .

The flow equation for the coupling constant $\beta$ introduced in Eq. (46) above,

[TABLE]

can be solved exactly and we obtain the solution

[TABLE]

VI.3 Main features of the RG flows

In Fig. 6 we briefly summarize the main features of the RG flow equations (49) - (53) to simplify the following analysis. To this end we analyze for which cut-offs the coupling constants diverge. The connection between the bound states derived from a MF analysis and these divergencies is explained in Sec. VI.5, and a detailed discussion of the polaron phase diagram is provided in Sec. VIII.

For $1/a_{\rm IB}<1/a_{\rm IB,-}^{\rm RG}$ (region I in Fig. 6) there are no divergencies. At $a_{\rm IB,-}^{\rm RG}$ the coupling constant $G_{-}$ diverges in the low energy (IR) limit for $\Lambda\to 0$. Beyond this point, for $1/a_{\rm IB}>1/a_{\rm IB,-}^{\rm RG}$ (regions II-IV in Fig. 6), this divergence is shifted to finite cut-offs $\Lambda_{\rm c}^{-}>0$. We will argue in subsection VI.5 that this divergence can be regularized and gives rise to an imaginary contribution to the polaron energy.

Similarly, at $a_{\rm IB,+}^{\rm RG}$ the coupling constant $G_{+}$ diverges at $\Lambda=0$. At the same point the effective mass $\mathcal{M}$ diverges for $\Lambda\to 0$. For all $1/a_{\rm IB}>1/a_{\rm IB,+}^{\rm RG}$ (regions III, IV in Fig. 6) the coupling constants $G_{\pm}$ diverge at finite cut-offs $\Lambda_{\rm c}^{\pm}>0$ and both divergencies can be regularized, see VI.5.

In the regime between $1/a_{\rm IB,+}^{\rm RG}<1/a_{\rm IB}<1/a_{\rm IB,+}^{\rm MF}$ (region III in Fig. 6) the coupling constant $\beta(\Lambda)$ diverges at a finite cut-off $\Lambda_{\rm c}^{\beta}>0$. This divergence can not be regularized and the RG breaks down in this regime.

VI.4 Dimensional analysis

Now we perform a dimensional analysis to understand when the RG is reliable and in order to identify the most relevant terms in the Hamiltonian. We will show that the Fröhlich terms in the Hamiltonian provide a universal description of the polaron branch in a regime not too close to the critical points $a_{\text{IB},\pm}^{\rm RG}$ where $G_{\pm}$ diverge in the IR limit, respectively. We keep the discussion in this section general and consider arbitrary dimensionality $d$. For the Hamiltonian (44) studied in this paper, $d=3$.

The basic idea is to rescale phonon operators $\hat{a}_{\bm{k}}\sim\Lambda^{\gamma}$ in such a way that the free part $\int^{\Lambda}d^{d}k~{}\Omega_{k}\hat{a}^{\dagger}_{\bm{k}}\hat{a}_{\bm{k}}\simeq\Lambda^{0}=1$ of the universal Hamiltonian (44) becomes of order unity, for discussion see e.g.Altland2010 . Then we compare the scaling of the remaining interaction terms with the cut-off $\Lambda$. We distinguish between two regimes of our theory which we treat separately: The high energy regime where $k\gg 1/\xi$ and the low energy regime where $k\ll 1/\xi$.

Our results are summarized in Tabs. 1 and 2, where we used engineering dimensions of the operators. We provide the derivation of these results in Appendix F. Our analysis shows that in $d=3$ dimensions the two-phonon terms $\sim G_{\pm}$ extending the Fröhlich Hamiltonian are irrelevant. This suggests that the Fröhlich Hamiltonian provides a universal description of the Bose polaron for $d=3$. Quantum fluctuations $\sim\mathcal{M}^{-1}:\hat{\Gamma}_{\bm{k}}\hat{\Gamma}_{\bm{k}^{\prime}}:$ around the Fröhlich MF polaron state become relevant in the regime where $k\gg 1/\xi$. For $k\ll 1/\xi$ they are marginal, and when the polaron mass $\mathcal{M}$ is sufficiently large we expect that MF polaron theory is valid. Therefore, as in the Fröhlich model, the most relevant phonons are those with momenta $k\simeq 1/\xi$.

Our simple analysis above suggests that in $d=3$ dimensions the two-phonon terms can not modify the fixed point of the Fröhlich Hamiltonian where $\mathcal{M}\to\infty$ diverges. Note however that $\hat{\Gamma}_{\bm{k}}$ depends explicitly on $G_{+}$, see Eqs.(45), (46) and (56). Therefore the RG flow of $G_{+}$ determines the position of the strong coupling fixed point of the Fröhlich Hamiltonian. Indeed, as shown in Fig. 6, the divergencies of $G_{+}$ and $\mathcal{M}$ are located at the same scattering length $a_{\rm IB,+}^{\rm RG}$. This is not in contradiction with the fact that the $G_{+}$ term appears to be irrelevant when engineering dimensions are used as in Tabs. 1 and 2.

In Appendix F we also include the scalings of the coupling constants with the cut-off $\Lambda$ in the analysis, i.e. we consider the anomalous scaling dimensions of the different terms. For $d=3$ dimensions the same conclusions can be drawn as from engineering dimensions, provided that $a_{\rm IB}$ is sufficiently far from $a_{\rm IB,+}^{\rm RG}$ where $G_{+}$ diverges. However, we also show in Appendix F that the anomalous scaling dimension of $G_{+}$ changes close to $a_{\rm IB,+}^{\rm RG}$. As a result, two-phonon scattering can become marginal close to $a_{\rm IB,+}^{\rm RG}$. In this regime a competition between two-phonon and Fröhlich type terms in the Hamiltonian is expected.

The second coupling constant $G_{-}$ does not enter the effective Fröhlich Hamiltonian. Because it is irrelevant according to Tabs. 1 and 2, we expect that it has only a small effect on the polaron wavefunction. As will be shown later however, the divergence of $G_{-}$ leads to a change of the ground state of the system. As a result, terms which are irrelevant in the RG can lead to couplings of the polaron to the new ground state, giving rise to a finite lifetime of the polaron and thus to broadening of the spectral function as discussed in the introduction. This strong influence of irrelevant terms in the Hamiltonian on qualitative properties of the polaron state is analogous to the more general concept of dangerously irrelevant operators Amit1982 .

As in the case of $G_{+}$, the two-phonon term with strength $G_{-}$ becomes marginal close to a divergence of $G_{-}$ when its anomalous scaling dimension is taken into account. For details see Appendix F.

VI.5 Regularization of $G_{\pm}$ divergencies

For some scattering lengths the coupling constants $G_{\pm}\to\infty$ diverge during the RG flow, as summarized in Fig. 6. When $G_{\pm}(\Lambda)$ diverge in the IR limit for $\Lambda\to 0$, this indicates a transition where the nature of the ground state changes abruptly. We will now examine these divergencies more closely and describe their physical meaning. To this end we relate them to the appearance of modes below the polaron energy which we discussed in Secs. IV, V.2 using MF theory. We will start from the well-understood case of two-particle scattering where $n_{0}=0$.

While the polaron wavefunction flowing in the RG remains well-defined even when $G_{\pm}$ diverge at some point during the RG, the corresponding polaron energy diverges. We will now show that this divergence of the polaron energy can be regularized by introducing an infinitesimal imaginary part to $G_{\pm}$. This allows us to make predictions for the polaron properties despite the fact that $G_{\pm}$ may be divergent during the RG flow.

VI.5.1 The two-particle case

We start by revisiting the two-particle scattering problem where $n_{0}=0$ and the Fröhlich part of the Hamiltonian becomes trivial. In this case $W_{k}=1$, $\Omega_{k}=k^{2}/2m_{\rm red}$ and the RG flows of $G_{+}(\Lambda)=G_{-}(\Lambda)=g_{\text{IB}}(\Lambda)/2(2\pi)^{3}$ coincide. They reproduce the familiar result of the Lippmann-Schwinger equation:

[TABLE]

Let us consider $a_{\text{IB}}>0$ corresponding to the repulsive side of the resonance, such that $g_{\text{IB}}(0)>0$. Because $g_{\text{IB}}^{(0)}<0$ (for sufficiently large $\Lambda_{0}$, as considered here), the interaction strength $g_{\text{IB}}(\Lambda)$ diverges at some cut-off $\Lambda_{c}\geq 0$ during the RG flow (i.e. $g_{\text{IB}}^{-1}(\Lambda_{c})=0$). Nevertheless the Lippmann-Schwinger equation is exactly reproduced and the divergence of the coupling constants $G_{\pm}\to\infty$ does not invalidate the RG procedure, despite the fact that the two-particle scattering terms are only treated to lowest order in $G_{\pm}$ in every RG step.

The divergence of $g_{\text{IB}}(0)\to\infty$ when $a_{\text{IB}}\to\infty$ is related to the appearance of a two-particle bound state on the repulsive side of the Feshbach resonance Chin2010 . Its binding energy can be calculated e.g. from MF theory, see Eq. (106), or from a calculation of the full T-matrix using a two-body RG formalism Mitra1998 . Instead of following these exact approaches, we will now use a simple argument to estimate the bound state energy $E_{b}$ and connect it to the divergence of the interaction strength $g_{\text{IB}}(\Lambda)$ in the RG.

When $a_{\text{IB}}=\infty$ the divergence of $g_{\text{IB}}(\Lambda)$ is located in the deep IR limit at $\Lambda_{c}=0$. For $0<a_{\text{IB}}<\infty$ on the other hand, the divergence of $g_{\text{IB}}(\Lambda_{c})\to\infty$ appears at a finite momentum cut-off $\Lambda_{c}$, which can be associated with a characteristic energy of the bound state, $\Omega_{\Lambda_{c}}\approx E_{b}$. From this equation, and using $g_{\text{IB}}^{-1}(\Lambda_{c})=0$, we obtain the following estimate for the binding energy

[TABLE]

Up to a numerical factor $\pi^{2}/4$ this result yields the correct expression for the universal binding energy of the dimer $E_{\rm dim}=1/2m_{\rm red}a_{\text{IB}}^{2}$. We can understand this factor by recalling the simplifications used in the RG and using the exact equation for the dimer energy,

[TABLE]

The RG procedure is perturbative in $1/\Omega_{k}=2m_{\rm red}/k^{2}$ and at the cut-off $\Lambda_{c}$ only bosons at momenta $k>\Lambda_{c}$ are taken into account. If we evaluate the last expression in this regime and use perturbation theory in $\Omega_{k}^{-1}$ we arrive at the following approximation:

[TABLE]

For $\Lambda_{0}\to\infty$ this yields the same result as the RG, namely $E_{\rm dim}\approx\pi^{2}/8m_{\rm red}a_{\text{IB}}^{2}$. This indicates that the RG gives the correct physics.

In the two-particle case, the ground state energy $E_{0}(\Lambda)$ is not flowing at all. The Fröhlich term in Eq. (52) vanishes because $n_{0}=0$ and the two-particle term vanishes because $G_{+}=G_{-}$ and $W_{k}=1$.

VI.5.2 The many-body case

In the many-body case we can also interpret the divergencies of $G_{\pm}$ as indicators for bound state formation of phonons with the impurity. In contrast to the two-particle case we obtain two separate divergencies associated with particle number ($G_{+}\to\infty$) and phase fluctuation scattering ($G_{-}\to\infty$) respectively. These two divergencies correspond to the instabilities of the MF polaron state with respect to particle-number and phase fluctuations, see our discussion in Sec. V.2.

Unlike in the two-particle case, the divergencies of $G_{+}\neq G_{-}$ affect the polaron energy $E_{0}$. From Eq. (52) we can see that corrections due to two-phonon terms diverge when $G_{\pm}^{-1}(\Lambda_{c}^{\pm})=0$, i.e. $E_{0}\to\infty$. Now we introduce a regularization scheme to deal with these divergencies. Their origin are resonant three-particle processes, where bosons recombine with the impurity to form a bound state. Unlike in the two-particle case such processes are explicitly allowed because the total phonon number is no longer conserved in the presence of the condensate.

To regularize the polaron energy we include an imaginary part in $G_{\pm}$ as suggested in Refs. Schmidt2010 ; Moroz2010 ,

[TABLE]

We start by discussing the case when only $G_{-}$ diverges but $G_{+}$ remains finite.

The contribution to the ground state energy from phase fluctuation scattering on the impurity becomes

[TABLE]

Now we use $-\partial_{\Lambda}(G_{-}^{-1}+i\epsilon)^{-1}=(G_{-}^{-1}+i\epsilon)^{-2}(\partial_{\Lambda}G_{-}^{-1})$ for partial integration. In combination with the RG flow equation (51) and the identity $(G_{-}^{-1}+i\epsilon)^{-1}=\mathcal{P}\frac{1}{G_{-}^{-1}}-i\pi\delta(G_{-}^{-1})$ valid for $\epsilon\to 0^{+}$, we obtain for $G_{-}^{-1}(\Lambda_{c}^{-})=0$:

[TABLE]

The principle-value integral in the first line is finite when $G_{-}$ diverges at $\Lambda_{c}^{-}$, and we obtain a well-defined energy. It has an imaginary part which can be interpret as a decay rate into the bound state. In this case we expect a broadened polaron peak in the spectral function, and additional weight appearing below.

When both couplings $G_{\pm}$ diverge at cut-offs $\Lambda_{c}^{\pm}$ during the RG, we can use a similar regularization scheme to obtain a finite energy with an imaginary part.

VI.6 UV log-divergence of the energy

Now we study how the polaron energy $E_{0}$ depends on the UV cut-off $\Lambda_{0}$. For Fröhlich polarons in $d=3$ dimensions we have shown that $E_{0}\simeq-\log(\Lambda_{0}\xi)$ diverges logarithmically with the UV cut-off $\Lambda_{0}$ Grusdt2015RG . Because the Fröhlich Hamiltonian provides a universal description of the polaron branch (see Sec. VI.4) we expect the same log-divergence for the Bose polaron with zero range contact interactions. Indeed we will show by an explicit calculation that the ground state energy diverges logarithmically with the UV cut-off,

[TABLE]

Here $\Lambda_{1}$ denotes a momentum scale at low energies which is given by $\Lambda_{1}\approx\xi^{-1}$ if the total mass renormalization is small, $\mathcal{M}-M\ll M$. If the total mass renormalization is large, $\mathcal{M}\gg M$, it is given by $\Lambda_{1}\approx n_{0}a_{\text{IB}}^{2}m_{\rm red}/M$.

The UV log-divergence is a direct consequence of the polaron mass renormalization caused by quantum fluctuations of the mobile impurity, for a discussion see Refs. Grusdt2015RG ; Grusdt2015Varenna . It was confirmed by diagrammatic Monte Carlo calculations for the Fröhlich model Vlietinck2015 . To regularize the polaron energy, we suggest taking into account the finite range $r_{e}\neq 0$ of the impurity-boson interactions. Effectively this leads to a UV cut-off $\Lambda_{0}\approx 1/r_{e}$ which is of the order of $\Lambda_{0}\approx 200/\xi$ (in units of $\xi$) for typical experimental parameters, see Refs. Tempere2009 ; Jorgensen2016PRL .

We believe that the log-divergence provides an important test-case for theories of the Bose polaron. Simpler approaches like MF theory yield UV convergent polaron energies and thus do not fully include the quantum fluctuations leading to the log-divergence. In the diffusive Monte Carlo calculations of Ref. Ardila2015 finite range interactions have been used, but the dependence of the energy on the effective range has not been studied in detail. In the self-consistent T-matrix calculations of Ref. Rath2013 zero-range interactions were considered, but the energy was UV convergent. This implies that the quantum fluctuations included in the self-consistent T-matrix approach are different from the ones captured by the RG.

Now we derive the log-divergence of the polaron energy from the Fröhlich term. We show in Appendix G that two-phonon terms only lead to a UV-convergent correction to the polaron energy. In Appendix H we discuss how the log-divergence is connected to similar logarithmic corrections predicted by different theoretical methods in closely related physical systems.

VI.6.1 Derivation of the Fröhlich type log-divergence

We start by solving the RG flow of the mass $\mathcal{M}$ in the UV limit. There, to leading order in $\Lambda^{-1}$, it holds $\Omega_{\Lambda}\simeq\Lambda^{2}/2m_{\rm red}$ with the bare reduced mass $m_{\rm red}=1/\left(m_{\rm B}^{-1}+M^{-1}\right)$, and $g_{\text{IB}}^{(0)}\beta(\Lambda)\simeq g_{\text{IB}}^{(0)}\beta_{\text{MF}}$. Hence

[TABLE]

Integrating this separable differential equation yields the UV behavior,

[TABLE]

showing that deviations of $\mathcal{M}$ from $M$ correspond to higher order terms in $\Lambda^{-1}$.

In the UV limit we further obtain,

[TABLE]

and upon integration we find that the first term in Eq. (52) leads to the UV divergence (64). The only difference to the result previously encountered in the Fröhlich model is the renormalization of the pre-factor $g_{\text{IB}}^{(0)}\to\beta_{\text{MF}}g_{\text{IB}}^{(0)}$ because we expand around the MF solution of the full model.

To understand how the log-divergence is cut-off in the IR limit at $\Lambda_{1}$, we consider two separate cases. When impurity-boson interactions are sufficiently weak the mass renormalization is small for all cut-offs, $\mathcal{M}(\Lambda)-M\ll M$. Therefore $\partial_{\Lambda}\mathcal{M}^{-1}\approx-M^{-2}\partial_{\Lambda}\mathcal{M}\simeq\Lambda^{-2}$ until the RG flow of $\mathcal{M}^{-1}$ stops around $\Lambda_{1}=1/\xi$. This cuts off the log-divergence at

[TABLE]

In the second case impurity-boson interactions are so strong that the mass renormalization becomes dramatic, $\mathcal{M}\gg M$. This happens when $\mathcal{M}(\Lambda_{1})-M\approx M$, which yields

[TABLE]

i.e. $\Lambda_{1}\sim n_{0}a_{\text{IB}}^{2}m_{\rm red}/M$, using Eq. (66) for large $\Lambda_{0}\to\infty$. For smaller cut-offs $\Lambda<\Lambda_{1}$ we find $\partial_{\Lambda}\mathcal{M}^{-1}=-\mathcal{M}^{-2}\partial_{\Lambda}\mathcal{M}\simeq\Lambda^{2}\Lambda^{-2}=1$. Hence the logarithmic flow of the ground state energy ends, $\partial_{\Lambda}E_{0}\simeq\Lambda$ for $1/\xi\leq\Lambda\lesssim\Lambda_{1}$, see Eq. (52).

VII Results: Polaron energies

In this section we use the RG approach for the calculation of the polaron energy $E_{0}$ using Eq. (52). In the regime I, where $1/a_{\text{IB}}<1/a_{\text{IB},-}^{\rm RG}$ and before $G_{\pm}\to\infty$ diverge, we expect no low lying states in the spectrum below $E_{0}$. Within our theory the polaron energy corresponds to the ground state energy in this regime. Note however that we did not include Efimov states Levinsen2015 ; Sun2017 which may have even lower energies in the case of short-range interactions Braaten2006 . On the other hand, when $G_{\pm}(\Lambda)\to\infty$ diverge during the RG, we expect additional states at lower energies (see also Ref. Shchadilova2016PRL ) and $E_{0}$ corresponds to the energy of the polaron branch in the spectrum. In this regime we make use of the regularization scheme for the polaron energy discussed in Sec. VI.5.

In Fig. 7 we show the polaron energy as a function of $a_{\text{IB}}^{-1}$ (solid line) and find an attractive and a repulsive polaron branch. Close to the two-particle Feshbach resonance (at $a_{\text{IB}}^{-1}=0$) we predict large deviations from MF theory due to quantum fluctuations.

On the repulsive side of the resonance, the polaron energy saturates before the strong-coupling phase is reached where the RG breaks down. Similar saturation of the polaron properties in this regime has been observed in Figs. 2, 3. Such behavior would be expected for bubble polarons which have been predicted in this regime by Gross-Pitaevskii MF calculations blinova2013single .

We have argued by a dimensional analysis in VI.4 that the Fröhlich terms in the Hamiltonian are most relevant in the RG. Sufficiently far from divergencies of $G_{\pm}$ we expect that the RG flows towards an effective Fröhlich Hamiltonian with renormalized parameters. To demonstrate this, we considered a renormalized Fröhlich model where we ignored the additional contributions of two-phonon terms in the energy, i.e. the second line of Eq. (52). As can be seen in Fig. 7 the polaron energy predicted by this renormalized Fröhlich model is in excellent agreement with the full calculation based on Eq. (52) for the extended Fröhlich model.

To understand how much the Fröhlich model is renormalized by two-phonon scattering, we also compare our results to the original Fröhlich model introduced in Ref. Tempere2009 , for which we developed the RG in Ref. Grusdt2015RG . In this case two-phonon terms are completely neglected. On the attractive side we find reasonable quantiative agreement with the full calculation, and qualitatively the correct behavior is predicted. On the repulsive side the original Fröhlich model does not break down at strong couplings, while the curvature of $E_{0}(a_{\text{IB}}^{-1})$ is correctly captured. For sufficiently weak repulsion the agreement is reasonable.

In Fig. 7 we also compare our results to second order perturbation theory in $a_{\text{IB}}$ (dotted line). To this order the original Fröhlich model is exact, but it was shown in Ref. Christensen2015 that it misses higher order terms. This has caused serious concerns that the Fröhlich Hamiltonian can not be used to make predictions beyond those of second order perturbation theory. In Ref. Bruderer2007 a condition was derived when two-phonon terms can be safely neglected, see also Ref. Grusdt2015Varenna . For the parameters used in Fig. 7 it becomes

[TABLE]

In this regime all theories basically coincide, and beyond it the perturbative result deviates substantially, in particular on the attractive side. Remarkably, the original Fröhlich model provides a reasonable description of the polaron energy derived from the full model even beyond the careful estimate above.

VIII Polaron phase diagram

Now we discuss the phase diagram predicted by the RG, including quantum fluctuations around the MF saddle point solution. Their most dramatic effect is to cause divergencies of the interaction strengths $G_{\pm}$ which we related to the appearance of new states below the polaron energy in Secs. V.2, VI.5. This gives rise to the two critical values $a_{\text{IB},\pm}$ where the features in the spectral function change as discussed in Sec. II.

In Fig. 8 we calculate the IR values of the coupling constants $G_{\pm}(\Lambda=0)$ as a function of the inverse two-particle scattering length $1/a_{\text{IB}}$. Let us start our discussion far on the attractive side, $1/a_{\text{IB}}\ll-1/\xi$, where the ground state is an attractive Fröhlich polaron Tempere2009 ; Rath2013 ; Shashi2014RF . Here $G_{\pm}(\Lambda)<0$ during the entire RG flow (note that $G_{\pm}(\Lambda_{0})<0$ is always attractive in the UV).

VIII.1 Phase fluctuations

When we approach the Feshbach resonance the interaction strength $G_{-}$ associated with phase-fluctuation diverges first, on the attractive side of the two-particle resonance. Note that because $W_{k}<1$ it always holds $G_{-}<G_{+}$. From the exact solution of the RG flow equation,

[TABLE]

we obtain the divergence $G_{-}^{-1}(0)=0$ for $a_{\text{IB}}=a_{\text{IB},-}^{\rm RG}$, with

[TABLE]

The last inequality follows from $W_{\Lambda}\leq 1$.

Eq. (72) is the analogue of $a_{\text{IB},-}^{\text{MF}}$, see Eq. (36), corrected for the effects of quantum fluctuations included in the RG. Because $\mathcal{M}\geq M$ it follows that

[TABLE]

For $a_{\text{IB}}^{-1}>1/a_{\text{IB},-}^{\rm RG}$ there exists a divergence of $G_{-}(\Lambda_{c}^{-})\to\infty$ in the RG at some cut-off $\Lambda_{c}^{-}$, because $G_{-}(0)>0$ and $\partial_{\Lambda}G_{-}^{-1}(\Lambda)<0$. A divergence of $G_{-}$ indicates the existence of new states below the polaron energy for $a_{\text{IB}}^{-1}\geq 1/a_{\text{IB},-}^{\rm RG}$, see Sec. VI.5.2. Because $G_{-}$ describes scattering of a phase fluctuation on the impurity, they contain predominantly phase fluctuations. This is confirmed by the mean-field theory presented in Sec. V.2. Because the impurity atom can only couple to particle number fluctuations due to the underlying density-density interactions, the polaron state is not much affected by the divergence of $G_{-}$.

VIII.2 Number fluctuations

For larger $a_{\text{IB}}^{-1}$ eventually also the interaction strength $G_{+}\to\infty$ diverges at $a_{\text{IB}}=a_{\text{IB},+}^{\rm RG}$. From the exact solution of the RG flow of $G_{+}^{-1}(\Lambda)$ we obtain

[TABLE]

This expression is closely related to the critical scattering length $a_{\text{IB},+}^{\text{MF}}$ where the MF polaron energy diverges, see Eq. (22). The only difference is the appearance of the renormalized mass $\mathcal{M}$ (instead of $M$) in the dispersion relation $\Omega_{\Lambda}$, because of quantum fluctuations included in the RG.

In this case we find that $a_{\text{IB},+}^{\rm RG}$ can become either positive or negative, but obeys the following inequality:

[TABLE]

$a_{\text{IB},+}^{\rm RG}<0$ can become negative only when the mass renormalization is large, $\mathcal{M}\gg M$.

As in the case of $G_{-}$, the divergence of $G_{+}$ indicates the existence of states below the polaron energy for $a_{\text{IB}}^{-1}\geq 1/a_{\text{IB},+}^{\rm RG}$. In this regime $G_{+}(\Lambda)$ diverges at some cut-off $\Lambda_{c}^{+}$ during the RG. Because $G_{+}$ describes the interaction of a particle number fluctuation with the impurity, we expect states bound to the impurity which contain number fluctuations. As discussed in Sec. V.2, in this regime molecular states form together with the phase fluctuations bound to the impurity.

The shift of the particle-number scattering resonance from $a_{\text{IB}}^{-1}=0$ in the two-particle case to $a_{\text{IB},+}^{\rm RG}$ in the many-body case has a simple physical interpretation. Let us first consider an infinitely heavy impurity $M=\infty$, where $a_{\text{IB},+}^{\rm RG}=a_{\text{IB},+}^{\text{MF}}$. In the two-particle case the coupling constants $G_{+}=G_{-}$ diverge in the IR when the energy of the two-particle bound state, which is universally given by $1/2m_{\rm red}a_{\text{IB}}^{2}$, reaches the continuum threshold. This happens for $a_{\text{IB}}\to\infty$. In the many-body case the molecular bound state acquires a density-induced shift to higher energies by an amount $\mu_{\rm BEC}=g_{\rm BB}n_{0}$ due to repulsive boson-boson interactions. Therefore, to first order, it reaches the two-particle scattering continuum when $g_{\rm BB}n_{0}-1/2m_{\rm red}a_{\text{IB}}^{2}=0$. From this simple argument we expect a shift of the resonance to $a_{\text{IB},c}^{-1}=\xi^{-1}>0$. On the other hand, the MF result for $M=\infty$ is $1/a_{\text{IB},+}^{\text{MF}}=\sqrt{2}\xi^{-1}$, which has the same scaling with the coherence length $\xi$.

When the impurity has a finite mass, it becomes increasingly heavy during the RG flow and its effective mass should be replaced by $\mathcal{M}>M$. Therefore the universal binding energy of the dimer $\sim 1/\mathcal{M}a_{\text{IB}}^{2}$ decreases in comparison to the two-particle case. This shifts the energy of the bound state upwards and therefore the bound state reaches the scattering continuum at smaller values of $1/a_{\text{IB}}$, explaining why $1/a_{\text{IB},\pm}^{\rm RG}<1/a_{\text{IB},\pm}^{\text{MF}}$.

VIII.3 Instability of the RG

When $|g_{\text{IB}}^{(0)}\beta(\Lambda)|\to\infty$ in the RG, the MF amplitude $|\alpha_{\bm{k}}|\to\infty$ diverges and the RG breaks down. In this regime quantum fluctuations lead to a collapse of all bosons onto the impurity, and MF theory is no longer justified. We expect that this strong-coupling regime is an artifact of the Bogoliubov approximation because an infinite number of phonons can only be sustained in the polaron cloud if their mutual interactions can be neglected.

We find that the RG breaks down in the extended regime (III in Fig. 1)

[TABLE]

where it runs into the strong coupling phase described above. To understand this, we first note that $\partial_{\Lambda}g_{\text{IB}}^{(0)}\beta(\Lambda)>0$, i.e. the impurity-boson coupling becomes more attractive during the RG. Next recall the initial condition for the RG, $g_{\text{IB}}^{(0)}\beta(\Lambda_{0})=E_{0}^{\text{MF}}(a_{\text{IB}})/n_{0}$. Therefore on the attractive side of the MF resonance, for $1/a_{\text{IB}}<1/a_{\text{IB},+}^{\text{MF}}$ where $E_{0}^{\text{MF}}<0$, we see that $g_{\text{IB}}^{(0)}\beta(\Lambda)$ becomes increasingly more negative during the RG. In fact we can show using Eq. (56) that already for $a_{\text{IB}}=a_{\text{IB},+}^{\rm RG}$ the strong-coupling fixed point is reached where $g_{\text{IB}}^{(0)}\beta(\Lambda=0)=-\infty$ diverges in the IR. In the entire regime (76) $g_{\text{IB}}^{(0)}\beta(\Lambda)$ always diverges during the RG flow. Only for $1/a_{\text{IB}}>1/a_{\text{IB},+}^{\text{MF}}$, where $g_{\text{IB}}^{(0)}\beta(\Lambda_{0})$ starts out positive, a finite value can be reached in the IR limit.

IX Discussion

We have investigated the problem of a mobile impurity interacting with a surrounding BEC in the presence of an inter-species Feshbach resonance. Our theoretical analysis is based on a renormalization group approach to go beyond mean-field calculations and describe quantum fluctuations around the polaron branch non-perturbatively. This gives rise to a logarithmic correction of the polaron energy and we predict new qualitative features appearing in the spectral function close to the Feshbach resonance.

Already on the attractive side, beyond a critical interaction strength $a_{\text{IB},-}^{\rm RG}$, the polaron peak in the spectrum broadens and the quasiparticle weight $Z=0$ vanishes. In this regime there exists a dynamically unstable mode bound to the impurity at energies below that of the polaron. At a second critical interaction strength $a_{\text{IB},+}^{\rm RG}$ the number of phonons in the polaron cloud diverges and beyond this point stable molecular states exist at low energies. We predict an extended strong coupling regime where the number of bosons bound to the impurity diverges as a consequence of quantum fluctuations. In reality this divergency is expected to be prevented by direct phonon-phonon interactions not included in our approach.

At the critical point $a_{\text{IB},+}^{\rm RG}$ we also predict a divergence of the effective polaron mass. It can be extracted experimentally from dipole oscillations of the impurity Catani2012 . Different theoretical approaches predict vastly different polaron masses. In Ref. Rath2013 it was suggested that the polaron hybridizes with the molecule, leading to an effective quasiparticle mass which is only slightly larger than the bare molecular mass. The RG presented in this paper, in contrast, suggests a diverging polaron mass at $a_{\text{IB},+}^{\rm RG}$. Therefore we expect that a measurement of the effective mass will provide valuable physical insights into the Bose polaron problem. More generally we expect that far-from-equilibrium experiments may provide additional insights, beyond the physics that can be probed by a measurement of the spectral function.

We have argued that the attractive polaron can be accurately described by a renormalized Fröhlich Hamiltonian, and the corresponding strong coupling regime can be realized. We suggest to prepare this polaronic state by adiabatically increasing the interaction strength. This scenario can be realized in experiments with multi-component atomic mixtures which have been experimentally realized in laboratories around the world Hadzibabic2002 ; Roati2002 ; Schreck2001 ; Catani2008 ; Spethmann2012 ; Scelle2013 ; Ferrier-Barbut2014 ; Hohmann2015 ; Moses2015 ; Park2015 ; Wacker2015 ; Hu2016PRL ; Jorgensen2016PRL . The mass ratio can be tuned by coupling the impurity atom to a light field Grusdt2016 , which can also be used to realize light polarons in an exciton-polariton condensate Sidler2016 . Such systems should be less vunerable to three-body losses because there are no additional low-lying molecular bound states into which the polaron can decay. This should furthermore allow for an accurate measurement of the three-body losses into the dynamically unstable mode and the molecular state bound to the impurity which we predicted in this paper.

In this work we employed the Bogoliubov approximation and assumed sufficiently weak boson-boson interactions. To describe molecular bound states, or polarons with a large phonon cloud, we expect that residual phonon-phonon interactions need to be considered. In the diffusive Monte Carlo calculations by Ardila and Giorgini Ardila2015 boson-boson interactions have been fully taken into account. The predictions of our RG and the MF method yield considerably lower energies than those predicted in Ardila2015 , and a detailed analysis will be required to understand whether this has to do with a shortcoming of the Bogoliubov approximation. In one dimension a detailed comparison of our RG method with diffusive Monte Carlo calculations has been performed Grusdt2016RG1D and the results are in agreement with recent calculations using flow-equation techniques Volosniev2017 .

Acknowledgements

We would like to thank Christoph Gohle, Lars Wacker, Artur Widera, Atac Imamoglu, Frauke Seeßelberg, Nikolaus Buchheim, Immanuel Bloch, Anatoli Polkovnikov, Michael Fleischhauer, Gregory Astrakharchik, Annabelle Bohrdt, Lode Pollet, Valentin Kasper, Tobias Lausch, Vladimir Stojanovic and Dries Sels for fruitful discussions. We acknowledge support from Harvard-MIT CUA, NSF Grant No. DMR-1308435, AFOSR Quantum Simulation MURI, AFOSR MURI Photonic Quantum Matter. F.G. gratefully acknowledges support from the Gordon and Betty Moore foundation. R.S. acknowledges support by the NSF through a grant for the Institute for Theoretical Atomic, Molecular, and Optical Physics at Harvard University and the Smithsonian Astrophysical Observatory.

Appendix A Comparison to experiments

In Fig. 9 we compare the polaron energies predicted by our RG approach to the peak positions measured recently in radio-frequency absorption spectra Jorgensen2016PRL ; Hu2016PRL .

In the first experiment Jorgensen2016PRL the BEC is very weakly interacting and the model which we use for the RG, based on the Bogoliubov approximation, should be accurate. The agreement of our data with the experiment is reasonable, although we would have expected smaller deviations for weak interactions. We notice however that the agreement with the experimental data is excellent when we shift the resonance position horizontally, but the reason for this is unclear. Our prediction from the RG is in excellent agreement with the theoretical calculations presented in Ref. Jorgensen2016PRL , except at very strong couplings where the spectrum is extremely broad and the comparison in Fig. 9 is not meaningful.

For the second experiment Hu2016PRL we obtain excellent agreement of the repulsive polaron energies. In this case shifting the resonance position does not improve the comparison. For weak attractive interactions the agreement is also excellent. For strong attractive interactions the RG predicts a divergence of the polaron energy which is not seen in the experiment. We attribute this to the fact that, in the experiment, the position of a broad feature with a large weight in the spectrum is recorded. As has been shown in Ref. Shchadilova2016PRL by a calculation of the spectral function, this position does not coincide with the polaron energy in the strong coupling regime. Hence a comparison of our results with the data is meaningless in this regime.

Appendix B Polaronic self-trapping

Investigating the polaron mass allows to study the self-trapping of the impurity predicted in the strong-coupling regime by Landau and Pekar Landau1946 ; Landau1948 . It can be characterized by a dramatic increase of the polaron mass, as observed in Fig. 3, which effectively localizes the impurity Casteels2011 . There has been considerable debate in the literature whether this self-trapping is possible with ultracold atoms, and if it is, whether there exists is a smooth cross-over or a true polaronic phase transition.

So far the discussion in the literature was restricted to the Fröhlich model, where Feynman’s variational path-integral method suggests a sharp self-trapping transition Tempere2009 for sufficiently light impurities Casteels2012 . The applicability of Feynman’s method in the intermediate coupling regime has been questioned, however Feranchuk2005 ; Vlietinck2015 ; Grusdt2015RG . Using more sophisticated methods to describe quantum fluctuations in this regime Tulub1962 ; Shchadilova2016 ; Grusdt2015RG ; Grusdt2016RG ; Kain2016 it has has been shown that self-trapping takes place by a cross-over through an extended intermediate coupling regime. This effect becomes particularly pronounced for light impurities Grusdt2016 .

In Ref. Rath2013 the self-consistent T-matrix approximation has been applied to go beyond the Fröhlich model, but no indications for self-trapping were found. Instead the T-matrix calculations suggest a hybridization of the attractive polaron with a molecule, with an effective mass only slightly exceeding the molecule mass. This is in stark contrast to the predictions of the RG method presented in Fig. 3, where we find a smooth Gerlach1991 self-localization in the spirit of Landau and Pekar on the attractive side of the Feshbach resonance. We expect the scenario described by T-matrix calculations to become relevant for sufficiently large boson-boson interactions however.

Appendix C RG at finite momentum

In this appendix we summarize the RG flow equations for finite polaron momentum $\bm{p}\neq 0$. We keep our discussion general and consider arbitrary dimensionality $d$. In this case the universal Hamiltonian reads

[TABLE]

Here we introduced the new coupling constant $b_{\mu}$, where $\mu,\nu=x,y,...$ are summed over cartesian coordinates using Einstein convention, and we consider a tensor-valued effective impurity mass $\mathcal{M}_{\mu\nu}$. The effective phonon dispersion in the polaron frame reads

[TABLE]

where $Q_{\mu}$ is an additional coupling constant. The universal form of the MF amplitude $\alpha_{\bm{k}}(\Lambda)$, appearing in $\hat{\Gamma}_{\bm{k}}$, is

[TABLE]

where two additional coupling constants $g_{1}(\Lambda)$ and $g_{2}^{\mu}(\Lambda)$ appear.

The initial conditions for the coupling constants are

[TABLE]

The MF polaron state at finite momentum is determined by

[TABLE]

where

[TABLE]

The MF polaron energy is given by

[TABLE]

RG flow equations.– The RG flow equation for the effective impurity mass reads

[TABLE]

For $G_{\pm}$ we obtain

[TABLE]

The flow of $b_{\mu}$, which is only non-vanishing for finite polaron momentum $\bm{p}\neq 0$, is given by

[TABLE]

and for $g_{2}^{\mu}$ we derive

[TABLE]

The flow of $g_{1}$ is determined from the flows of $b_{\mu}$, $g_{2}^{\mu}$ and $\mathcal{M}_{\mu\mu}^{-1}$ (sum over the index $\mu$ is implied),

[TABLE]

Here, and in the following, we assume $\bm{p}=p\bm{e}_{x}$, such that $b_{\mu}=b_{x}\delta_{\mu x}$, $g_{2}^{\mu}=g_{2}^{x}\delta_{\mu x}$, $Q_{\mu}=Q_{x}\delta_{\mu x}$ and $\mathcal{M}_{\mu\nu}^{-1}={\rm diag}(\mathcal{M}_{xx}^{-1},\mathcal{M}_{yy}^{-1},\mathcal{M}_{yy}^{-1},...)$. Moreover we introduced the following integrals (where $J_{5}^{\mu}=I^{x}$ in the spherically symmetric case, see Eq. (55))

[TABLE]

The coupling constant $Q_{\nu}$ has an RG flow given by

[TABLE]

The RG flow of the ground state energy is given by

[TABLE]

Appendix D Quantum theory of molecule formation

In this appendix we introduce a quantum theory of molecular bound states on the repulsive side of the Feshbach resonance, for $1/a_{\text{IB}}>1/a_{\text{IB},+}^{\text{MF}}$. One approach would be to diagonalize the full Hamiltonian (26), which can be done in the limit $M\to\infty$ using Gaussian states Shchadilova2016PRL . Here, instead, we introduce an effective Hamiltonian describing only particle-number and phase fluctuation modes bound to the impurity. To this end we introduce the following two operators,

[TABLE]

which represent two conjugate variables,

[TABLE]

Because particle-number and phase fluctuation modes $\beta_{k}^{\pm}$ are both normalized, $\int d^{3}k~{}(\beta_{k}^{\pm})^{2}=1$, we obtain $\lambda\leq 1$. For sufficiently small boson-boson interactions $\mu_{+}\approx\mu_{-}$ and we thus obtain $\lambda\approx 1$.

In the regime $1/a_{\text{IB},-}^{\text{MF}}<1/a_{\text{IB}}<1/a_{\text{IB},+}^{\text{MF}}$ (II in Fig. 5) only $\mu_{-}<0$ exists. The unstable mode of phase fluctuations can be described by the following effective Hamiltonian,

[TABLE]

This corresponds to a dynamically unstable harmonic oscillator where the potential has negative curvature.

In the regime $1/a_{\text{IB}}>1/a_{\text{IB},+}^{\text{MF}}$ (IV in Fig. 5) both $\mu_{\pm}<0$ exist. To describe both particle-number and phase fluctuations, the following effective Hamiltonian can be used,

[TABLE]

Its variational energy landscape has the same qualitative properties as derived in the proceeding section and summarized in Fig. 5. Using Eq. (102) we can easily diagonalize this Hamiltonian and obtain

[TABLE]

where $[\hat{b},\hat{b}^{\dagger}]=1$ are ladder operators constructed from $\hat{\vartheta}_{b}/\sqrt{\lambda}$ and $\hat{n}_{b}/\sqrt{\lambda}$. The binding energy is given by

[TABLE]

Eq. (105) describes a series of molecular bound states at integer multiples of the binding energy $E_{b}$ below the polaron energy. In Ref. Shchadilova2016PRL we have performed dynamical simulations of the polaron spectrum, where we found a series of equally spaced peaks corresponding to these bound states. Our expression for their energies matches their positions very accurately. Note that this series of molecular bound states has been observed experimentally Schmidt2016 ; Schlagmuller2016 for finite range potentials.

Appendix E Derivation of the RG equations

For the derivation of the RG flow equations presented in Appendix C we choose slightly different representation of the two phonon terms in the universal Hamiltonian. We define

[TABLE]

for $\alpha,\beta=0,1$. Recall that $\hat{d}_{\bm{k}}$ denote the bare boson operators before applying the Bogoliubov transformation. The two-phonon terms can be written as

[TABLE]

and the coupling constants are related by

[TABLE]

E.1 Sketch of the RG procedure

Now we will explain how the RG procedure works, and derive the universal form of the Hamiltonian in Eq. (77). The basic idea is to consider separately slow phonons (S) at momenta $p\leq\Lambda-\delta\Lambda$ and fast phonons (F) from a thin shell $\Lambda-\delta\Lambda<k\leq\Lambda$. Interactions between fast phonons can be ignored, because they are $\mathcal{O}(\delta\Lambda)$, and the fast-phonon Hamiltonian reads

[TABLE]

Slow phonons, on the other hand, are described by the Hamiltonian in Eq. (77), but with all integrals restricted to momenta $|\bm{p}|\leq\Lambda-\delta\Lambda$. We denote the slow-phonon Hamiltonian by $\tilde{\mathcal{H}}_{\text{S}}$. Finally, the coupling between slow and fast phonons is described by

[TABLE]

We made use of the symmetries $\mathcal{M}_{\mu\nu}^{-1}=\mathcal{M}_{\nu\mu}^{-1}$ and $G_{\alpha\beta}=G_{\beta\alpha}$, which can always be enforced by proper ordering of operators.

Next we decouple fast and slow degrees of freedom, by applying a unitary transformation

[TABLE]

The operator $\hat{F}_{\bm{k}}$ depends only on slow variables. To calculate the transformation properties of slow-phonon operators $\hat{O}_{\text{S}}$, we make use of a perturbative expansion in the inverse of the fast-phonon frequency $\Omega_{\bm{k}}^{-1}$, which provides the largest energy scale. It holds

[TABLE]

up to terms of order $\mathcal{O}(\Omega_{\bm{k}}^{-2})$.

The resulting Hamiltonian becomes

[TABLE]

Our basic strategy is to choose $\hat{F}_{\bm{k}}$ such that all terms which linear in $\hat{a}_{\bm{k}}$ (in $\hat{a}^{\dagger}_{\bm{k}}$) vanish, to order $\mathcal{O}(\Omega_{\bm{k}}^{-1},\delta\Lambda)$. In $\hat{U}^{\dagger}_{\Lambda}\tilde{\mathcal{H}}_{\text{MIX}}\hat{U}_{\Lambda}$ the following terms appear: From $\hat{\Gamma}_{\bm{k}}\hat{\Gamma}_{\bm{p}}$ we get $\hat{a}^{\dagger}_{\bm{k}}p_{\mu}\hat{\Gamma}_{\bm{p}}$ because $\hat{\Gamma}_{\bm{k}}$ contains $\alpha_{\bm{k}}\hat{a}^{\dagger}_{\bm{k}}$; Moreover there is a term $\hat{a}^{\dagger}_{\bm{k}}[p_{\mu}\hat{\Gamma}_{\bm{p}},\hat{F}_{\bm{k}}]$ originating from the $\hat{a}^{\dagger}_{\bm{k}}\hat{a}_{\bm{k}}$-term in $\hat{\Gamma}_{\bm{k}}$ and using

[TABLE]

Note the importance of normal-ordering, $\int_{\text{F}}d^{d}kd^{d}k^{\prime}:\hat{a}^{\dagger}_{\bm{k}}\hat{a}_{\bm{k}}\hat{a}^{\dagger}_{\bm{k}^{\prime}}:=\mathcal{O}(\delta\Lambda^{2})$ whereas $\int_{\text{F}}d^{d}kd^{d}k^{\prime}~{}\hat{a}^{\dagger}_{\bm{k}}\hat{a}_{\bm{k}}\hat{a}^{\dagger}_{\bm{k}^{\prime}}=\int_{\text{F}}d^{d}k~{}\hat{a}^{\dagger}_{\bm{k}}+\mathcal{O}(\delta\Lambda^{2})$. Next, from $\hat{D}^{\dagger}~{}_{\bm{k}}\hat{D}_{\bm{p}}$ we derive $\hat{a}^{\dagger}_{\bm{k}}\hat{D}_{\bm{p}}$; The hermitian conjugate $\hat{D}^{\dagger}~{}_{\bm{k}}\hat{D}_{\bm{p}}$ yields $\hat{a}^{\dagger}_{\bm{k}}\hat{D}^{\dagger}~{}_{\bm{p}}$. From the last two terms in Eq. (111) we obtain $\hat{a}^{\dagger}_{\bm{k}}p_{\mu}\hat{\Gamma}_{\bm{p}}$ and $\hat{a}^{\dagger}_{\bm{k}}W_{p}\left(\hat{a}^{\dagger}_{\bm{p}}+\hat{a}_{\bm{p}}\right)$. This final term can be simplified by noting that $W_{p}\left(\hat{a}^{\dagger}_{\bm{p}}+\hat{a}_{\bm{p}}\right)=\hat{D}^{\alpha}_{\bm{p}}+\hat{D}^{\dagger}~{}_{\bm{p}}^{\alpha}$ for arbitrary $\alpha=0,1$.

Therefore, we obtain the following list of terms linear in $\hat{a}^{\dagger}_{\bm{k}}$, ordered in powers of $\Omega_{\bm{k}}^{-1}$,

[TABLE]

where we included the full $\bm{p}$-dependency. We conclude that $\hat{F}_{\bm{k}}$ is of order $\Omega_{\bm{k}}^{-1}$, and given by a certain linear combination of $p_{\mu}\hat{\Gamma}_{\bm{p}}$, $\hat{D}^{\alpha}_{\bm{p}}$ and $\hat{D}^{\dagger}~{}^{\alpha}_{\bm{p}}$. Note that the last two terms in Eq. (117) only contribute higher-order corrections $\mathcal{O}(\Omega_{\bm{k}}^{-2})$ to $\hat{F}_{\bm{k}}$.

Next, we confirm that newly generated terms in the transformed Hamiltonian are all included in the universal Hamiltonian (77). Again, we make a list of such terms. $\Omega_{\bm{k}}\hat{F}^{\dagger}_{\bm{k}}\hat{F}_{\bm{k}}$ contains the following cross-products:

[TABLE]

Note that in the second term, expressions like $\hat{D}^{\dagger}~{}_{\bm{p}}^{\alpha}\hat{D}_{\bm{p}^{\prime}}^{\beta}$ can also be included. For the term in the middle we employed that final expression needs to be hermitian, allowing only the combination $\hat{D}^{\dagger}~{}_{\bm{p}^{\prime}}^{\alpha}+\hat{D}_{\bm{p}^{\prime}}^{\alpha}$ to appear.

Now we check that all terms generated in $\hat{U}^{\dagger}_{\Lambda}\tilde{\mathcal{H}}_{\text{MIX}}\hat{U}_{\Lambda}$ are also included in the list (118). From $\hat{\Gamma}_{\bm{k}}p_{\mu}\hat{\Gamma}_{\bm{p}}$ we obtain terms like $\hat{F}_{\bm{k}}p_{\mu}\hat{\Gamma}_{\bm{p}}$, included in (118). Moreover, from normal-ordering as described above, we get $[p_{\mu}\hat{\Gamma}_{\bm{p}},\hat{F}_{\bm{k}}]$. This yields only terms which are at maximum of order $\mathcal{O}(\hat{a}_{\bm{k}})$, however, and which hence can be eliminated by a MF shift later. In addition there are terms like $\hat{a}^{\dagger}_{\bm{k}}\hat{a}_{\bm{k}}p_{\mu}\hat{\Gamma}_{\bm{p}}$, which renormalize the fast-phonon frequency. We will return to such contributions below. The remaining terms in Eq. (111) lead to similar contributions.

Finally, all terms in Eq. (118) need to be brought into a normal-ordered form. One can check that this only produces terms which are at maximum of order $\mathcal{O}(\hat{a}_{\bm{k}})$. As mentioned previously, these linear terms can be eliminated by a unitary MF-shift, $\hat{a}_{\bm{k}}\to\hat{a}_{\bm{k}}+\delta\alpha_{\bm{k}}$. This in turn leaves the universal Hamiltonian (77) invariant, as will be shown in detail below. Hence we conclude that the RG procedure described above leads to a closed set of RG flow equations for the coupling constants introduced in the beginning. We present the detailed results of our calculations in the following subsection.

E.2 Derivation of the RG

As described above, in every RG step we apply the unitary transformation $\hat{U}_{\Lambda}$ to decouple fast from slow phonons, see Eq. (113). Eliminating terms linear in the fast-phonon operators $\hat{a}_{\bm{k}}$ yields

[TABLE]

E.2.1 Quantum fluctuations

We begin by deriving the properties of quantum fluctuations in the transformed Hamiltonian, which has the form

[TABLE]

The direct renormalization of the slow-phonon Hamiltonian is described by

[TABLE]

The renormalization of the mass is given by

[TABLE]

and the renormalization of $b_{\mu}$ is given by

[TABLE]

For $\delta G_{\alpha\beta}$ we obtain (note the symmetries $G_{00}=G_{11}$ and $G_{01}=G_{10}$)

[TABLE]

The flow of the groundstate energy is determined by

[TABLE]

Furthermore we introduced the following renormalization of the fast-phonon frequency,

[TABLE]

and the MF-terms are characterized by

[TABLE]

E.2.2 Mean-field shift

Next we eliminate the terms linear in the slow phonon operators, $\sim\delta\nu_{\bm{p}}\left(\hat{a}^{\dagger}_{\bm{p}}+\hat{a}_{\bm{p}}\right)$, which give rise to an RG flow of the MF amplitude $\alpha_{\bm{k}}(\Lambda)$. To this end we apply a second unitary transformation, of the form

[TABLE]

Note that $\delta\alpha_{\bm{p}}=\mathcal{O}(\delta\Lambda)$, because the new displacement $\delta\nu_{\bm{p}}$ is also $\mathcal{O}(\delta\Lambda)$.

The saddle-point is determined by eliminating all terms linear in $\hat{a}_{\bm{p}}$ in the (normal-ordered) transformed Hamiltonian $\hat{V}_{\text{MF}}^{\dagger}\hat{U}^{\dagger}_{\Lambda}\tilde{\mathcal{H}}_{\Lambda}\hat{U}_{\Lambda}\hat{V}_{\text{MF}}$. Collecting such terms yields the self-consistency equations for $\delta\alpha_{\bm{p}}$,

[TABLE]

Here we defined the following integrals,

[TABLE]

Next we calculate the transformed Hamiltonian, making use of the fact that the MF shift commutes with normal ordering, $\hat{V}_{\text{MF}}^{\dagger}:\hat{O}:\hat{V}_{\text{MF}}=:\hat{V}_{\text{MF}}^{\dagger}\hat{O}\hat{V}_{\text{MF}}:$. We introduce new $\hat{\Gamma}$-operators by

[TABLE]

Collecting all terms in the transformed Hamiltonian yields

[TABLE]

which, finally, is of the universal form (77) again. As a consequence of the MF shift, the phonon frequency acquires additional renormalization,

[TABLE]

The first term corresponds to mass renormalization, whereas the second term renormalizes $Q_{\mu}$, see Eq. (78).

Before solving the self-consistency equations for $\delta\alpha_{\bm{p}}$, we derive the universal form of $\alpha_{\bm{p}}(\Lambda)$ anticipated in Eq. (80). Comparison of Eqs.(134), (130) directly yields

[TABLE]

Here we can read off the RG flows $\delta g_{1}$ and $\delta g_{2}^{\mu}$, leading to the RG flow equations presented in the follow subsection. While $\delta g_{1}$ depends on $\delta I_{1}$, $\delta I_{2}^{\nu}$, the coupling constant $\delta g_{2}^{\mu}$ is independent of these integrals,

[TABLE]

We can formulate a set of linear equations for $\delta I_{1}$ and $\delta I_{2}^{\mu}$ by plugging Eq. (130) into the definitions of the integrals in Eq. (131). For simplicity, we will now assume, without loss of generality, that the polaron momentum $\bm{p}=p\bm{e}_{x}$ points along $x$-direction. Hence, from symmetry, it follows that $b_{\mu}=\delta_{\mu x}b_{x}$ and thus $\delta I_{2}^{\nu}=\delta_{\nu x}\delta I_{2}^{x}$. We arrive at

[TABLE]

where the integrals were defined in Eqs.(95) - (98). Now it is easy to solve for $\delta I_{1}$ and $\delta I_{2}^{x}$.

Appendix F Dimensional analysis

In this appendix we provide the details of the dimensional analysis which we summarized in Sec. VI.4. After deriving the results for engineering dimensions, we include the anomalous scaling dimensions of the coupling constants.

F.1 Engineering dimensions

We start by deriving the engineering dimensions of all operators appearing in the Hamiltonian, i.e. we neglect the cut-off dependencies of coupling constants generated during the RG flow. This gives a good idea about the importance of the various terms in the Hamiltonian.

High energies, $k\gg 1/\xi$.– In this high energy regime we have $\Omega_{k}\simeq k^{2}/2m_{\rm red}$, such that $\Omega_{\Lambda}\simeq\Lambda^{2}$. Therefore we obtain $\hat{a}_{\Lambda}\simeq\Lambda^{-1-d/2}$. Because $V_{\Lambda}\simeq\Lambda^{0}$ it follows that $\alpha_{\Lambda}\simeq\Lambda^{-2}$. The operator $\hat{\Gamma}_{\bm{k}}$ is a sum of two terms which scale like $\hat{\Gamma}_{\Lambda}\simeq\Lambda^{-d-2}+\Lambda^{-3-d/2}$. Depending on dimensionality, the first or the second term dominates the behavior of $\hat{\Gamma}_{\bm{k}}$ at small cut-offs,

[TABLE]

Because for $k\gg 1/\xi$ it holds $W_{\Lambda},W_{\Lambda}^{-1}\simeq 1$, the scalings of $\delta\hat{n}_{\Lambda}\simeq\hat{\vartheta}_{\Lambda}\simeq\Lambda^{-d/2-1}$ coincide.

Low energies, $k\ll 1/\xi$.– In this low energy regime we have $\Omega_{k}\simeq ck$, such that $\Omega_{\Lambda}\simeq\Lambda^{1}$. Therefore we obtain $\hat{a}_{\Lambda}\simeq\Lambda^{-1/2-d/2}$. Because $V_{\Lambda}\simeq\Lambda^{1/2}$ it follows that $\alpha_{\Lambda}\simeq\Lambda^{-1/2}$. The operator $\hat{\Gamma}_{\bm{k}}$ is a sum of two terms of which $\hat{a}^{\dagger}_{\bm{k}}\hat{a}_{\bm{k}}$ dominates, and it scales as $\hat{\Gamma}_{\Lambda}\simeq\Lambda^{-d-1}$. In the low energy regime we obtain $W_{\Lambda}\simeq\Lambda^{1/2}$. The phase-fluctuations have the same scaling as before, $\hat{\vartheta}_{\Lambda}\simeq\Lambda^{-d/2-1}$, while the particle-number scattering is suppressed, $\delta\hat{n}_{\Lambda}\simeq\Lambda^{-d/2}$.

F.2 Analysis in the absence of divergencies

In the analysis so far we ignored the RG flows of the coupling constants $\mathcal{M}(\Lambda)$ and $G_{\pm}(\Lambda)$. If they change with sufficiently large powers, however, the relative importance of various terms in the Hamiltonian may change. Below we derive the scaling of the coupling constants from the RG flow equations (49) - (51), i.e. their anomalous dimensions.

Here we will consider the regime when neither of the two couplings $G_{\pm}$ diverges during the RG flow. The resulting dimensional analysis is presented in Tab.3 for small energies $k\ll 1/\xi$ and in Tab.4 for large energies $k\gg 1/\xi$.

Scalings of the RG coupling constants.– We begin with the analysis of the effective polaron mass. Its RG flow in the high energy regime ($k\gg 1/\xi$) yields $\partial_{\Lambda}\mathcal{M}\simeq\Lambda^{d-5}$. In $d\leq 3$ dimensions this leads to a monotonic growth of $\mathcal{M}(\Lambda)$, and we obtain $\mathcal{M}\simeq\Lambda^{d-4}$ and thus $\mathcal{M}^{-1}\simeq\Lambda^{4-d}$. At small energies, $k\ll 1/\xi$, the situation is different. There we obtain $\partial_{\Lambda}\mathcal{M}\simeq\Lambda^{d-1}$, meaning that additional renormalization $\Delta\mathcal{M}\simeq\Lambda^{d}$ at small cut-offs becomes irrelevant. We should set $\mathcal{M}(\Lambda)\simeq\Lambda^{0}$ because $\mathcal{M}(\Lambda)\approx\mathcal{M}(1/\xi)$, and accordingly $\mathcal{M}^{-1}(\Lambda)\simeq\Lambda^{0}$.

For $G_{\pm}$ we can perform a similar analysis. For $k\gg 1/\xi$ we obtain that $G_{\pm}^{-1}(\Lambda)=G_{\pm}^{-1}(\Lambda_{0})+\Delta G_{\pm}^{-1}(\Lambda)$, where $\Delta G_{\pm}^{-1}(\Lambda)\simeq\Lambda^{d-2}$. In one dimension $G_{\pm}^{-1}$ grows indefinitely and the term $G_{\pm}^{-1}(\Lambda_{0})$ can be discarded. Here we obtain $G_{\pm}(\Lambda)\simeq\Lambda$. In $d=2,3$ dimensions $\Delta G^{-1}_{\pm}(\Lambda)$ stops to flow (we ignored a slow logarithmic divergence in two dimensions). As long as we are not in the vicinity of a resonance where $G_{\pm}^{-1}=0$ we can thus assume that $G_{\pm}^{-1}(\Lambda)\approx G_{\pm}^{-1}(\Lambda_{0})\simeq\Lambda^{0}$. A similar analysis for $k\ll 1/\xi$ shows that $G_{+}(\Lambda)\simeq\Lambda^{0}$ in any dimension, whereas $G_{-}(\Lambda)$ scales exactly as it does for $k\gg 1/\xi$.

Finally we need to study the RG flow of $\beta(\Lambda)$, which can be done using Eq. (56). In $d=2,3$ dimensions the changes of both $G_{+}^{-1}(\Lambda)$ and $J_{1}(\Lambda)$ become irrelevant at high and low energies (again discarding a peculiar log-divergence in two dimensions). In $d=1$ this is also true for $k\ll 1/\xi$. In these cases we obtain $\beta(\Lambda)\simeq\Lambda^{0}$. The same is true even in one dimension for $k\gg 1/\xi$ as long as we can assume that $\Omega_{\Lambda}\simeq\Lambda^{2}$, which is the case when $m_{\rm B}\leq\mathcal{M}$. The situation where $\mathcal{M}\ll m_{\rm B}$ requires a separate treatment.

F.3 Analysis in the presence of divergencies

Finally we briefly discuss the physics in the vicinity of divergencies of $G_{\pm}$. The behavior of $\mathcal{M}$ is not modified by this. As before we obtain the following scaling for $k\gg 1/\xi$, $\Delta G_{\pm}^{-1}(\Lambda)\simeq\Lambda^{d-2}$. Because at the resonance (at a cut-off $\Lambda_{c}$, say) $G_{\pm}^{-1}(\Lambda_{c})=0$, it follows that $G_{\pm}(\Lambda)\simeq\Lambda^{2-d}$. Accordingly we obtain for $k\ll 1/\xi$ that $G_{-}(\Lambda)\simeq\Lambda^{2-d}$ and $G_{+}(\Lambda)\simeq\Lambda^{-d}$.

In the dimensional analysis of Tabs. 3 and 4 the scalings of two-phonon terms are modified when the respective couplings $G_{\pm}$ diverge. For $d=1,2,3$ the terms become marginal both for $k\ll 1/\xi$ and $k\gg 1/\xi$. I.e. when $G_{+}\to\infty$ diverges,

[TABLE]

and a similar result holds when $G_{-}\to\infty$ diverges. In this case there can be a competition between two-phonon and Fröhlich terms.

Appendix G Two-phonon terms and

the logarithmic UV divergence

In this appendix we show that two-phonon terms do not modify the Fröhlich type UV log-divergence in the order of the RG that we consider. To this end we show that the second term in Eq. (52) for the polaron energy is UV convergent.

We will make use of the asymptotic expression

[TABLE]

Furthermore we note that $m_{\rm red}^{*}=1/\left(m_{\rm B}^{-1}+\mathcal{M}^{-1}\right)=m_{\rm red}\left(1+\mathcal{O}(\Lambda^{-1},\Lambda_{0}^{-1})\right)$, see Eq. (66). Hence we can write to leading order in $\Lambda^{-1}$:

[TABLE]

Now we solve for the RG flows of $G_{\pm}(\Lambda)$. Because the corresponding RG flow equations are separable, we obtain exact expressions in terms of integrals,

[TABLE]

again to leading order in $\Lambda$, $\Lambda_{0}$. Note that $G_{\pm}^{-1}(\Lambda_{0})$ is UV divergent itself,

[TABLE]

Using these expressions, we find that

[TABLE]

Therefore, the second term in the second line of Eq. (141) gives rise to energy contributions scaling like

[TABLE]

which are UV convergent. For the first term in the second line of Eq. (141) we use that

[TABLE]

giving rise to contributions

[TABLE]

which are UV convergent as $\Lambda_{0}\to\infty$.

Appendix H Relation to other logarithmic corrections to the ground state energy

The log-divergence discussed in Eq. (64) is closely related to the logarithmic corrections of the ground state energy in an interacting Bose gas Wu1959 ; Sawada1959 ; Hugenholtz1959 . By considering every particle as an impurity interacting with the surrounding bosons we expect a similar correction $\sim-a^{4}n_{0}^{2}\log\left(\Lambda_{0}/\Lambda_{1}\right)$ to the energy per particle $E_{0}/N$ in this case, where $a$ is the scattering length of the bosons. When the UV cut-off is provided by $\Lambda_{0}\sim 1/a$ as assumed in Ref. Hugenholtz1959 , and considering the case when the mass renormalization is large, we obtain a correction to the energy density $\Delta E_{0}/L^{3}\sim a^{4}n_{0}^{3}\log\left(n_{0}a^{3}\right)$. Indeed this is the correct scaling with $a$ and $n_{0}$ of the leading order correction to the ground state energy due to quantum fluctuations Wu1959 ; Sawada1959 ; Hugenholtz1959 . When the mass renormalization is small the same scaling is obtained by taking into account the dependence of the healing length on the scattering length $a_{\rm BB}=a=a_{\text{IB}}$, see Ref. Christensen2015 . Because $\xi=1/\sqrt{2m_{\rm B}g_{\rm BB}n_{0}}\sim\sqrt{n_{0}a}$ we obtain $\log(\Lambda_{0}\xi)\sim\log n_{0}a^{3}$, again for $\Lambda_{0}\sim 1/a$.

In Ref. Christensen2015 corrections to the Bose polaron energy were calculated exactly up to order $\mathcal{O}(a_{\text{IB}}^{3})$ using diagrammatic perturbation theory. The authors derived a logarithmic UV divergence which scales like $\sim n_{0}a_{\text{IB}}^{3}\xi^{-2}\log(\Lambda_{0}\xi)$. This term is not included in our RG and we expect that it corresponds to a contribution of order $\mathcal{O}(\Omega_{k}^{-2})$ or higher, neglected in our analysis. The authors of Christensen2015 pointed out that this term also gives rise to the correct scaling $\Delta E_{0}/L^{3}\sim n_{0}^{3}a^{4}\log(n_{0}a^{3})$ for the interacting Bose gas (where $a_{\text{IB}}=a_{\rm BB}=a$) Wu1959 ; Sawada1959 ; Hugenholtz1959 .

Our RG approach suggests a unified picture of all logarithmic corrections discussed above. We showed in Grusdt2015Varenna that the polaronic mass renormalization, which is generically given by $\mathcal{M}(k)=M+\mathcal{O}(1/k)$ in the UV limit $k\lesssim\Lambda_{0}$, gives rise to a correction of the impurity energy that diverges logarithmically with the UV cut-off $\Lambda_{0}$. On generic grounds we expect that interaction terms in the Hamiltonian lead to some mass renormalization, connected with the appearance of generic log-divergent terms in the energy. While the log-divergence found by Christensen et al. is related to polaron mass generated by two-phonon terms Christensen2015 , we study mass renormalization due to the Fröhlich terms in the Hamiltonian of the Bose polaron problem. These different mechanisms of mass renormalization explain why the scaling with $a_{\text{IB}}$ in front of the log-divergencies is different for the two situations.