Feynman path-integral treatment of the Bose polaron beyond the Fr\"ohlich model

TL;DR

This paper extends the Feynman path-integral approach to study Bose polarons by including two-phonon interactions beyond the Fr"ohlich model, revealing new effects on ground-state properties and resolving previous discrepancies.

Contribution

It introduces a beyond-Fr"ohlich model for Bose polarons using Feynman's variational approach, capturing two-phonon effects and clarifying previous theoretical inconsistencies.

Findings

No self-trapping for repulsive polarons with two-phonon interactions.

Divergence of energy and effective mass at weaker couplings for attractive polarons.

Discrepancy with renormalization group results is resolved in the extended model.

Abstract

An impurity immersed in a Bose-Einstein condensate is no longer accurately described by the Fr\"ohlich Hamiltonian as the coupling between the impurity and the boson bath gets stronger. We study the dominant effects of the two-phonon terms beyond the Fr\"ohlich model on the ground-state properties of the polaron using Feynman's variational path-integral approach. The previously reported discrepancy in the effective mass between the renormalization group approach and this theory is shown to be absent in the beyond-Fr\"ohlich model on the positive side of the Feshbach resonance. Self-trapping, characterized by a sharp and dramatic increase of the effective mass, is no longer observed for the repulsive polaron once the two-phonon interactions are included. For the attractive polaron we find a divergence of the ground-state energy and effective mass at weaker couplings than previously observed within the Fr\"ohlich model.

Equations

$$\hat{H}=\sum _{\mathbf{j}}^{N_I}\frac{{\hat{\mathbf{p}}}_j^2}{2m_I}+\sum _{\mathbf{k}}{E}_{\mathbf{k}}{\hat{a}}_{\mathbf{k}}^{†}{\hat{a}}_{\mathbf{k}}+\frac{1}{2}\sum _{\mathbf{q}}{V}_{\mathbf{q}}^{BB}\sum _{\mathbf{k},{\mathbf{k}}^{\prime }}{\hat{a}}_{\mathbf{k}+\mathbf{q}}^{†}{\hat{a}}_{{\mathbf{k}}^{\prime }-\mathbf{q}}^{†}{\hat{a}}_{{\mathbf{k}}^{\prime }}{\hat{a}}_{\mathbf{k}}+\sum _{\mathbf{q}}{V}_{\mathbf{q}}^{IB}{\hat{\rho }}_{\mathbf{q}}\sum _{\mathbf{k}}{\hat{a}}_{\mathbf{k}-\mathbf{q}}^{†}{\hat{a}}_{\mathbf{k}}+\frac{1}{2}\sum _{i\ne j}^{N_I}U({\hat{\mathbf{r}}}_i-{\hat{\mathbf{r}}}_j).$$

$$\hat{H}=$$

$$+\frac{g_{ib}}{V}\sum _{\mathbf{k},\mathbf{s}}{\hat{\rho }}_{\mathbf{k}-\mathbf{s}}{W}_{\mathbf{k},\mathbf{s}}^{(1)}{\hat{\alpha }}_{\mathbf{s}}^{†}{\hat{\alpha }}_{\mathbf{k}}+\frac{1}{2}\frac{g_{ib}}{V}\sum _{\mathbf{k},\mathbf{s}}{\hat{\rho }}_{\mathbf{k}-\mathbf{s}}{W}_{\mathbf{k},\mathbf{s}}^{(2)}\left({\hat{\alpha }}_{\mathbf{s}}^{†}{\hat{\alpha }}_{-\mathbf{k}}^{†}+{\hat{\alpha }}_{\mathbf{k}}{\hat{\alpha }}_{-\mathbf{s}}\right)$$

$$+\sum _{\mathbf{j}}^{N_I}\frac{{\hat{\mathbf{p}}}_j^2}{2m_I}+\frac{1}{2}\sum _{i\ne j}^{N_I}U({\mathbf{r}}_{\mathbf{i}}-{\mathbf{r}}_{\mathbf{j}}).$$

$$Q_{\mathbf{k}}=\sqrt{\frac{\hslash }{2M{\omega }_{\mathbf{k}}}}\left({\hat{\alpha }}_{\mathbf{k}}+{\hat{\alpha }}_{-\mathbf{k}}^{†}\right).$$

$$L=$$

$$-\frac{g_{ib}}{V}\frac{M}{2}\sum _{\mathbf{k},\mathbf{s}}{\rho }_{\mathbf{k}-\mathbf{s}}{V}_{\mathbf{k}}{V}_{\mathbf{s}}\frac{\sqrt{{\omega }_{\mathbf{k}}{\omega }_{\mathbf{s}}}}{\hslash }{Q}_{\mathbf{s}}^{\ast }{Q}_{\mathbf{k}}-\frac{g_{ib}}{V}\frac{M\eta }{2}\sum _{\mathbf{k},\mathbf{s}}\frac{V_{\mathbf{k}}^{-1}{V}_{\mathbf{s}}^{-1}}{\hslash \sqrt{{\omega }_{\mathbf{k}}{\omega }_{\mathbf{s}}}}{\rho }_{\mathbf{k}-\mathbf{s}}{\dot{Q}}_{\mathbf{k}}{\dot{Q}}_{\mathbf{s}}^{\ast },$$

$$\eta ={\left(1+\frac{g_{ib}}{V}\sum _{\mathbf{k}}\frac{V_{\mathbf{k}}^{-2}}{\hslash {\omega }_{\mathbf{k}}}\right)}^{-1},$$

$$S=S^F+\frac{g_{ib}}{V}\frac{M}{2}\sum _{\mathbf{k},\mathbf{s}}{\int }_0^{\hslash \beta }d\tau {\rho }_{\mathbf{k}-\mathbf{s}}(\tau )V_{\mathbf{k}}{V}_{\mathbf{s}}\frac{\sqrt{{\omega }_{\mathbf{k}}{\omega }_{\mathbf{s}}}}{\hslash }{Q}_{\mathbf{s}}^{\ast }(\tau )Q_{\mathbf{k}}(\tau ),$$

$$S^F={\int }_0^{\hslash \beta }\left(\frac{m_I{\dot{\mathbf{r}}}^2}{2}+\frac{M}{2}\sum _{\mathbf{k}}{\dot{Q}}_{\mathbf{k}}^{\ast }{\dot{Q}}_{\mathbf{k}}+\sum _{\mathbf{k}}\frac{M{\omega }_{\mathbf{k}}^2}{2}{Q}_{\mathbf{k}}^{\ast }{Q}_{\mathbf{k}}+\frac{\sqrt{N_0}{g}_{ib}}{V}\sum _{\mathbf{k}}{\rho }_{\mathbf{k}}\sqrt{\frac{2M{\omega }_{\mathbf{k}}}{\hslash }}{V}_{\mathbf{k}}{Q}_{\mathbf{k}}\right)d\tau$$

$$e^{-\beta F}=\int \mathcal{D}\mathbf{r}\int \mathcal{D}\{Q_{\mathbf{k}}\}{e}^{-\mathcal{S}\left[\mathbf{r},Q_{\mathbf{k}}\right]/\hslash }=\int \mathcal{D}\mathbf{r}{e}^{-{\mathcal{S}}_{\text{eff}}\left[\mathbf{r}\right]/\hslash },$$

$$F\le F_0+\frac{1}{\hslash \beta }\text{expectationvalue}{\mathcal{S}-{\mathcal{S}}_0}_0,$$

