Diagrammatic Monte Carlo study of the acoustic and the BEC polaron

We consider two large polaron systems that are described by a Fr\"{o}hlich type of Hamiltonian, namely the Bose-Einstein condensate (BEC) polaron in the continuum and the acoustic polaron in a solid. We present ground-state energies of these two systems calculated with the Diagrammatic Monte Carlo (DiagMC) method and with a Feynman all-coupling approach. The DiagMC method evaluates up to very high order a diagrammatic series for the polaron Green's function. The Feynman all-coupling approach is a variational method that has been used for a wide range of polaronic problems. For the acoustic and BEC polaron both methods provide remarkably similar non-renormalized ground-state energies that are obtained after introducing a finite momentum cutoff. For the renormalized ground-state energies of the BEC polaron, there are relatively large discrepancies between the DiagMC and the Feynman predictions. These differences can be attributed to the renormalization procedure for the contact interaction.


I. INTRODUCTION
By virtue of the Coulomb interaction the presence of a charge carrier in a charged lattice induces a polarization. This effect is well-known from the description of an electron or a hole in a polar or ionic semiconductor. The term polaron was coined by Landau in 1933 1 to denote the quasiparticle comprised of a charged particle coupled to a surrounding polarized lattice. For latticedeformation sizes of the order of the lattice parameter, one refers to the system as a small or Holstein polaron 2,3 . For lattice-deformation sizes that are large compared to the lattice parameter, the lattice can be treated as a continuum. This system is known as a large polaron for which Fröhlich proposed the model Hamiltonian 4 Here, theĉ † k (ĉ k ) are the creation (annihilation) operators of the charge carriers with band mass m and momentum k. The second term in the above Hamiltonian gives the energy of the phonons which carry the polarization. Thereby, the operatorb † k (b k ) creates (annihilates) a phonon with wave vector k and energy ω(k). The last term in Eq. (1) denotes the interaction between the charge carrier and the phonons. A plethora of physical phenomena can be described by the above Fröhlich type of Hamiltonian by varying the dispersion ω(k) and the interaction strength V (q). Fröhlich considered the special situation of longitudinal optical (LO) phonons which are dispersionless ω(k) = ω LO . In the LO limit, the interaction amplitude V (q) in Eq. (1) adopts the form Here, V is the volume of the crystal and α LO the dimensionless coupling parameter: arXiv:1406.6506v1 [cond-mat.quant-gas] 25 Jun 2014 with ε ∞ (ε 0 ) the electronic (static) dielectric constants of the crystal and e the charge of the electron. The Fröhlich polaron which is defined by the Eqs. (1)- (2) and the dispersion ω(k) = ω LO , has no analytical solution. More generally, solutions to the Eq. (1) describe a quasiparticle interacting with a bath of non-interacting bosons with energies ω(k) through the mediation of the interaction V (q). One example is the acoustic polaron which corresponds to the interaction of a charge carrier with acoustic phonons 5 . Another example is the BEC polaron consisting of an impurity atom interacting with the Bogoliubov excitations of an atomic Bose-Einstein condensate (BEC) [6][7][8] . Other examples are an electron on a helium film ("ripplopolaron") 9-11 and a charge carrier in a piezoelectric semiconductor ("piezopolaron") 12 .
Due to the relative simplicity of the model Hamiltonian of Eq. (1) it is an ideal testing ground for conducting comparative studies with various many-body techniques (see for example Refs. 13,14 for an overview). The weak coupling regime (small α LO ) was described by Fröhlich with second-order perturbation theory 4 which is equivalent to the Lee-Low-Pines scheme using a canonical transformation 15 . For the strong coupling regime (large α LO ) Landau and Pekar developed a variational technique which predicts the formation of a bound state of the charge carrier in his self-induced potential 16,17 . Feynman developed a superior all-coupling approach 18,19 which captures all the coupling regimes.
A numerical solution of the Fröhlich Hamiltonian of Eq. (1) with the interaction of Eq. (2) has been proposed in Refs. 20,21 . Thereby, a series expansion for the polaron Green's function was evaluated with the aid of a Diagrammatic Monte Carlo (DiagMC) method. The method is "exact" in the sense that the series expansion is convergent and sign-definite and therefore it can be stochastically evaluated with a controllable error. The polaron's energy is extracted from the asymptotic behavior of its Green's function.
Polaron systems are ideal for comparative studies of many-body techniques. Examples of such studies for the Fermi polaron are reported in Refs. [22][23][24] . For the Fermi polaron, a comparison has been made between the Di-agMC method and the variational technique which includes a limited number of particle-hole excitations. It was demonstrated that a variational one particle-hole calculation is already a good approximation, even for strong interactions between the impurity and the particles in the Fermi sea 23,24 . Recently a comparative study of the neutron polaron has been conducted with quantum Monte Carlo and effective field theories 25 . For the ground-state energy of the Fröhlich polaron of Eqs. (1) and (2) it has been shown in Ref. 20 that Feynman's approach reproduces the DiagMC results to a remarkable accuracy. We have reproduced those numerical results. As can be appreciated from Fig. 1 the deviations between the variational Feynman and DiagMC predictions for the ground-state energies of the Fröhlich polaron, are of the order of a few percent, even for the large cou- pling strengths. It is not clear, however, how accurate the Feynman technique is for polaron systems described by a Hamiltonian of the type of Eq. (1) with alternate dispersions ω(k) and interaction amplitudes V (q). Indeed, Feynman's approach is based on a variational action functional that models the coupling to the phonons by a single phononic degree of freedom with a variationally determined mass and harmonic coupling to the electron. This is a rather natural choice for LO phonons, which are dispersionless. However, it seems intuitively less suitable in situations that the phonons' energies cover a finite range of values. Thornber 26 has argued that in those situations, Feynman's model is unlikely to yield accurate results for the system's dynamical properties, but that the system's ground-state energy can still be captured accurately. To our knowledge, this assertion has not yet been sufficiently confirmed. In order to remedy this situation, in this work we compare polaron groundstate energies calculated with the Feynman variational approach against DiagMC results. This will allow us to test the robustness of the Feynman approach. The two prototypical polaron problems considered in this work are the BEC polaron and the acoustic polaron. These problems have been selected because they highlight complementary aspects. The effect of broadening the range of phonon energies is captured by the acoustic polaron. The BEC polaron problem allows one to additionally cover the issues related to renormalizing V (q).
The structure of this manuscript is as follows. In Sec. II the Hamiltonians for the BEC and acoustic polaron are introduced. In Sects. III A and III B the adopted manybody methods for obtaining the ground-state energies of those Hamiltonians are sketched. Results of the two techniques for the ground-state energies of the BEC and acoustic polaron are contained in Sec. IV.