$${\mathcal{S}}_{\text{eff}}^F={\int }_0^{\hslash \beta }\frac{m_I{\dot{\mathbf{r}}}^2}{2}dt-\frac{1}{V}\sum _{\mathbf{k}}\frac{g_{ib}^2{n}_0}{2\hslash }{V}_{\mathbf{k}}^2{\int }_0^{\hslash \beta }d\tau {\int }_0^{\hslash \beta }d\sigma \mathcal{G}\left(\mathbf{k},|\tau -\sigma |\right){\rho }_{\mathbf{k}}(\tau ){\rho }_{\mathbf{k}}^{\ast }(\sigma ),$$

$$\mathcal{G}(\mathbf{k},u)=\frac{\cosh \left[{\omega }_{\mathbf{k}}\left(|u|-\hslash \beta /2\right)\right]}{\sinh \left({\omega }_{\mathbf{k}}\hslash \beta /2\right)}.$$

$${\mathcal{S}}_0={\int }_0^{\hslash \beta }\frac{m_I{\dot{\mathbf{r}}}^2}{2}dt+\frac{M{W}^3}{8}{\int }_0^{\hslash \beta }d\tau {\int }_0^{\hslash \beta }d\sigma \frac{\cosh \left[W\left(|\tau -\sigma |-\hslash \beta /2\right)\right]}{\sinh \left(W\hslash \beta /2\right)}{\left[\mathbf{r}(\tau )-\mathbf{r}(\sigma )\right]}^2.$$

$$e^{-{\mathcal{S}}_{\text{eff}}/\hslash }=\int \mathcal{D}\{Q_{\mathbf{k}}\}\exp \left(-\frac{g_{ib}}{V}\frac{M}{2}\sum _{\mathbf{k},\mathbf{s}}{V}_{\mathbf{k}}{V}_{\mathbf{s}}\frac{\sqrt{{\omega }_{\mathbf{k}}{\omega }_{\mathbf{s}}}}{{\hslash }^2}{\int }_0^{\hslash \beta }d\tau {\rho }_{\mathbf{k}-\mathbf{s}}(\tau )Q_{\mathbf{k}}(\tau )Q_{\mathbf{s}}^{\ast }(\tau )\right)e^{-{\mathcal{S}}^F/\hslash }.$$

$$S^F[J_{\mathbf{k}}]=S^F+\frac{1}{2}\frac{\sqrt{N_0}{g}_{ib}}{V}{\int }_0^{\hslash \beta }\sum _{\mathbf{k}}\sqrt{\frac{2M{\omega }_{\mathbf{k}}}{\hslash }}{V}_{\mathbf{k}}\left[Q_{\mathbf{k}}(\tau )J_{\mathbf{k}}(\tau )+Q_{\mathbf{k}}^{\ast }(\tau )J_{\mathbf{k}}^{\ast }(\tau )\right]d\tau .$$

$$S_{\text{eff}}^F\left[J_{\mathbf{k}}\right]={\int }_0^{\hslash \beta }\frac{m_I{\dot{\mathbf{r}}}^2}{2}dt-\sum _{\mathbf{k}}\frac{g_{ib}^2{n}_0}{2\hslash V}{V}_{\mathbf{k}}^2{\int }_0^{\hslash \beta }d\tau {\int }_0^{\hslash \beta }d\sigma \mathcal{G}\left(\mathbf{k},|\tau -\sigma |\right)\left[{\rho }_{\mathbf{k}}(\tau )+J_{\mathbf{k}}(\tau )\right]\left[{\rho }_{\mathbf{k}}^{\ast }(\sigma )+J_{\mathbf{k}}^{\ast }(\sigma )\right].$$

$$e^{-{\mathcal{S}}_{\text{eff}}/\hslash }={\sum _{n=0}^{\infty }\frac{(-1{)}^n}{n!}{\left[\frac{\hslash }{g_{ib}{n}_0}\sum _{\mathbf{k},\mathbf{s}}{\int }_0^{\hslash \beta }d\tau {\rho }_{\mathbf{k}-\mathbf{s}}(\tau )\frac{\delta }{\delta J_{\mathbf{k}}(\tau )}\frac{\delta }{\delta J_{\mathbf{s}}^{\ast }(\tau )}\right]}^n{e}^{-{\mathcal{S}}_{\text{eff}}^F\left[J_{\mathbf{k}}\right]/\hslash }|}_{J_{\mathbf{k}}=0}.$$

$$e^{-{\mathcal{S}}_{\text{eff}}/\hslash }=\left[1-\frac{g_{ib}}{2\hslash V}\sum _{\mathbf{k},\mathbf{s}}{V}_{\mathbf{s}}^2{\int }_0^{\hslash \beta }d\tau {\rho }_{\mathbf{k}-\mathbf{s}}(\tau )\frac{\delta }{\delta J_{\mathbf{k}}(\tau )}{\int }_0^{\hslash \beta }\mathcal{G}\left(\mathbf{s},\tau -\sigma \right)\left[{\rho }_{\mathbf{s}}(\sigma )+J_{\mathbf{s}}(\sigma )\right]d\sigma \right]{e^{-{\mathcal{S}}_{\text{eff}}^F\left[J_{\mathbf{k}}\right]/\hslash }|}_{J_{\mathbf{k}}=0}.$$

$$e^{-{\mathcal{S}}_{\text{eff}}/\hslash }=\left[1-\frac{g_{ib}\beta }{2V}\sum _{\mathbf{k}}{V}_{\mathbf{k}}^2\mathcal{G}\left(\mathbf{k},0\right)\right]{e^{-{\mathcal{S}}_{\text{eff}}^F\left[J_{\mathbf{k}}\right]/\hslash }|}_{J_{\mathbf{k}}=0}.$$

$$\delta S_{\text{eff}}=\frac{g_{ib}\hslash \beta }{2V}\sum _{\mathbf{k}}{V}_{\mathbf{k}}^2\mathcal{G}(\mathbf{k},0).$$

$$\Delta E_{vac}=\frac{g_{ib}}{2V}\sum _{\mathbf{k}}{V}_{\mathbf{k}}^2.$$

$$e^{-{\mathcal{S}}_{\text{eff}}/\hslash }=\left(1-g_{ib}{O}_1\left[J_{\mathbf{k}}\right]\right){e^{-{\mathcal{S}}_{\text{eff}}^F\left[J_{\mathbf{k}}\right]/\hslash }|}_{J_{\mathbf{k}}=0},$$

$$O_1\left[J_{\mathbf{k}}\right]=\frac{g_{ib}^2{n}_0}{\hslash {\left(2\hslash V\right)}^2}$$

$${\rho }_{{\mathbf{k}}_1}({\tau }_2){\rho }_{{\mathbf{k}}_2}^{\ast }({\tau }_2)\left[{\rho }_{{\mathbf{k}}_2}({\tau }_3)+J_{{\mathbf{k}}_2}({\tau }_3)\right]\mathcal{G}\left({\mathbf{k}}_1,{\tau }_1-{\tau }_2\right)\mathcal{G}\left({\mathbf{k}}_2,{\tau }_2-{\tau }_3\right).$$

$$e^{-{\mathcal{S}}_{\text{eff}}/\hslash }=\left[\sum _{n=0}^{\infty }\frac{(-1{)}^n{g}_{ib}^n{O}_1^n}{n!}\right]e^{-S_{\text{eff}}^F/\hslash }.$$

$$\frac{1}{2!}\frac{\hslash }{n_0}\sum _{\mathbf{k},\mathbf{s}}{\int }_0^{\hslash \beta }d\tau {\rho }_{\mathbf{k}-\mathbf{s}}(\tau )\frac{\delta }{\delta J_{\mathbf{k}}(\tau )}\frac{\delta }{\delta J_{\mathbf{s}}^{\ast }(\tau )}{O}_1\left[J_{\mathbf{k}}\right]{e^{-{\mathcal{S}}_{\text{eff}}^F\left[J_{\mathbf{k}}\right]/\hslash }|}_{J_{\mathbf{k}}=0}.$$

$$g_{ib}^2{O}_2\left[J_{\mathbf{k}}\right]{e^{-{\mathcal{S}}_{\text{eff}}^F\left[J_{\mathbf{k}}\right]/\hslash }|}_{J_{\mathbf{k}}=0},$$

$$O_2\left[J_{\mathbf{k}}\right]=$$

$${\rho }_{{\mathbf{k}}_2}^{\ast }({\tau }_2){\rho }_{{\mathbf{k}}_2}({\tau }_3){\rho }_{{\mathbf{k}}_3}^{\ast }({\tau }_3)\left[{\rho }_{{\mathbf{k}}_3}({\tau }_4)+J_{{\mathbf{k}}_3}({\tau }_4)\right]\mathcal{G}\left({\mathbf{k}}_1,{\tau }_1-{\tau }_2\right)\mathcal{G}\left({\mathbf{k}}_2,{\tau }_2-{\tau }_3\right)\mathcal{G}\left({\mathbf{k}}_3,{\tau }_3-{\tau }_4\right).$$

$$e^{-{\mathcal{S}}_{\text{eff}}/\hslash }=\left[\sum _{n=0}^{\infty }\frac{(-1{)}^n{g}_{ib}^n{O}_1^n}{n!}\right]\left(1+g_{ib}^2{O}_2\right)e^{-S_{\text{eff}}^F/\hslash }.$$

$$g_{ib}^2\frac{{\tilde{O}}_1}{2}{e}^{-{\mathcal{S}}_{\text{eff}}^F/\hslash },$$

Full text

I INTRODUCTION

The generic problem of an impurity interacting with a bath of bosonic excitations has been studied for nearly a century. The concept was first introduced by Landau to describe an electron interacting with an ionic lattice in a solid Landau (1933). Here, the electron induces a polarization cloud in the lattice with which it combines to form a quasiparticle called a polaron. After much progress in the description of the polaron in various coupling regimes Pekar (1946); Landau and Pekar (1948); Fröhlich (1954); Lee et al. (1953) Feynman proposed a variational path-integral description that interpolated previous weak and strong coupling results Feynman (1955). Feynman’s all-coupling description starts from the Fröhlich Hamiltonian Fröhlich (1954) and only incorporates processes of phonon emission or absorption by the electron. In many crystals this is a good approximation and Feynman’s method is considered a very successful description of polarons in solids, although for anharmonic phonons, additional processes need to be considered as well KUSSOW (2009).

Polaronic effects are not limited to polar crystals but have also been observed in ultracold atomic gases where an impurity immersed in the quantum gas becomes dressed by the excitations of the gas. Quantum gases in general have a large experimental tunability, certainly compared to solids, and provide an ideal ground for the study of polaronic physics throughout various interaction regimes Tempere et al. (2009). In experiments impurities can be generated by transferring a small fraction of the gas atoms to a different hyperfine state. The impurity-gas interaction can be tuned from attractive interactions across unitarity towards effective repulsive interactions with Feshbach resonances. The first observations Schirotzek et al. (2009); Koschorreck et al. (2012); Kohstall et al. (2012); Zhang et al. (2012) focused on the Fermi polaron, a single fermionic atom in a sea of opposite-spin fermions. Impurities immersed in a Bose-Einstein condensate (BEC) have also been the subject of a number of experimental studies Heinze et al. (2011); Hohmann et al. (2015); Rentrop et al. (2016) and recently the energy of Bose polarons has been measured in two experiments at lower temperatures Hu et al. (2016); Jørgensen et al. (2016) and near criticality Yan et al. (2019).

The last decade has seen a significant amount of theoretical work towards understanding Bose polarons. At weaker impurity-boson coupling strengths, when the density of excitations in the gas is small, the Hamiltonian describing the problem reduces to the Fröhlich Hamiltonian Tempere et al. (2009); Casteels et al. (2011a, b, 2013a); Blinova et al. (2013); Dasenbrook and Komnik (2013); Vlietinck et al. (2015); Grusdt et al. (2015); Shchadilova et al. (2016a); Grusdt (2016); Ardila and Pohl (2018). However, it has been pointed out using T-matrix calculations Rath and Schmidt (2013) and direct perturbative calculations Christensen et al. (2015) that processes beyond those included in the Fröhlich model cannot be neglected at stronger coupling. In particular two-phonon processes, that we will refer to as “extended Fröhlich interactions” come into play. These are characterized in terms of Feynman diagrams by a vertex joining two impurity lines and two phonon lines. Such processes have also been shown to be crucial for the description of Rydberg polarons Schmidt et al. (2018). This realization has lead to intense theoretical efforts to explicitly incorporate the extended Fröhlich interactions in various analytical methods that were developed for and applied to the Fröhlich model Levinsen et al. (2015); Shchadilova et al. (2016b); Grusdt et al. (2017a); Van Loon et al. (2018). Quantum Monte Carlo (QMC) methods have also been utilized on the full impurity-boson Hamiltonian for the single polaron Ardila and Giorgini (2015); Peña Ardila et al. (2019) and for the bipolaron Camacho-Guardian et al. (2018) as well. In several of these studies, the excitations of the Bose condensate are described within the framework of the Bogoliubov approximation and it remains an open question as to the range of validity of this approximation.

Feynman’s path-integral method, the most successful approach for the solid state polaron, has thus far not been extended to incorporate the two-phonon processes beyond the Fröhlich model. Moreover, while it is known to give extremely accurate results for the polaronic ground-state properties in crystals, its accuracy for the Bose polaron has been questioned as it consistently displays quantitative discrepancies from diagrammatic Monte Carlo (MC) calculations within the Fröhlich-Hamiltonian model at high-momentum cutoff Vlietinck et al. (2015). This concern is largely addressed in Grusdt et al. (2015), where a renormalization group (RG) approach is used to point out that logarithmic divergences that have not been regularized in the MC results in Vlietinck et al. (2015) are not being captured within Feynman’s approach nor within the mean-field (MF) approach, explaining the large discrepancies. Note that beyond weak coupling Feynman’s method still provides a much lower bound for the polaronic contribution to the energy than the mean-field approach. The absence of this logarithmic divergence however suggests that Feynman’s method does not fully capture the same quantum fluctuations as the ones studied in the RG theory at high momentum. At smaller cutoff scales, where the logarithmic regularization is presumably of less importance, Feynman’s approach still yields an energy very close to the MC for the Bose polaron. The other concern pointed out in Grusdt et al. (2015); Grusdt and Demler ; Grusdt and Fleischhauer (2016); Grusdt (2016) is that within the Fröhlich model, Feynman’s method predicts a very sudden and sharp transition to the strong coupling regime, in particular for the effective mass, while this is not expected in the mean-field, RG or other variational models Shchadilova et al. (2016a). In this paper we will show that for the repulsive polaron this artifact is no longer present in Feynman’s approach once extended Fröhlich interactions are added and the transition is smoothed out.