A. BEC polaron
The Hamiltonian of an impurity immersed in a bath of interacting bosons 8 is given by a sum of two termŝ H =Ĥ B +Ĥ I with, The operatorsâ † k (â k ) create (annihilate) bosons with momentum k, mass m and energy k = 2 k 2 /2m. Further, V is the volume of the system. The operatorsĉ † k (ĉ k ) create (annihilate) the impurity with momentum k and mass m I . The boson-boson and impurity-boson interactions in momentum space are V BB (q) and V IB (q). These potentials are replaced by the pseudopotentials g BB and g IB . These constants are chosen such that the two-body scattering properties in vacuum are correctly reproduced. The sum of all vacuum ladder diagrams, given by the T -matrix, represents all possible ways in which two particles can scatter in vacuum. For zero momentum and frequency the T -matrix is given by T (0): with m r = (1/m I + 1/m) −1 the reduced mass. For lowenergy collisions the first-order Born approximation can be applied to model the boson-boson and boson-impurity collisions. As a result, g IB = 2πa IB 2 mr , with a IB the boson-impurity scattering length and g BB = 4πa BB 2 m , with a BB the boson-boson scattering length.
In the Bogoliubov approximation 27 , the Hamiltonian H B of Eq. (4) is written in the diagonal form where the operatorsb † k (b k ) create (annihilate) Bogoliubov quasi-particles. The quasi-particle vacuum energy is with n = N/V the total density and n 0 = N 0 /V the density of the condensed bosons. The average total particle number N = N is fixed, witĥ and N 0 the number of bosons in the condensate. The collective Bogoliubov excitations obey the dispersion relation At long wavelengths, the spectrum becomes ω(k) = |k|c, which is characteristic of a sound wave with velocity c = n 0 g BB /m. The excitation spectrum is conveniently written in the form with k = |k| and ξ = 1/ √ 2mn 0 g BB the healing length of the Bose condensate.
Application of the Bogoliubov transformation to the impurity partĤ I of Eq. (4) gives 6-8 in which we have defined For g IB = 2πa IB 2 mr a dimensionless coupling constant α IB can be defined 8 The final expression for the Hamiltonian for the BEC polaron is given bŷ Obviously, theĤ BP has the format of a Fröhlich-type of Hamiltonian defined in Eq. (1). When presenting numerical results for the BEC polaron, lengths will be expressed in units of ξ, energies in units of 2 mξ 2 and phonon wave vectors in units of 1/ξ. In this way, all quoted variables are dimensionless. In the numerical calculations, we consider an 6 Li impurity in a Na condensate for which m I /m B = 0.263158 8 .