The main aim of this work is to improve the Feynman path-integral method for Bose polarons Tempere et al. (2009) so as to take the extended Fröhlich interactions into account. In section II we start by outlining the general problem in second quantization and find the corresponding Lagrangian required for the path-integral approach. The extended Fröhlich interactions lead to additional quadratic but non-diagonal contributions in the phonon position and velocity variables in the Lagrangian. We perform the path integration over the phonon variables within the Bogoliubov approximation in section III by writing the additional position-dependent extended interaction terms as a full perturbative series with respect to the Fröhlich model. Applying the Jensen-Feynman inequality yields an expression for the variational free energy containing impurity density operator correlations at different times up to an arbitrarily large order. To retrieve an analytic result a random phase-approximation (RPA) of these correlations can be made. In the limit of weak coupling and zero temperature the approximated variational free energy and the effective mass reduce to the extended Fröhlich mean-field results at zero polaron momentum Shchadilova et al. (2016b); Van Loon et al. (2018). This provides a justification of the RPA at weak to intermediate coupling. We find that the addition of the extended Fröhlich interactions allows us to fully regularize the contact potential in a manner similar as in Shchadilova et al. (2016b). In section IV we compare our results for the repulsive polaron to those obtained within the Fröhlich model Tempere et al. (2009) and find significant differences at strong coupling in the ground-state energy, effective mass and polaron radius. The effective mass no longer exhibits a sudden sharp increase, but slowly transitions to the strong coupling regime, accompanied by a finite non-zero polaron radius. This provides an indication against the self-trapping of the repulsive polaron that was present in the Fröhlich model. We also directly compare our results to other recent theoretical models across the resonance. For the energy of the repulsive polaron we find good agreement with the QMC data points presented in Peña Ardila et al. (2019) but quantitative discrepancies with the RG approach Grusdt et al. (2017a), especially in the limit of strong coupling. Nevertheless, the effective mass of the repulsive branch is in good agreement with the predictions of the RG approach. On the attractive branch we find that at weak coupling the variational landscape contains a local minimum that is identified as the polaronic state, but also separated poles. At a critical coupling strength the local minimum vanishes and the polaron energy and effective mass diverge. Similar behavior is observed within RG but not within QMC and our results provide another indication that this early divergence on the attractive branch could be related to the shortcomings of the Bogoliubov approximation Grusdt et al. (2017a), or the importance of correctly capturing Efimov physics in this regime Levinsen et al. (2015); Sun et al. (2017); Yoshida et al. (2018).

II The Hamiltonian and Lagrangian of an impurity in a condensate

The full Hamiltonian describing $N_{I}$ impurities with mass $m_{I}$ immersed in a gas of bosons with mass $m_{b}$ confined to a box with volume $V$ is given by:

[TABLE]

The bosonic creation and annihilation operators $\hat{a}^{\dagger}_{\mathbf{k}},\hat{a}_{\mathbf{k}}$ are written in second quantization, while the impurities are considered in first quantization with impurity position and momentum operators $\mathbf{\hat{r}}_{j},\mathbf{\hat{p}}_{j}$. The first term describes the kinetic energy of the impurities and the second term describes the kinetic energy of the bosonic atoms with $E_{\mathbf{k}}=(\hbar k)^{2}/(2m_{b})$. The third term represents the interactions between the bosonic atoms, where $V^{BB}_{\mathbf{q}}$ is the Fourier representation of the boson-boson potential. The fourth term represents the impurity-boson interaction and similarly contains the impurity-boson potential $V_{\mathbf{q}}^{IB}$ and the impurity density operator $\rho_{\mathbf{q}}=\sum_{j}\exp\left(i\mathbf{q}\cdot\mathbf{\hat{r}}_{j}\right)$. The final term is the impurity-impurity interaction potential $U(\mathbf{\hat{r}}_{i}-\mathbf{\hat{r}}_{j})$. In what follows we will consider the interparticle interactions to be contact interactions $V_{\mathbf{q}}^{BB}=g_{bb}/V$ and $V_{\mathbf{q}}^{IB}=g_{ib}/V$.

Following the Bogoliubov approximation for a weakly interacting BEC we assume the $\mathbf{k}=0$ mode to be macroscopically occupied with $N_{0}$ bosons and rewrite the Hamiltonian in terms of the Bogoliubov-transformed operators $\hat{\alpha}^{\dagger}_{\mathbf{k}},\hat{\alpha}_{\mathbf{k}}$ that create resp. annihilate a Bogoliubov excitation. Keeping up to quadratic order in the Bogoliubov operators one obtains Levinsen et al. (2015); Shchadilova et al. (2016b); Grusdt et al. (2017a); Schmidt et al. (2018):

[TABLE]

The first term, $E_{0}$, is the ground-state energy of the BEC, the ‘vacuum’ energy for the Bogoliubov excitations. The second term is the first-order impurity-boson interaction energy where $N_{B}$ is the total number of bosons. The third term is the kinetic energy of the Bogoliubov excitations with the dispersion relation $\epsilon(\mathbf{k})=\hbar\omega_{\mathbf{k}}=\sqrt{E_{\mathbf{k}}\left(E_{\mathbf{k}}+2n_{0}g_{bb}\right)}$ where $n_{0}=N_{0}/V$. The fourth term is the Fröhlich interaction term, characterized by an effective interaction potential $V_{\mathbf{k}}=\left[E_{\mathbf{k}}/\left({E_{\mathbf{k}}+2g_{bb}n_{0}}\right)\right]^{1/4}$. The next two terms represent the extended Fröhlich interactions where the impurity interacts with two excitations simultaneously, described by effective potentials $W_{\mathbf{k,k^{\prime}}}^{(1)}=\frac{1}{2}\left(V_{\mathbf{k}}V_{\mathbf{k^{\prime}}}+V_{\mathbf{k}}^{-1}V_{\mathbf{k^{\prime}}}^{-1}\right)$ and $W_{\mathbf{k,k^{\prime}}}^{(2)}=\frac{1}{2}\left(V_{\mathbf{k}}V_{\mathbf{k^{\prime}}}-V_{\mathbf{k}}^{-1}V_{\mathbf{k^{\prime}}}^{-1}\right)$. In the rest of this work we will not consider depletion effects and approximate $N_{B}$ by the number of condensed atoms $N_{0}$.

To apply the path-integral formalism, the Lagrangian corresponding to (2) is needed. For this, the Bogoliubov creation and annihilation operators are combined into effective phonon coordinates $Q_{\mathbf{k}}$ following the standard prescription:

[TABLE]

Here, $M$ is an arbitrary phonon mass which will not appear in the effective action of the polaron, not to be confused with the variational mass of the model system. The Lagrangian, written in terms of these phonon degrees of freedom is derived in appendix (A) for a single impurity and is given by:

[TABLE]

where

[TABLE]

and $\mathbf{r}$ is the impurity coordinate. It is clear that if we do not perform any regularization procedures of the contact interaction, $\eta$ vanishes at infinite cutoff. The first four terms of (4) correspond to the Lagrangian of the Fröhlich model and yield the Fröhlich action $\mathcal{S}^{F}$. The additional terms in (4) take into account the extended interactions beyond the Fröhlich model and consist of a part depending on the phonon coordinates $Q_{\mathbf{k}}$ multiplied by $V_{\mathbf{k}}$, and a part depending on the phonon velocities $\dot{Q}_{\mathbf{k}}$ along with $V_{\mathbf{k}}^{-1}$. Within a mean-field approach where the excitation operators acquire a polaronic shift $\hat{\alpha}_{\mathbf{k}}\rightarrow\hat{\alpha}_{\mathbf{k}}-f_{\mathbf{k}}$, these velocity-dependent terms can be shown to arise due to a non-zero imaginary contribution from $f_{\mathbf{k}}$ and vanish for the saddle-point solution Shchadilova et al. (2016b); Van Loon et al. (2018), while the terms containing $V_{\mathbf{k}}$ arise due to the real part of $f_{\mathbf{k}}$ and have a non-negligible contribution to the ground-state energy resulting in a resonance shift. In the RG approach Grusdt et al. (2017a) it is pointed out that the RG coupling constant corresponding to the $V_{\mathbf{k}}^{-1}$ terms has a small effect on the polaron wavefunction, but is expected to be important when considering other qualitative properties such as the lifetime of the polaron due to the appearance of bound states at lower energies. These considerations are, however, beyond the scope of this work and we will only consider the position-dependent terms of the extended interactions. In this way we capture the same effects as the mean-field treatment but treat them beyond the mean-field level. Hence, in the remainder of this work, we consider the Euclidean polaron action functional,

[TABLE]

where

[TABLE]

is the action of the Frohlich model.

III Feynman’s variational path-integral approach

The free energy $F$ of the polaron can be expressed as a path integral over the impurity and phonon degrees of freedom, weighted by the exponent of the action functional $\mathcal{S}\left[\mathbf{r},Q_{\mathbf{k}}\right]$ corresponding to Lagrangian (4):

[TABLE]

where $\beta=(k_{B}T)^{-1}$ is the inverse temperature and $\mathcal{S}_{\textrm{eff}}\left[\mathbf{r}\right]$ is the effective imaginary-time action where the phonon degrees of freedom have been integrated out. For a variational model system with action $\mathcal{S}_{0}$ and free energy $F_{0}$ the Jensen-Feynman inequality provides a variational upper bound for the free energy of the polaron Feynman (1955) (see Feynman (1998); Kleinert (2009) for details):

[TABLE]

where $\expectationvalue{...}_{0}$ is the expectation value with respect to the variational model system.

III.1 Outline within the Fröhlich model

Within the Fröhlich model (so, excluding extended interactions) the effective action for a Bose polaron is given by Tempere et al. (2009):

[TABLE]

with $n_{0}$ the BEC density and $\rho_{\mathbf{k}}(\tau)=\exp\left[i\mathbf{k}\cdot\mathbf{r}(\tau)\right]$. As a consequence of integrating out the phonon variables the effective action (10) now contains a retarded interaction mediated by the Green’s function of the Bogoliubov excitations:

[TABLE]

We will consider the same variational model system $\mathcal{S}_{0}$ as in Feynman (1955); Tempere et al. (2009). This system physically corresponds to a particle with the same mass as the impurity $m_{I}$ coupled to a second mass $M$ by a spring constant $MW^{2}$, $M$ and $W$ being variational parameters. Hence, the variational model action is given by:

[TABLE]

The free energy of the model system $F_{0}$ can be straightforwardly calculated and the variational upper bound for the polaron energy (9) can be minimized as a function of $M$ and $W$. For the Fröhlich model this has been done by Feynman for a polaron in a crystal Feynman (1955) and more recently applied to a Bose polaron as well Tempere et al. (2009). In the following subsection we will go beyond the Fröhlich model for the Bose polaron by including the extended Fröhlich interactions in the effective action of the polaron system $\mathcal{S}_{\textrm{eff}}$.

III.2 Perturbative expansion for the beyond Fröhlich terms

The effective action corresponding to (6) is obtained by integrating out the phonon degrees of freedom. First, $\exp({-\mathcal{S}}/\hbar)$ is factorized in the Fröhlich contribution $\exp({-\mathcal{S}^{F}}/\hbar)$ and the beyond-Fröhlich part:

[TABLE]

The idea is to take into account the exponential of the beyond-Fröhlich terms in (13) perturbatively through a series expansion of the exponential and a subsequent integration over the phonon degrees of freedom. The terms in the resulting perturbation series can be obtained more straightforwardly with the generating functional formalism. The generating functional is obtained by adding source terms (and a prefactor that will simplify the algebra) to the Fröhlich action:

[TABLE]

The source terms resemble the Fröhlich impurity-phonon interaction term, and can be added to it. Hence, including the source terms in the Fröhlich action amounts to shifting $\rho_{\mathbf{k}}$ to $\rho_{\mathbf{k}}+J_{\mathbf{k}}$. The generating functional is then obtained by integrating out the phonon degrees of freedom. The resulting effective action of the Fröhlich model including source terms then becomes:

[TABLE]

In the series expansion of the exponential in (13) the phonon position variables $Q_{\mathbf{k}}$ can be replaced by functional derivatives with respect to $J_{\mathbf{k}}$, which can be brought out of the functional integral over $Q_{\mathbf{k}}$. After performing the path integral over the phonon variables one is left with the following expression:

[TABLE]

We will now provide an overview of the structure of the various terms appearing in the generating functional series (16) and argue which terms can be neglected. They can be classified in three categories.

III.2.1 Vacuum energy terms

It is illustrative to consider the $n=1$ order in the expansion. After the first $\frac{\delta}{\delta J^{*}_{\mathbf{s}}(\tau)}$ in (16) is applied to the exponential, one obtains

[TABLE]

Now there is a choice whether to apply the second operator $\frac{\delta}{\delta J_{\mathbf{k}}(\tau)}$ to the exponential again or to the $J_{\mathbf{s}}(t)$ in front. The former option leads to terms which will be discussed in the next subsection. The latter option fully eliminates the impurity variable and results in:

[TABLE]

In every order $n$ of the expansion there will be $0\leq j\leq n$ operator pairs in which each individual pair is applied in the same way as in this example to eliminate the impurity variable and merely yield the term in (18) to the power $j$ multiplied with a combinatorial factor. This allows to separate these terms and perform their complete series summation, given that they are not to be counted from this point on. The summation results in the following contribution to the effective action:

[TABLE]

At zero temperature the corresponding energy shift is given by:

[TABLE]

This contribution is precisely of the same type as the divergent terms arising from non-commuting variables in the derivation of the Lagrangian in appendix (A). First-order corrections in $g_{ib}$ to the ground-state energy are not observed in a rigorous perturbative calculation Christensen et al. (2015) and we do not expect these terms to be of physical significance. Furthermore note that (20) is UV divergent and can not be regularized by taking the cutoff dependence of $g_{ib}$ into account. Therefore we will discard contribution (19) in the rest of our calculations.

III.2.2 Scattering terms

In the previous example of expansion order $n=1$ we could have also applied the second functional derivative to the exponential again in (17) to find:

[TABLE]

where $O_{1}\left[J_{\mathbf{k}}\right]$ is given by

[TABLE]

This term can be interpreted in relation to a process where an impurity creates an excitation out of the BEC at time $\tau_{1}$, scatters with this excitation at time $\tau_{2}$, and finally returns it to the BEC at time $\tau_{3}$. Hence we will refer to this term as the first-order scattering term. Every higher order term in the expansion will contain precisely one combination where every pair of functional derivatives is applied only to the exponential and contributes a power of $O_{1}$. A short calculation shows that these terms form the exponential power series:

[TABLE]

Here, we use $O_{1}$ as a notation for $O_{1}\left[J_{\mathbf{k}}=0\right]$, i.e. where the source terms are set to zero. We can proceed to derive the second-order scattering term $O_{2}$. The $n=2$ term in the expansion of (16) can be written as:

[TABLE]

Applying one of the two functional derivatives in (24) to the exponential and the other to $O_{1}\left[J_{\mathbf{k}}\right]$ and vice versa will result in two terms that are combined in:

[TABLE]

where

[TABLE]

For every term of order $n>2$ in the expansion there will be a combination of functional derivatives that will result in a contribution $\sim(n)(n-1)O_{1}^{n-2}O_{2}$. Performing the explicit calculation shows that these terms can be combined with the series of first-order scattering terms in (23) in the following way:

[TABLE]

We will return to the factorization pattern appearing in (27) later. First we have to address the terms that are not included in this reasoning.

III.2.3 Excitation bath terms

In the $n=2$ order of the expansion we could have also chosen to apply both functional derivatives to $O_{1}\left[J_{\mathbf{k}}\right]$ in expression (24), which would result in:

[TABLE]

where

[TABLE]

This term can be related to a process where an impurity exchanges momentum with the excitation bath without coupling to the BEC, i.e. without creating an excitation from the condensate or scattering it back to the condensate. Whereas the scattering terms discussed in the previous subsection, such as (22), are proportional to $n_{0}$ (the number of atoms in the condensate), the enhancement factor $n_{0}$ is absent in the excitation bath terms such as (29).