B. Acoustic polaron
In a crystal with two or more atoms per primitive cell, the dispersion relation ω(k) for the phonons develops acoustic as well as optical branches. The acoustic polaron comprises a charge carrier interacting with the longitudinal acoustic phonons and can be described by the Fröhlich type of Hamiltonian of Eq. (1) with the dispersion ω(k) = sk, with s the sound velocity 5 . For the acoustic polaron, the interaction V AC (q) in the Fröhlich Hamilonian adopts the form 5 : with V the volume of the crystal and α AC a dimensionless coupling parameter. When discussing results concerning the acoustic polaron, lengths will be expressed in units of /(ms), energies in units of ms 2 and phonon wave vectors in units of ms/ . The summations over the phonon momenta | k | have a natural cut-off at the boundary k 0 of the first Brillouin zone. At strong coupling, the Feynman approach to the acoustic polaron predicts the emergence of a self-induced binding potential for the impurity ("self-trapped state"). For a system with both Fröhlich and acoustic phonons, the Feynman approach predicts that the dominant mechanism for this transition is the interaction with the acoustic phonons 28 . Only considering the acoustic phonons results in a transition of the first order for k 0 > 18 and a critical point at k 0 ≈ 18 and α AC ≈ 0.151 5 . This transition was also predicted by the path integral Monte Carlo method 29 .

A. Feynman variational path integral
The Feynman approach is based on the Jensen-Feynman inequality for the free energy F of a system with action S 19 : Here, F 0 is the free energy of a trial system with action S 0 , ... S0 denotes the expectation value with respect to the trial system and β = (k B T ) −1 is the inverse temperature. Feynman proposed a variational trial system of a charge carrier harmonically coupled with spring frequency W to a fictitious particle with mass M . For T = 0 the Jensen-Feynman inequality of Eq. (16) applied to this system produces an upper bound E F p for the pola-ronic ground-state energy 18,19 : with Ω = W 1 + M/m . The function D (k, u) is the phonon Green's function in momentum-imaginary-time representation (k, τ ) where θ(τ ) is the Heaviside function. The memory function M (k, u) is: The u-integral in Eq. (17) is of the following form: The Green's function of the polaron in the (k, τ ) representation is defined as: withĉ the annihilation operator in the Heisenberg representation and |vac the vacuum state. The BEC polaron HamiltonianĤ BP of Eq. (14) contains a vacuum energy E 0 + n 0 g IB which we choose as the zero of the energy scale. Accordingly,Ĥ BP |vac = 0. We define {|ν(k) } as those eigenfunctions ofĤ BP with energy eigenvalue E ν (k) and with one impurity with momentum k. Inserting a complete set of eigenstates in Eq. (21) gives Under the conditions that the polaron is a stable quasiparticle in the ground state (in the sense that it appears as a δ-function peak in the spectral function), one can extract its energy E p (k) and Z-factor Z 0 by studying the long imaginary time behavior of the polaron's Green's function: where the parameter µ is introduced to render a descending exponential tail and with Ψ(k) the fully interacting ground state. The asymptotic behavior of Eq. (24) is associated with a pole singularity for the Green's function in imaginary-frequency representation. For (E p (k) − µ) > 0 one has The one-body self-energy Σ(k, ω) is related to the Green's function by means of the Dyson equation with G 0 (k, ω) the free impurity Green's function. Since the Eqs. (26) and (27) possess the same pole structure, the following expression for the polaronic ground-state energy E p = E p (k = 0) can be obtained 20 : with Σ(τ ) = Σ(0, τ ). Calculating the Green's function boils down to summing a series of Feynman diagrams over all topologies and orders, thereby integrating over all internal variables (like momentum and imaginary time). It is shown in 20 that DiagMC is very suitable to accurately compute the Green's function through a series expansion. We consider irreducible diagrams (an example is shown in Fig. 2) and evaluate a large number of diagrams D in order to numerically compute the Σ(p, τ ) where ξ n represents the topology, n the diagram order, q i is the internal momentum and τ i is the imaginary time. The DiagMC technique allows one to sample over all topologies, all orders and all values of the internal variables. In Fig. 2 some Feynman diagrams for Σ(τ ) are shown. The algebraic expression for these diagrams is given in terms of free propagators and interaction vertices: (i) The free-impurity propagator in imaginary time is determined by (ii) The propagator for an elementary phonon excitation, either of the Bogoliubov type for the BEC polaron, or acoustic phonons for the acoustic polaron is defined in Eq. (18).
(iii) A vertex factor V (q) whenever an elementary excitation carrying momentum q is created or annihilated.
The diagram order is defined by the number of elementary excitations.