The higher-power contributions of this term will present themselves as $\tilde{O}_{1}^{n/2}$ in all the even $n>2$ expansion terms and together with (28) they form an exponential power series as well:

[TABLE]

Just like for the scattering terms, this reasoning can be extended to higher order terms $\tilde{O}_{n}$ or to the combination of these terms with the scattering terms such as $\tilde{O}_{1}O_{1}^{n-2}$. In general, terms uncoupled from the condensate, arise when both functional derivatives $\frac{\delta}{\delta J_{\mathbf{k}}(\tau)}\frac{\delta}{\delta J^{*}_{\mathbf{s}}(\tau)}$ in a pair in (16) are applied to the two source terms contained in a scattering term $O_{n}\left[J_{\mathbf{k}}\right]$. As mentioned above, the main difference between the scattering terms and the excitation bath terms is that the former contain a leading order $n_{0}$, while the latter do not and are therefore suppressed by a relative factor of $a_{bb}/\xi$, with $a_{bb}$ the boson-boson scattering length and $\xi$ the coherence length of the BEC. Within the range of validity of the Bogoliubov approximation, i.e. $(n_{0}a_{bb}^{3})\ll 1$, they are negligible. Note that a similar argument has been made in a perturbative calculation Christensen et al. (2015) to ignore diagrams where the impurity couples to bosons outside of the BEC. We will therefore not include these excitation bath terms in the rest of our calculations.

III.2.4 Result

In the discussion of the scattering terms above we have found that the power series in the first-order terms $O_{1}$ and all the product terms $O_{2}O_{1}^{n-2}$ compactly factorize in expression (27). An explicit calculation shows that this factorization pattern extends to higher-order scattering terms $O_{n}$:

[TABLE]

Here $O_{n}$ represents an $n$-th order scattering process where an impurity creates an excitation out of the BEC and scatters with it $n$ times before scattering it into the condensate again.

[TABLE]

The factorization pattern appearing in (31) suggests that the effective action can be written as:

[TABLE]

It might be illustrative to point out that even within the conventional Fröhlich model a similar structure can be observed. Performing a perturbative expansion of the Fröhlich contribution in (4) with respect to the free impurity yields:

[TABLE]

where $S_{\textrm{free}}$ is the action functional of a non-interacting impurity and $O_{0}$ is a “zeroth-order” scattering term characterizing an impurity that creates an excitation and absorbs it a time later, without any interaction in between. This series can be recombined in the effective action to obtain precisely the action in expression (10).

Applying the Jensen-Feynman inequality (9) to the effective action (33) will yield a variational free energy that contains impurity density correlation functions at different times, of an arbitrarily large order corresponding to the number of scattering events in the scattering terms. This can explicitly be seen from the product in (32).

III.3 Random phase approximation

To proceed analytically, an approximation of the impurity density correlations can be made. Relative to the model system $\mathcal{S}_{0}$, the impurity density correlation between the creation of an excitation at time $\tau_{j}$ and its absorption at time $\tau_{j+1}$ depends only on the absolute value of the time step Tempere et al. (2009):

[TABLE]

$\mathcal{F}_{\mathbf{k}}(u)$ is the memory function of the impurity:

[TABLE]

We consider a random phase approximation (RPA) where the dominant contribution to the correlation of a number of subsequent scattering events is given by the correlations within one scattering event (35):

[TABLE]

Using $\mathcal{G}_{\mathbf{k}_{j}}\left(\hbar\beta-u\right)=\mathcal{G}_{\mathbf{k}_{j}}\left(u\right)$ and $\mathcal{F}_{\mathbf{k}_{j}}\left(\hbar\beta-u\right)=\mathcal{F}_{\mathbf{k}_{j}}\left(u\right)$ it is not difficult to show that within this approximation the additional contributions to the effective action constitute a power series:

[TABLE]

After substituting the expectation value of the Fröhlich effective action (10) with respect to the variational model system and adding the first-order energy contribution $g_{ib}n_{0}$ the expectation value of the full effective action becomes:

[TABLE]

To regularize the series we substitute the full Lippmann-Schwinger equation:

[TABLE]

where $\mu=\left(m_{b}^{-1}+m_{I}^{-1}\right)^{-1}$ is the reduced impurity-boson mass and $a_{ib}$ is the impurity-boson scattering length. The free energy $F_{0}$ of the model system and expectation value of the action of the model system $\expectationvalue{\mathcal{S}_{0}}_{0}$ can be computed. Substitution into the Jensen-Feynman inequality (9) yields the following variational free energy:

[TABLE]

The variational parameters are $\Omega$ and $M$ and the relation between $\Omega$ and the original oscillator frequency in the model system is given by $\Omega=W\sqrt{1+M/m_{I}}$, see Feynman (1955, 1998); Tempere et al. (2009) for a detailed description. The free energy (41) is written in a suggestive form to make the analogy with the resonance shift observed in Shchadilova et al. (2016b); Grusdt et al. (2017a). The resonance shift is UV convergent and in our case depends on both the temperature and the variational parameters:

[TABLE]

The free energy (41) contains the first-order energy contribution $2\pi\hbar^{2}a_{ib}n_{0}/\mu$ as well, which is not included in the expression given in Tempere et al. (2009). As a consistency check we consider the limit of weak coupling with a simplified model system where the phonon mass of the Feynman model approaches zero, $M\rightarrow 0$, while the spring constant $MW^{2}$ remains fixed. At zero temperature ($\beta\rightarrow\infty$), the energy in the weak coupling limit $E_{\text{weak}}$ is independent of the variational parameters and given by:

[TABLE]

In this limit the $u$-integral in $a^{-1}_{0,\text{weak}}$ can be analytically performed:

[TABLE]

Expression (43) with the resonance shift (44) is precisely the mean-field result including extended interactions at zero polaron momentum $\mathbf{P}=0$ Shchadilova et al. (2016b); Grusdt et al. (2015); Van Loon et al. (2018).

Feynman’s path-integral formalism allows us to calculate an effective mass for the polaron $m_{\textrm{pol}}$ and a root mean square (RMS) estimate of the polaron size $\sqrt{\expectationvalue{\mathbf{r}^{2}}}$. The expression for the polaron radius depends only on the model system and remains the same as in Tempere et al. (2009):

[TABLE]

An expression for the effective mass can be derived by introducing a boost to the memory function of the system $\expectationvalue{\exp\left(i\mathbf{k}\left[\mathbf{r}(\tau)-\mathbf{r}(\sigma)\right]\right)}_{0}\rightarrow\expectationvalue{\exp\left(i\mathbf{k}\left[\mathbf{r}(\tau)-\mathbf{r}(\sigma)\right]\right)}_{0}\times\exp\left(i\mathbf{k}\cdot\mathbf{v}\left(\tau-\sigma\right)\right)$. This method is used in Feynman’s seminal work on polarons within the Fröhlich model Feynman (1955, 1998). However, to incorporate the effects on the effective mass of the extended Fröhlich contributions we only apply this boost to the effective action after the RPA contributions have been separated in (39). After deriving the energy as a function of $\mathbf{v}$ and expanding it up to $\mathbf{v}^{2}$, the factor in front of $\mathbf{v}^{2}/2$ can be identified as the polaron effective mass:

[TABLE]

where $\Gamma$ is given by:

[TABLE]