A. BEC polaron
For the Fröhlich polaron for which the ground-state energies are displayed in Fig. 1, the one-body selfenergy Σ(τ ) can be computed by means of the procedure sketched in Sec. III B. For the BEC polaron, on the other hand, one encounters ultraviolet divergences when evaluating Σ(τ ) and its energy cannot be extracted. Renormalization/regularization of the impurity-boson pseudopotential is required to obtain physically relevant results for the energies. As a first step in the renormalization procedure, we introduce a momentum cutoff q c upon replacing the momentum summations in Eq. (14) by integrals: This allows us to calculate Σ(τ ) and the accompanying ground-state energy E M C p . From now on we will make the distinction between the polaron energy calculated by DiagMC (E M C p ) and calculated by the Feynman approach (E F p ). Obviously, E M C,F p depends on q c and in order to stress this dependence we use the notation E M C,F p (q c ). In Fig. 3 we show an example of the time dependence of the one-body self-energy Σ(τ ) for the BEC polaron for q c = 200. As can be noticed, after introducing a momentum cutoff q c , the τ dependence is well behaved and the asymptotic regime of Σ(τ ) can be identified. The ∞ n=0 in Eq. (29) implies a summation over an infinite number of diagram orders. In practice, we set a cutoff N max for n in evaluating Σ(τ ). For each N max we can find a corresponding imaginary time τ max for which Σ(τ ) converges. Upon increasing N max we can choose a larger value for τ max . An optimal N max is reached when we can find a τ max in the asymtotic regime that allows us to fit the tail of Σ(τ ). In this way we make an extrapolation for τ → ∞ which determines the value N max . Typical values of N max are of the order 10 4 for large values of α IB . With the aid of the Eq. (28), E M C p (q c ) can be extracted from the computed Σ(τ ). The error on E M C p (q c ) contains a statistical error and a systematic error stemming from the fitting procedure. As can be appreciated from Fig. 3, the grid in imaginary time has to be chosen carefully, since the short-time behavior of Σ(τ ) is strongly peaked. The Σ(p, τ ) for these short times delivers a large contribution to the energy.
In Fig. 4, results for the non-renormalized energies E F p (q c ) and E MC p (q c ) are presented as a function of the dimensionless coupling parameter α IB defined in Eq. (13). The α IB and q c dependence of the DiagMC energies is remarkably similar to those of the Feynman energies. We observe that E MC p (q c ) lies a few percent below E F p (q c ) for all combinations of α IB and q c considered.
In Ref. 8 a renormalization procedure to eliminate the q c dependence of the computed polaron energy is outlined. When determinig the T -matrix of Eq. (5) up to second order, the following relation between the scattering length a IB and the coupling strength g IB is obtained: Using this expression, the n 0 g IB term in Eq. (14) can be replaced by : whereby we have defined E ren (q c ) : This renormalization procedure was developed in the context of the Feynman approach 8 . The same procedure can also be applied in the DiagMC framework. In both frameworks, the renormalized polaron ground-state energy can be found by evaluating the sum In order to illustrate the convergence of the Eq. (35) in both approaches, in Fig. 5 the energies [E M C p (q c ) + E ren (q c )] and [E F p (q c )+E ren (q c )] are plotted as a function of q c for a representative value α IB = 3 of the coupling strength. We notice that the DiagMC and the Feynman approach display an analogous q c dependence. Convergence is reached for q c 3000. Fig. 6 shows that the Feynman path-integral predictions for the BEC-polaron ground-state energies overshoot the DiagMC ones. The relative difference between the two predictions increases with growing values of q c . The very good agreement between the two methods that was found in Fig. 4 for the non-renormalized energies, is no longer observed for the renormalized energies. Indeed, the latter are obtained with Eq. (35), which amounts to substracting two numbers of almost equal magnitude. Accordingly, the final result for the renormalized BEC-polaron ground-state energy is highly sensitive to the adopted many-body technique and renormalization procedure. Fig. 7 illustrates that for small α IB both methods reproduce the result from second-order perturbation theory. The DiagMC method samples diagrams according to their weight and it can be recorded how many times a specific diagram is sampled. In this way, one can identify those diagrams with the largest weight in the self-energy Σ(τ ). At fixed diagram order, we have observed that the number of first-order subdiagrams-the definition of which is explained in the caption of Fig. 8-plays a crucial role in the weight of the diagram. Our studies indicate that for q c > 50 the most important diagram is the one with the highest number of first-order subdiagrams. We have considered many combinations of α IB and q c and could draw this conclusions in all those situations. The dominance of this diagram becomes more explicit with increasing values of q c .
FIG. 8. A diagram of order five for the one-body self-energy. Line conventions as in Fig. 2. Imaginary time runs from left to right. A first-order subdiagram occurs whenever a firstorder diagram drops out from the full diagram by cutting the solid line at two selected times. For example, the considered diagram contains four first-order subdiagrams.