As far as we are aware, Feynman’s prescription is only valid in the low temperature limit, so the limit $\beta\rightarrow\infty$ in expression (46) must be taken. The effect of temperature is then estimated through the implicit temperature dependence of the variational parameters $M$ and $\Omega$, as has been done in Casteels et al. (2012). Note that in the limit of $M\rightarrow 0$ in the model system, our result for the effective mass reduces to the extended Fröhlich interactions mean-field result Van Loon et al. (2018) which suggests that corrections beyond the mean-field level are captured within this method.

IV Results

IV.1 Comparison with the Fröhlich model for the repulsive polaron

First, we make a direct comparison between the results obtained with Feynman’s variational description within the Fröhlich model in Tempere et al. (2009) and the results including extended interactions derived in the previous section. Because the polaronic contribution to the free energy within the Fröhlich model is the same on both sides of the resonance we will only consider the repulsive polaron in this section. It is important to note that on this side of the resonance various shallow bound states do exist at lower energies Rath and Schmidt (2013); Shchadilova et al. (2016b) and we are only retrieving the energy of the repulsive branch in our approach. The rich physics of Efimov bound states for an impurity in a BEC Levinsen et al. (2015); Sun et al. (2017); Yoshida et al. (2018) is not expected to be captured in this approach. Within the Fröhlich model, the results at a given temperature can be expressed as a function of a single dimensionless polaronic coupling constant $\alpha=a_{ib}^{2}/(\xi a_{bb})$. However, for (41) this is no longer the case, as the results depend also explicitly on $a_{bb}$. Nevertheless, for a fixed $a_{bb}$ we can still plot our results as a function of $\alpha$ at the repulsive side of the resonance for the purpose of the comparison. A physical cutoff corresponding to the range of the interatomic interaction is used, given by $\Lambda_{c}\approx 200\xi^{-1}$ for the current system.

Figure 1 shows the results for the polaronic contribution to the free energy $F_{p}=F-2\pi\hbar^{2}a_{ib}n_{0}/\mu$, in polaronic units ($\xi=m_{I}=\hbar=1$) at various temperatures. At weak coupling both results coincide but they start to significantly differ around $\alpha\approx 3.5$ where the Fröhlich model predicts a very steep decrease in energy, indicative of self-trapping. The extended interactions appear to moderate this into a much slower linear decrease of the free energy. The decrease of the polaronic contribution is even slower than the increase of the first-order contribution $2\pi\hbar^{2}a_{ib}n_{0}/\mu$, and the full polaron energy for the extended interactions model never becomes negative.

Figure 2 presents a comparison of the polaron RMS radius (45) between the two models. The first noticeable difference is that the sharp kink within the Fröhlich model, previously identified with the transition into the strong coupling regime around $\alpha=3.5$, is replaced by a smoother non-monotonic transition due to the extended interactions. Most significant is the difference at extremely strong coupling however. The inclusion of the extended interactions disproves previous predictions of the asymptotically shrinking Bose polaron within this method, which was suggestive of self-trapping as well, and shows that the polaron radius approaches a finite non-zero value around $\approx 0.35\xi$. Comparable conclusions follow for the effective mass of the polaron (46), shown in fig. 3. The effective mass no longer exhibits a sudden and steep transition into the strong coupling regime. A period of faster increase of the effective mass is still observed around $\alpha=5\text{-}10$, but flattens out towards a value of roughly $\approx 20m_{I}$ at even stronger coupling. Furthermore we can see that in the case of a light impurity such as considered here, the effective mass is more sensitive to temperature differences than the energy and radius. It has been pointed out that measurements of the effective mass of the polaron are expected to be particularly useful to discern between various theoretical models Grusdt et al. (2015); Grusdt and Demler . Based on our results we expect this to be even more the case when the temperature dependence is measured as well.

The converging effective mass and polaron radius together with the positive free energy suggest that self-trapping does not take place for the repulsive polaron when the extended interactions are included. This is qualitatively in agreement with the findings of the RG approach Grusdt et al. (2017a), where no self-trapping is observed for the repulsive polaron.

IV.2 Comparison with other theoretical results

In this subsection we provide a comparison with other recent results in the literature, specifically with the mean-field approach Shchadilova et al. (2016b), the RG approach Grusdt et al. (2017a) and Quantum Monte Carlo calculations Peña Ardila et al. (2019). As mentioned above, the mean-field expressions for the energy and effective mass Shchadilova et al. (2016b); Van Loon et al. (2018) can also be obtained from the weak coupling limit of Feynman’s model.

Before proceeding to the discussion, one aspect of the Feynman model has to be addressed. As can be seen from expression (41), the variational landscape can contain poles where the free energy diverges to negative (or positive) infinity, accompanied by a divergence of the effective mass (46). On the negative side of the resonance, even at weaker coupling, these poles are present. However, below a critical coupling strength there exists a separated local minimum that corresponds to the polaronic state. To plot the polaron energy of the attractive branch we follow this local minimum starting from weak interactions up to the point where it merges with one of the aforementioned poles, at which both the energy and effective mass diverge.

To better understand the physical significance of these divergences, it is illustrative to observe that the same type of pole is present in the extended mean-field treatment Shchadilova et al. (2016b); Van Loon et al. (2018), where it is independent of any additional variational parameters. In the MF model this divergence can be shown to be accompanied by a rapid depletion of the BEC, which is no longer accurately described within Bogoliubov theory. We therefore believe that the poles observed within our treatment can be interpreted as a runaway pathway related to the shortcomings of the Bogoliubov approximation. A detailed discussion of similar divergences, observed in RG theory, is presented in Grusdt et al. (2017a). Note that in Yoshida et al. (2018) no divergences are observed for the polaron at unitarity, which indicates that a correct treatment of Efimov physics is of importance here as well.

For the repulsive polaron no runaway pathways exist at weak coupling and we simply follow the global minimum of the variational landscape. Only at extremely strong coupling, separated divergences start to appear and the polaronic state becomes a local minimum. This local minimum continues to exist across the resonance towards negative scattering lengths, and it is not clear if we can interpret it as the repulsive polaron state from this point on. For the purpose of the comparison with RG in this subsection, we will restrict our study of the repulsive branch to couplings below the critical coupling presented in Grusdt et al. (2017a). At this point the aforementioned transition into a local minimum has not yet taken place.

Figures 4 and 5 show a comparison of the polaron energies obtained with various methods, across the resonance for impurity-condensate parameters used in the experiments Jørgensen et al. (2016) and Hu et al. (2016), respectively. To provide an accurate comparison with RG, the same respective finite values of the momentum cutoff were used as in Grusdt et al. (2017a). As shown on the inset of fig. 4, further convergence of a few percent is expected at infinite cutoff in this case. The cutoff used for fig. 5 is larger and the results are much closer to convergence. The cutoff dependence of the Feynman approach within the Fröhlich BEC model has been discussed in Tempere et al. (2009) and Vlietinck et al. (2015). For both calculations finite temperatures, at which the energy has converged beyond any noticeable difference in the figures, were used to represent zero temperature. This convergence is shown on the inset of fig. 5.

For the repulsive polaron we observe a relatively good quantitative agreement with QMC data in fig. 4 and an excellent agreement in fig. 5. Note, however, that the QMC calculation does not rely on the Bogoliubov approximation. Our results predict no divergence of the repulsive branch energy in contrast to the mean-field treatment, but towards stronger coupling a quantitative discrepancy with the RG approach appears. However, as shown in fig. 6, a much better agreement exists for the effective mass of the repulsive polaron between the two methods. One possible explanation for the discrepancy in energy is the previously discussed logarithmic divergence captured in the RG theory. The QMC study Peña Ardila et al. (2019) does not elaborate on the cutoff dependence so the status of the logarithmic divergences in this method is unclear. On the negative side of the resonance we see qualitative agreement with RG where the polaron energy and effective mass diverge at a weaker interaction strength than predicted by the mean-field description or QMC.

The theoretical results can also be compared to the experimental data points from Hu et al. (2016); Jørgensen et al. (2016), which we have not explicitly added to the figures for the purpose of clarity. At weak coupling all theoretical approaches are in excellent agreement with experiments. For the repulsive branch as the coupling gets stronger, the data points of Jørgensen et al. Jørgensen et al. (2016) lie at higher energies than QMC, even after the non-homogeneity of the three-body decay processes is taken into account in the spectroscopic signal Peña Ardila et al. (2019). As both RG and our results lie below QMC it follows that the agreement with this experiment is not close in this regime. Several reasons for this discrepancy are suggested in Peña Ardila et al. (2019). In contrast we find excellent agreement with the experiment of Hu et al. Hu et al. (2016) on the repulsive branch, which is also in much closer agreement with QMC than the Jørgensen et al. experiment. The early divergence for the attractive polaron that is found in our results and the RG approach is observed in neither experiment, and as mentioned above, its understanding requires a further study of the validity of the Bogoliubov approximation and the effects of Efimov physics in that regime.

V Conclusion

We have studied the ground-state properties of the Bose polaron beyond the Fröhlich paradigm using Feynman’s variational path-integral formalism. For this purpose we derived the Lagrangian of an impurity immersed in a condensate within the Bogoliubov approximation. The extended Fröhlich interactions take the form of quadratic position- and velocity-dependent terms in the phonon variables. By expanding the position-dependent terms as a full perturbative series the path integral over the phonon variables can be performed to obtain an effective action. This is done within the Bogoliubov approximation, neglecting perturbative terms that contain no coupling to the condensate. We do not expect the velocity-dependent terms to contribute significantly to the ground-state properties based on other theoretical studies in the literature. The Jensen-Feynman inequality provides a variational expression for the upper bound on the free energy. Due to the extended interactions it contains a series of impurity density correlations that, as far as we know, does not reduce to an analytic expression. To proceed analytically a random phase approximation is made that decomposes the higher order impurity-excitation scattering correlations as a product of subsequent scattering correlations. The RPA yields simple variational expressions for the polaron energy and effective mass that reduce to the extended Fröhlich mean-field results in the weak coupling limit. For the repulsive polaron we compared the predictions with those of the Fröhlich model and found that the sharp transition to the strong coupling regime, which was interpreted as a possible shortcoming of the path-integral approach for the Fröhlich model, is now replaced by a smooth crossover suggestive of the absence of self-trapping. For the attractive polaron we observed an abrupt divergence of the energy and effective mass at a certain critical coupling strength. This is related to the local polaronic minimum getting absorbed by a runaway pathway in the variational landscape, and is interpreted as a breakdown of the Bogoliubov approximation within our approach.

Various future perspectives for this method exist. While Feynman’s method in theory captures the full effect of the excitations at the level of the effective action, it relies on a simple two-parameter model system to capture their influence on the impurity at the level of the free energy and effective mass. Moreover, we invoked an additional approximation by using the random phase approximation. One future perspective would be to consider different model systems with more degrees of freedom, as has been already proposed in Grusdt et al. (2015); Grusdt and Demler . It would also be interesting to extend our study to an impurity in a one-dimensional BEC. In Grusdt et al. (2017b) it has been shown that in this regime Feynman’s method and the RG approach are in excellent agreement within the Fröhlich model, and the observed discrepancies in three dimensions are not present. Finally, at the impurity densities created in current experiments, many-polaron effects are expected to be non-negligible already on a mean-field level Van Loon et al. (2018). The variational path-integral approach has been used to study these effects for Fröhlich polarons in solids Verbist et al. (1991); Klimin et al. (2004) and was applied to the study of the Bose bipolaron within the Fröhlich model Casteels et al. (2013b). Combining the inclusion of the extended interactions and the approximations made in this work with these methods would open a possible avenue towards the study of bipolarons and many-polaron effects in Bose gases beyond the mean-field level with the path-integral formalism.

Acknowledgements.

T. Ichmoukhamedov acknowledges financial support in the form of a Ph.D. fellowship of the Research Foundation - Flanders (FWO), project 1135519N. This research was supported by the University Research Fund (BOF) of the University of Antwerp and by the Flemish Research Foundation (FWO-Vl), projects G.0429.15.N and GOG66.16.N.

Appendix A Derivation of the classical Lagrangian

Starting from expression (2), relative to $E_{0}$ and the first-order energy shift, the Hamiltonian for a single impurity in the BEC can be written as:

[TABLE]

The next step is to introduce position and momentum operators defined by:

[TABLE]

which obey $\left[\hat{Q}_{\mathbf{k}},\hat{P}_{\mathbf{k^{\prime}}}\right]=i\hbar\delta_{\mathbf{k,k^{\prime}}}$ and $\hat{Q}^{\dagger}_{\mathbf{k}}=\hat{Q}_{\mathbf{-k}}$, $\hat{P}^{\dagger}_{\mathbf{k}}=\hat{P}_{\mathbf{-k}}$. Here, $M$ is an arbitrary phonon mass and the frequency $\hbar\omega_{\mathbf{k}}=\epsilon(\mathbf{k})$ corresponds to the Bogoliubov energy dispersion. The Hamiltonian in terms of $\hat{Q}_{\mathbf{k}}$ and $\hat{P}_{\mathbf{k}}$ becomes:

[TABLE]

Two types of diverging terms containing the commutator $\left[\hat{Q}_{\mathbf{k}},\hat{P}_{\mathbf{k}}\right]$ arise in the derivation of (51). The first one corresponds to the ground-state energy of the introduced harmonic oscillators $-\sum_{\mathbf{k}}\hbar\omega_{\mathbf{k}}/2$. The other one arises from the cross terms in the extended interactions and is given by $-\frac{g_{ib}}{2V}\sum_{\mathbf{k}}W_{\mathbf{k,k}}^{(1)}$. This term contains a UV divergence that can not be regularized by taking the cutoff-dependence of $g_{ib}$ into account. Neither of these terms contains the impurity coordinate and we will not include them in further discussion. The classical Hamiltonian corresponding to (51) is obtained by replacing the operators with complex scalar variables that obey $Q_{\mathbf{k}}^{*}=Q_{\mathbf{-k}}$:

[TABLE]

The Legendre transformation:

[TABLE]

results in the classical Lagrangian:

[TABLE]

where the impurity coordinate $\mathbf{r}$ has been introduced. The Lagrangian (54) still has to be written as a function of:

[TABLE]

To simplify the algebra we can multiply the RHS of (55) by $P_{\mathbf{k}}/(2M)$ and perform the summation over $\mathbf{k}$. The two momentum-dependent terms in the Lagrangian can then be compactly written as:

[TABLE]

Next we shall look for an explicit expression for $P_{\mathbf{k}}$. Expression (55) can be equivalently written as:

[TABLE]

Note that for a single impurity $\rho_{\mathbf{k-q}}=\rho_{\mathbf{k}}\rho_{\mathbf{-q}}$. By multiplying (57) with $\frac{V_{\mathbf{k}}^{-1}\rho_{\mathbf{-k}}}{\sqrt{\hbar\omega_{\mathbf{k}}}}$ and performing the summation over $\mathbf{k}$ we find:

[TABLE]

Equation (58) can be algebraically solved to obtain:

[TABLE]

where $\eta=\left(1+\frac{g_{ib}}{V}\sum_{\mathbf{k}}\frac{V_{\mathbf{k}}^{-2}}{\hbar\omega_{\mathbf{k}}}\right)^{-1}$. After substituting (59) into (57) the expression for $P_{\mathbf{k}}$ becomes:

[TABLE]

Finally we can substitute (60) in the Lagrangian:

[TABLE]