B. Acoustic polaron
We now discuss the numerical results for the groundstate energy of the acoustic polaron introduced in Sec. II B. In Figs  From a detailed analysis of the DiagMC results for k 0 = 50 we find that the class of diagrams of the type sketched in Fig. 8 plays a dominant role for α AC < α c . For α AC > α c we observe a dramatic change in the importance of those diagrams, and we can no longer identify a class of a diagrams that provides the major contribution to the self-energy Σ(τ ).
The knowledge of a certain class of dominant diagrams can be exploited to develop approximate schemes. Indeed, one can set up a self-consistent scheme thereby summing over an important class of diagrams, including the observed dominant ones. In practice, the procedure can be realized by introducing bold (or dressed) propa-gators Σ (i−1) (p, ω) = dω dq (2π) 3 × G (i−1) (p − q, ω − ω )D(q, ω ) with ω and ω the imaginary frequencies. The self-energy Σ (i−1) and the dressed Green's function G (i) (p, ω) are calculated for subsequent values of i, starting from i = 1, until G (i) (p, ω) is converged. In this way Σ (i) (p, ω) will contain all diagrams for which the lines of the phonon propagators do not cross.

V. CONCLUSIONS
We have studied the ground-state energies of the BEC polaron and the acoustic polaron, two large polaron systems that can be described by a Fröhlich type of Hamiltonian. When calculating energies for the BEC polaron with the DiagMC and the Feynman variational technique, we encounter similar ultraviolet divergences. For the acoustic polaron, the ultraviolet regularization is achieved by a hard momentum cutoff which is naturally set at the edge of the first Brillouin zone. In this case, the DiagMC and Feynman predictions for the ground-state energies agree within a few percent. The largest deviation between the predictions of both methods, was found at a coupling strength that marks the transition between a quasifree and a self-trapped state. For the BEC polaron, a more involving two-step renormalization procedure is required. The first step is the introduction of a hard momentum cutoff. In line with the results for the acoustic polaron, the DiagMC and Feynman non-renormalized ground-state energies of the BEC polaron which are produced in this step are remarkably similar. Therefore, one can infer that the Feynman variational method reproduces the "exact" DiagMC non-renormalized polaron ground-state energies at a finite momentum cutoff.
In order to obtain the physical, or renormalized, BECpolaron energies from the non-renormalized ones, an additional procedure is required. Thereby, the contact interaction is renormalized with the aid of the lowest-order correction obtained from the Lippmann-Schwinger equation (34). Despite the fact that the absolute difference between the Feynman and DiagMC BEC-polaron energies remains unaffected by this procedure, the final result for the physical energies displays a large discrepancy.
S.N. Klimin and L.A. Pena-Ardila are gratefully acknowledged. The computational resources (Stevin Supercom-puter Infrastructure) and services used in this work were provided by Ghent University, the Hercules Foundation, and the Flemish Government.