# Diagrammatic Monte Carlo study of the acoustic and the Bose-Einstein condensate polaron

Jonas Vlietinck, Wim Casteels, Kris Van Houcke, Jacques Tempere, Jan

Ryckebusch, Jozef T Devreese

► To cite this version:

Jonas Vlietinck, Wim Casteels, Kris Van Houcke, Jacques Tempere, Jan Ryckebusch, et al.. Diagrammatic Monte Carlo study of the acoustic and the Bose-Einstein condensate polaron. New Journal of Physics, 2015, 17 (3), pp.033023. 10.1088/1367-2630/17/3/033023. [hal-01250087](https://hal.science/hal-01250087)

HAL Id: hal-01250087

https://hal.science/hal-01250087

Submitted on 4 Jan 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Diagrammatic Monte Carlo study of the acoustic and the Bose-Einstein condensate polaron

This content has been downloaded from IOPscience. Please scroll down to see the full text.

2015 New J. Phys. 17 033023

(http://iopscience.iop.org/1367-2630/17/3/033023)

View the table of contents for this issue, or go to the journal homepage for more

Download details:

IP Address: 134.157.80.136

This content was downloaded on 04/01/2016 at 11:03

Please note that terms and conditions apply.
<figure><img src="image_2.png" /></figure>

CrossMark

OPEN ACCESS

RECEIVED
29 August 2014

REVISED
20 January 2015

ACCEPTED FOR PUBLICATION
12 February 2015

PUBLISHED
12 March 2015

Content from this work
may be used under the
terms of the Creative
Commons Attribution 3.0
licence.

Any further distribution of
this work must maintain
attribution to the
author(s) and the title of
the work, journal citation
and DOI.

<figure><img src="image_3.png" /></figure>

PAPER

# Diagrammatic Monte Carlo study of the acoustic and the Bose-Einstein condensate polaron

Jonas Vlietinck¹, Wim Casteels²,³, Kris Van Houcke¹,⁴, Jacques Tempere²,⁵, Jan Ryckebusch¹ and Jozef T Devreese²

¹ Department of Physics and Astronomy, Ghent University, Proeftuinstraat 86, B-9000 Gent, Belgium

² TQC, Universiteit Antwerpen, Universiteitsplein 1, B-2610 Wilrijk, Belgium

³ Laboratoire Matériaux et phénomènes Quantiques, Université Paris Diderot Paris-7, 10 rue Alice Domon et Léonie Duquet, F-75013 Paris, France

⁴ Laboratoire de Physique Statistique, Ecole Normale Supérieure, UPMC, Université Paris Diderot, CNRS, 24 rue Lhomond, F-75231 Paris Cedex 05, France

⁵ Lyman Laboratory of Physics, Harvard University, Cambridge, Massachusetts 02138, USA

E-mail: kvhoucke@lps.ens.fr

Keywords: polaron systems, diagrammatic Monte Carlo, diagrammatic series

## Abstract

We consider two large polaron systems that are described by a Frö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's self-energy. 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.

## 1. 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 lattice-deformation 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]

$$ \hat{H}_{\text{pol}} = \sum_{\mathbf{k}} \frac{\hbar^2 \mathbf{k}^2}{2m} \hat{c}_{\mathbf{k}}^{\dagger} \hat{c}_{\mathbf{k}} + \sum_{\mathbf{k}} \hbar\omega(\mathbf{k}) \hat{b}_{\mathbf{k}}^{\dagger} \hat{b}_{\mathbf{k}} + \sum_{\mathbf{k}, \mathbf{q}} V(\mathbf{q}) \hat{c}_{\mathbf{k}+\mathbf{q}}^{\dagger} \hat{c}_{\mathbf{k}} \hat{b}_{-\mathbf{q}}^{\dagger} + \text{h.c.} \quad (1) $$

Here, the $\hat{c}_{\mathbf{k}}^{\dagger}$ ($\hat{c}_{\mathbf{k}}$) are the creation (annihilation) operators of the charge carriers with band mass $m$ and momentum $\mathbf{k}$. The second term in the above Hamiltonian gives the energy of the phonons which carry the polarization. Thereby, the operator $\hat{b}_{\mathbf{k}}^{\dagger}$ ($\hat{b}_{\mathbf{k}}$) creates (annihilates) a phonon with wave vector $\mathbf{k}$ and energy $\hbar\omega(\mathbf{k})$. The last term in equation (1) denotes the interaction between the charge carrier and the phonons. A
<figure><img src="image_4.png" /><figcaption>Figure 1. Ground-state energies for the Fröhlich polaron are shown as a function of the coupling strength $\alpha_{LO}$ of equation (3). The inset shows the relative difference $\Delta E = \frac{E^{MC} - E^{F}}{E^{MC}}$, with $E^{F}$ ($E^{MC}$) the computed energy from the Feynman (DiagMC) approach.</figcaption></figure>

plethora of physical phenomena can be described by the above Fröhlich type of Hamiltonian by varying the dispersion $\omega(\mathbf{k})$ and the interaction strength $V(\mathbf{q})$. Fröhlich considered the special situation of longitudinal optical (LO) phonons which are dispersionless $\omega(\mathbf{k}) = \omega_{LO}$. In the LO limit, the interaction amplitude $V(\mathbf{q})$ in equation (1) adopts the form

$$V_{LO}(\mathbf{q}) = -i \frac{\hbar \omega_{LO}}{q} \left( \frac{4\pi\alpha_{LO}}{V} \right)^{1/2} \left( \frac{\hbar}{2m\omega_{LO}} \right)^{1/4} . \qquad (2)$$

Here, $V$ is the volume of the crystal and $\alpha_{LO}$ the dimensionless coupling parameter:

$$\alpha_{LO} = \frac{e^2}{\hbar} \sqrt{\frac{m}{2\hbar\omega_{LO}}} \left( \frac{1}{\varepsilon_{\infty}} - \frac{1}{\varepsilon_0} \right), \qquad (3)$$

with $\varepsilon_{\infty}$ ($\varepsilon_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 equations (1)–(2) and the dispersion $\omega(\mathbf{k}) = \omega_{LO}$, has no analytical solution.

More generally, solutions to the equation (1) describe a quasiparticle interacting with a bath of non-interacting bosons with energies $\hbar\omega(\mathbf{k})$ through the mediation of the interaction $V(\mathbf{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–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 equation (1) it is an ideal testing ground for conducting comparative studies with various many-body techniques (see for example [13, 14] for an overview). The weak coupling regime (small $\alpha_{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 $\alpha_{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 an approximation scheme [18, 19] which should capture all the coupling regimes.

A numerical solution of the Fröhlich Hamiltonian of equation (1) with the interaction of equation (2) has been proposed in [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 [22–24]. For the Fermi polaron, a comparison has been made between the DiagMC 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 equations (1) and (2) it has been shown in [20] that Feynman’s approach reproduces the DiagMC results to a remarkable accuracy. We have reproduced those numerical results. As can be appreciated from figure 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 coupling strengths. It is not clear, however, how accurate the Feynman technique is for polaron
systems described by a Hamiltonian of the type of equation (1) with alternate dispersions $\omega(\mathbf{k})$ and interaction amplitudes $V(\mathbf{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 ground-state 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(\mathbf{q})$.

The structure of this manuscript is as follows. In section 2 the Hamiltonians for the BEC and acoustic polaron are introduced. In sections 3.1 and 3.2 the adopted many-body 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 section 4.

## 2. Large polaron models

### 2.1. BEC polaron

The Hamiltonian of an impurity immersed in a bath of interacting bosons [8] is given by a sum of two terms $\hat{H} = \hat{H}_B + \hat{H}_I$ with

$$ \begin{aligned} \hat{H}_B &= \sum_{\mathbf{k}} \epsilon_{\mathbf{k}} \hat{a}_{\mathbf{k}}^{\dagger} \hat{a}_{\mathbf{k}} + \frac{1}{2V} \sum_{\mathbf{k}, \mathbf{k}', \mathbf{q}} V_{BB}(\mathbf{q}) \hat{a}_{\mathbf{k}^{-}\mathbf{q}}^{\dagger} \hat{a}_{\mathbf{k}^{+}\mathbf{q}} \hat{a}_{\mathbf{k}} \hat{a}_{\mathbf{k}'} \\ \hat{H}_I &= \sum_{\mathbf{k}} \frac{\hbar^2 \mathbf{k}^2}{2m_I} \hat{c}_{\mathbf{k}}^{\dagger} \hat{c}_{\mathbf{k}} + \frac{1}{V} \sum_{\mathbf{k}, \mathbf{k}', \mathbf{q}} V_{IB}(\mathbf{q}) \hat{c}_{\mathbf{k}+\mathbf{q}}^{\dagger} \hat{c}_{\mathbf{k}} \hat{a}_{\mathbf{k}^{-}\mathbf{q}}^{\dagger} \hat{a}_{\mathbf{k}}. \end{aligned} \quad (4) $$

The operators $\hat{a}_{\mathbf{k}}^{\dagger}(\hat{a}_{\mathbf{k}})$ create (annihilate) bosons with momentum $\mathbf{k}$, mass $m$ and energy $\epsilon_{\mathbf{k}} = \hbar^2 \mathbf{k}^2 / 2m$. Further, $V$ is the volume of the system. The operators $\hat{c}_{\mathbf{k}}^{\dagger}(\hat{c}_{\mathbf{k}})$ create (annihilate) the impurity with momentum $\mathbf{k}$ and mass $m_I$. The boson-boson and impurity-boson interactions in momentum space are $V_{BB}(\mathbf{q})$ and $V_{IB}(\mathbf{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)$:

$$ T(0) = g_{IB} - g_{IB} \sum_{\mathbf{k}} \frac{2m_r}{\hbar^2 \mathbf{k}^2} T(0), \quad (5) $$

with $m_r = (1/m_I + 1/m)^{-1}$ the reduced mass. For low-energy collisions the first-order Born approximation can be applied to model the boson-boson and boson-impurity collisions. As a result, $g_{IB} = \frac{2\pi a_{IB} \hbar^2}{m_r}$, with $a_{IB}$ the boson-impurity scattering length and $g_{BB} = \frac{4\pi a_{BB} \hbar^2}{m}$, with $a_{BB}$ the boson-boson scattering length.

In the Bogoliubov approximation [27], the Hamiltonian $\hat{H}_B$ of equation 4 is written in the diagonal form

$$ \hat{H}_B \approx E_0 + \sum_{\mathbf{k} \neq 0} \hbar \omega(\mathbf{k}) \hat{b}_{\mathbf{k}}^{\dagger} \hat{b}_{\mathbf{k}}, \quad (6) $$

where the operators $\hat{b}_{\mathbf{k}}^{\dagger}(\hat{b}_{\mathbf{k}})$ create (annihilate) Bogoliubov quasi-particles. The quasi-particle vacuum energy is

$$ E_0 = \frac{V}{2} n^2 g_{BB} + \frac{1}{2} \sum_{\mathbf{k} \neq 0} \left( \hbar \omega(\mathbf{k}) - \epsilon_{\mathbf{k}} - n_0 g_{BB} \right), \quad (7) $$

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 = \langle \hat{N} \rangle$ is fixed, with

$$ \hat{N} = N_0 + \sum_{\mathbf{k} \neq 0} \hat{a}_{\mathbf{k}}^{\dagger} \hat{a}_{\mathbf{k}}, \quad (8) $$
and N₀ the number of bosons in the condensate. The collective Bogoliubov excitations obey the dispersion relation

$$ \hbar\omega(\mathbf{k}) = \sqrt{(\epsilon_{\mathbf{k}} + n_0 g_{BB})^2 - (n_0 g_{BB})^2}. \quad (9) $$

At long wavelengths, the spectrum becomes $\omega(\mathbf{k}) = |\mathbf{k}|c$, which is characteristic of a sound wave with velocity $c = \sqrt{n_0 g_{BB}/m}$. The excitation spectrum is conveniently written in the form

$$ \omega(\mathbf{k}) = k c \sqrt{1 + \frac{(k \xi)^2}{2}}, \quad (10) $$

with $k = |\mathbf{k}|$ and $\xi = 1/\sqrt{2mn_0g_{BB}}$ the healing length of the Bose condensate.

Application of the Bogoliubov transformation to the impurity part $\hat{H}_I$ of equation (4) gives [6–8]

$$ \hat{H}_I \approx \sum_{\mathbf{k}} \frac{\hbar^2 \mathbf{k}^2}{2m_I} \hat{c}_{\mathbf{k}}^{\dagger} \hat{c}_{\mathbf{k}} + n_0 g_{IB} + \sum_{\mathbf{q} \neq 0, \mathbf{k}} V_{BP}(\mathbf{q}) \hat{c}_{\mathbf{k}+\mathbf{q}}^{\dagger} \hat{c}_{\mathbf{k}} (\hat{b}_{-\mathbf{q}}^{\dagger} + \hat{b}_{\mathbf{q}}), \quad (11) $$

in which we have defined

$$ V_{BP}(\mathbf{q}) = \frac{g_{IB}}{\mathcal{V}} \sqrt{\frac{N_0 \epsilon_{\mathbf{q}}}{\omega(\mathbf{q})}} = \frac{g_{IB} \sqrt{N_0}}{\mathcal{V}} \left( \frac{(\xi q)^2}{(\xi q)^2 + 2} \right)^{1/4}, \quad (12) $$

with $g_{IB} = \frac{2\pi a_{IB} \hbar^2}{m_I}$. A dimensionless coupling constant $\alpha_{IB}$ can be defined [8]

$$ \alpha_{IB} = \frac{a_{IB}^2}{a_{BB}\xi}. \quad (13) $$

The final expression for the Hamiltonian for the BEC polaron is given by

$$ \hat{H}_{\text{BP}} = E_0 + n_0 g_{\text{IB}} + \sum_{\mathbf{k}} \frac{\hbar^2 \mathbf{k}^2}{2m_I} \hat{c}_{\mathbf{k}}^{\dagger} \hat{c}_{\mathbf{k}} + \sum_{\mathbf{k} \neq 0} \hbar \omega(\mathbf{k}) (\hat{b}_{\mathbf{k}}^{\dagger} \hat{b}_{\mathbf{k}} + \sum_{\mathbf{q} \neq 0, \mathbf{k}} V_{\text{BP}}(\mathbf{q}) \hat{c}_{\mathbf{k}+\mathbf{q}}^{\dagger} \hat{c}_{\mathbf{k}} (\hat{b}_{-\mathbf{q}}^{\dagger} + \hat{b}_{\mathbf{q}})). \quad (14) $$

Obviously, the $\hat{H}_{\text{BP}}$ has the format of a Fröhlich-type of Hamiltonian defined in equation (1). When presenting numerical results for the BEC polaron, lengths will be expressed in units of $\xi$, energies in units of $\frac{\hbar^2}{m\xi^2}$ and phonon wave vectors in units of $1/\xi$. 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].

The Fröhlich Hamiltonian of equation (14) provides an effective low-energy description of an impurity atom in a BEC. Its accurateness depends on the validity of the Bogoliubov approximation. In the limit of weak polaronic coupling, this approximation is expected to be accurate.

At strong polaronic coupling a collapse of the impurity wave function is expected for an attractive impurity–boson interaction [28], making the polaron picture inaccurate. For a repulsive impurity–boson interaction, the local depletion of the boson density around the impurity can render the Bogoliubov approximation invalid. This results in the formation of a bubble state, reminiscent of an electron in a condensed helium superfluid [29]. It was shown in [29] that for a repulsive impurity–boson interaction the Fröhlich Hamiltonian remains valid in the strong coupling regime as long as the mass scaled gas parameter $\frac{m}{m_I} \sqrt{n_0 a_{BB}^3}$ is sufficiently small. In the following a system with a repulsive impurity–boson interaction within the validity regime of the Fröhlich Hamiltonian is considered.

## 2.2. Acoustic polaron

In a crystal with two or more atoms per primitive cell, the dispersion relation $\omega(\mathbf{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 equation (1) with the dispersion $\omega(\mathbf{k}) = sk$, with s the sound velocity [5]. For the acoustic polaron, the interaction $V_{AC}(\mathbf{q})$ in the Fröhlich Hamiltonian adopts the form [5]:

$$ V_{AC}(\mathbf{q}) = \left( \frac{4\pi\alpha_{AC}}{\mathcal{V}} \right)^{1/2} \frac{\hbar^2}{m} \sqrt{\mathbf{q}}, \quad (15) $$

with $\mathcal{V}$ the volume of the crystal and $\alpha_{AC}$ a dimensionless coupling parameter. When discussing results concerning the acoustic polaron, lengths will be expressed in units of $\hbar/(\text{ms})$, energies in units of $\text{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₀ 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 [30]. Only considering the acoustic phonons results in a transition of the first order for k₀ > 18 and a critical point at k₀ ≈ 18 and αₐₐₐ ≈ 0.151 [5]. This transition was also predicted by the path integral Monte Carlo method [31].

### 3. Numerical methods

#### 3.1. 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]:

$$ \mathcal{F} \le \mathcal{F}_0 + \frac{1}{\hbar\beta} \langle S - S_0 \rangle_{S_0}. \qquad (16) $$

Here, $\mathcal{F}_0$ is the free energy of a trial system with action $S_0$, $\langle \dots \rangle_{S_0}$ denotes the expectation value with respect to the trial system and $\beta = (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 equation (16) applied to this system produces an upper bound $E_p^F$ for the polaronic ground-state energy [18, 19]:

$$ E_p \le \frac{3\hbar\Omega}{4} \frac{\left(\sqrt{\left(1 + M/m_I\right)} - 1\right)^2}{1 + M/m_I} + \sum_k \frac{|V_k|^2}{\hbar} \int_0^\infty du \mathcal{D}(\mathbf{k}, u) \mathcal{M}(\mathbf{k}, u), \qquad (17) $$

with $\Omega = W\sqrt{1 + M/m}$. The function $\mathcal{D}(\mathbf{k}, u)$ is the phonon Green's function in momentum-imaginary-time representation $(\mathbf{k}, \tau)$

$$ \mathcal{D}(\mathbf{k}, \tau) = -\theta(\tau) \exp[-\omega(\mathbf{k})\tau], \qquad (18) $$

where $\theta(\tau)$ is the Heaviside function. The memory function $\mathcal{M}(\mathbf{k}, u)$ is:

$$ \mathcal{M}(\mathbf{k}, u) = \exp\left[-\frac{\hbar k^2}{2(m_I + M)} \left(u + \frac{M}{m_I} \frac{1 - \exp[-\Omega u]}{\Omega}\right)\right]. \qquad (19) $$

The u-integral in equation (17) is of the following form:

$$ \int_0^\infty du \exp[-au + be^{-u}] = (-b)^{-a} \Gamma(a, -b, 0), \qquad (20) $$

with $\Gamma(a, z_0, z_1) = \int_{z_0}^{z_1} t^{a-1}e^{-t}dt$ the generalized incomplete gamma function. The parameters M and $\Omega$ are used to minimize the upper bound for the ground state energy of equation (17). This approach captures the different coupling regimes.

#### 3.2. One-body propagator and DiagMC

The Green's function of the polaron in the $(\mathbf{k}, \tau)$ representation is defined as:

$$ G(\mathbf{k}, \tau) = -\theta(\tau) \langle \text{vac} | \hat{c}_{\mathbf{k}}(\tau) \hat{c}_{\mathbf{k}}^{\dagger}(0) | \text{vac} \rangle, \qquad (21) $$

with

$$ \hat{c}_{\mathbf{k}}(\tau) = e^{\hat{H}_{\text{BP}}} \hat{c}_{\mathbf{k}} e^{-\hat{H}_{\tau}}, \qquad (22) $$

the annihilation operator in the Heisenberg representation and $|\text{vac}\rangle$ the vacuum state. The BEC polaron Hamiltonian $\hat{H}_{\text{BP}}$ of equation (14) contains a vacuum energy $E_0 + n_0 g_{\text{IB}}$ which we choose as the zero of the energy scale. Accordingly, $\hat{H}_{\text{BP}}|\text{vac}\rangle = 0$. We define $\{|\nu(\mathbf{k})\rangle\}$ as those eigenfunctions of $\hat{H}_{\text{BP}}$ with energy eigenvalue $E_\nu(\mathbf{k})$ and with one impurity with momentum $\mathbf{k}$. Inserting a complete set of eigenstates in equation (21) gives

$$ G(\mathbf{k}, \tau) = -\theta(\tau) \sum_{\nu} |\langle \nu(\mathbf{k}) | \hat{c}_{\mathbf{k}}^{\dagger} | \text{vac} \rangle|^2 e^{-E_{\nu}(\mathbf{k})\tau}. \qquad (23) $$

Under the conditions that the polaron is a stable quasi-particle in the ground state (in the sense that it appears as a $\delta$-function peak in the spectral function), one can extract its energy $E_p(\mathbf{k})$ and Z-factor $Z_0$ by studying the long
<figure><img src="image_5.png" /><figcaption>Figure 2. Irreducible diagrams for the polaron’s self-energy $\Sigma(\mathbf{k}, \tau, \mu)$. Imaginary time runs from left to right. A solid line represents a free-impurity propagator and a dashed line stands for an elementary excitation. The interaction vertices are denoted by dots.</figcaption></figure>

imaginary time behavior of the polaron’s Green’s function:

$$G(\mathbf{k}, \tau, \mu) \equiv G(\mathbf{k}, \tau) e^{\mu\tau} \stackrel{\tau \to +\infty}{\sim} -Z_0(\mathbf{k}) e^{-(E_p(\mathbf{k})-\mu)\tau}, \quad (24)$$

where

$$Z_0(\mathbf{k}) = |\langle \Psi(\mathbf{k}) | \hat{c}_{\mathbf{k}}^{\dagger} | \text{vac} \rangle|^2, \quad (25)$$

with $\Psi(\mathbf{k})$ the fully interacting ground state. The unphysical parameter $\mu$ is introduced to control the exponential tail of $G$ in imaginary time, and to ensure that it is always descending. The particular choice of $\mu$ has no impact on the final results. The asymptotic behavior of equation (24) is associated with a pole singularity for the Green’s function in imaginary-frequency representation. For $(E_p(\mathbf{k}) - \mu) > 0$ one has

$$G(\mathbf{k}, \omega, \mu) = \int_{0}^{+\infty} d\tau e^{i\omega\tau} G(\mathbf{k}, \tau, \mu) = \frac{Z_0(\mathbf{k})}{i\omega + \mu - E_p(\mathbf{k})} + \text{regular part.} \quad (26)$$

The one-body self-energy $\Sigma(\mathbf{k}, \omega, \mu)$ is related to the Green’s function by means of the Dyson equation

$$G(\mathbf{k}, \omega, \mu) = \frac{1}{\frac{1}{G^0(\mathbf{k}, \omega, \mu)} - \Sigma(\mathbf{k}, \omega, \mu)}, \quad (27)$$

with $G^0(\mathbf{k}, \omega, \mu)$ the free impurity Green’s function (see equation (31)). Since the equations (26) and (27) possess the same pole structure, the following expression for the polaronic ground-state energy $E_p = E_p(\mathbf{k} = \mathbf{0})$ can be obtained [20]:

$$E_p = \Sigma(\mathbf{k} = \mathbf{0}, \omega = 0, \mu = E_p). \quad (28)$$

The $\mu$ dependence of the self-energy in the imaginary-time representation $\Sigma(\mathbf{k}, \tau, \mu)$ adopts a simple form

$$\Sigma(\mathbf{k}, \tau, \mu) = \frac{1}{2\pi} \int_{0}^{\infty} d\omega e^{-i\omega\tau} \Sigma(\mathbf{k}, \omega, \mu) = \Sigma(\mathbf{k}, \tau, \mu = 0) e^{\mu\tau}. \quad (29)$$

With the above expressions one immediately sees that

$$E_p = \int_0^\infty d\tau \Sigma(\tau, \mu = E_p) = \int_0^\infty d\tau \Sigma(\tau, \mu) e^{(E_p - \mu)\tau}, \quad (30)$$

with $\Sigma(\tau, \mu) \equiv \Sigma(\mathbf{0}, \tau, \mu)$. The $\tau$ dependence of the quantity $\Sigma(\tau, \mu)$ at a fixed value of $\mu$ can be numerically evaluated in DiagMC. Thereby, it is beneficial to perform the calculation at a value of $\mu$ approaching the magnitude of $E_p$. Indeed, from equation (30) it is clear that for $\mu \lesssim E_p$ the statistical noise at large values of $\tau$ in $\Sigma(\tau, \mu)$ can be better kept under control. Although $E_p$ is a priori unknown, its value can be estimated in an early stage of the simulation. After performing the DiagMC simulation, equation (30) can be solved iteratively.

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.

In figure (2) some Feynman diagrams for $\Sigma$ 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

$$G^0(\mathbf{k}, \tau, \mu) = -\theta(\tau) e^{-(\epsilon_k - \mu)\tau}. \quad (31)$$

(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 equation (18).

(iii) A vertex factor $V(\mathbf{q})$ whenever an elementary excitation carrying momentum $\mathbf{q}$ is created or annihilated. We consider irreducible diagrams and evaluate a large number of diagrams $D$ in order to numerically compute the $\Sigma(\mathbf{k}, \tau, \mu)$
<figure><img src="image_6.png" /><figcaption>**Figure 3.** The one-body self-energy $\Sigma(\tau, \mu)$ for $\mu = -790$ for the BEC polaron plotted as a function of imaginary time $\tau$. Results are obtained for $\alpha_{IB} = 5$ and $q_c = 200$ and exclude the first-order contribution to $\Sigma(\tau, \mu)$ which can be easily computed analytically. The inset shows $\Sigma(\tau, \mu)$ for small imaginary times.</figcaption></figure>

$$ \Sigma(\mathbf{k}, \tau, \mu) = \sum_{n=1}^{\infty} \sum_{\xi_n} \sum_{\mathbf{q}_{1}=1,...,n} \int \int \int d\tau_1...d\tau_i...d\tau_{2n-2} \\ \times D(\xi_n, \mathbf{k}, \tau, \mu, \tau_1, ..., \tau_i, ..., \tau_{2n-2}, \mathbf{q}_1, ..., \mathbf{q}_i, ..., \mathbf{q}_n), \quad (32) $$

where $\xi_n$ represents the topology, $n$ the diagram order, $\mathbf{q}_i$ are the independent internal momentum and $\tau_i$ is the internal imaginary times. We define the diagram order by counting the number of elementary excitations (the number of dashed propagators in figure 2). The DiagMC technique allows one to sample over all topologies, all orders and all values of the internal variables, and thus to determine $\Sigma$.

# 4. Results and discussion

## 4.1. BEC polaron

For the Fröhlich polaron for which the ground-state energies are displayed in figure 1, the one-body self-energy $\Sigma(\tau, \mu)$ can be computed by means of the procedure sketched in section 3.2. For the BEC polaron, on the other hand, one encounters ultraviolet divergences when evaluating $\Sigma(\tau, \mu)$ 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 equation (14) by integrals:

$$ \sum_{\mathbf{k}} \rightarrow \frac{\nu}{(2\pi)^3} \int_{|\mathbf{k}|<q_c} \mathrm{d}\mathbf{k}. \quad (33) $$

This allows us to calculate $\Sigma(\tau, \mu)$ and the accompanying ground-state energy $E_p^{MC}$. From now on we will make the distinction between the polaron energy calculated by DiagMC ($E_p^{MC}$) and calculated by the Feynman approach ($E_p^F$). Obviously, $E_p^{MC,F}$ depends on $q_c$ and in order to stress this dependence we use the notation $E_p^{MC,F}(q_c)$. In figure 3 we show an example of the time dependence of the one-body self-energy $\Sigma(\tau, \mu)$ for the BEC polaron for $q_c = 200$. As can be noticed, after introducing a momentum cutoff $q_c$, the $\tau$ dependence is well behaved and the asymptotic regime of $\Sigma(\tau, \mu)$ can be identified. The $\sum_{n=0}^{\infty}$ in equation (32) implies a summation over an infinite number of diagram orders. In practice, we set a cutoff $N_{max}$ for $n$ in evaluating $\Sigma(\tau, \mu)$. For each $N_{max}$ we can find a corresponding imaginary time $\tau_{max}$ for which we observe that $\sum_{n=0}^{N_{max}} \Sigma^{(n)}(\tau \le \tau_{max}, \mu) = \sum_{n=0}^{\infty} \Sigma^{(n)}(\tau \le \tau_{max}, \mu)$ within the numerical noise. Hereby, $\Sigma^{(n)}(\tau, \mu)$ is the contribution from the $n$th order diagrams to the self energy. Upon increasing $N_{max}$ the value of $\tau_{max}$ increases accordingly. An optimal $N_{max}$ is reached when we can find a $\tau_{max}$ in the asymptotic regime that allows us to fit the tail of $\Sigma(\tau, \mu)$. In this way we make an extrapolation for $\tau \to \infty$. Typical values of $N_{max}$ are of the order $10^4$ for large values of $\alpha_{IB}$. With the aid of the equation (30), $E_p^{MC}(q_c)$ can be extracted from the computed $\Sigma(\tau, \mu)$. The error on $E_p^{MC}(q_c)$ contains a statistical error and a systematic error stemming from the fitting procedure. As can be appreciated from figure 3, the grid in imaginary time has to be chosen carefully, since the short-time behavior of $\Sigma(\tau, \mu)$ is strongly peaked. The $\Sigma(p, \tau, \mu)$ for these short times delivers a large contribution to the energy.

In figure 4, results for the non-renormalized energies $E_p^F(q_c)$ and $E_p^{MC}(q_c)$ are presented as a function of the dimensionless coupling parameter $\alpha_{IB}$ defined in equation (13). The $\alpha_{IB}$ and $q_c$ dependence of the DiagMC
<figure><img src="image_7.png" /><figcaption>Figure 4. The non-renormalized BEC-polaron energy $E_p$ as a function of the coupling strength $\alpha_{IB}$ as computed with the DiagMC (symbols) and with the Feynman (lines) approaches. Results are shown for four values of the cutoff momentum.</figcaption></figure>

energies is remarkably similar to those of the Feynman energies. We observe that $E_p^{MC}(q_c)$ lies a few percent below $E_p^F(q_c)$ for all combinations of $\alpha_{IB}$ and $q_c$ considered.

In [8] a renormalization procedure to eliminate the $q_c$ dependence of the computed polaron energy is outlined. When determining the T-matrix of equation (5) up to second order, the following relation between the scattering length $a_{IB}$ and the coupling strength $g_{IB}$ is obtained:

$$ \frac{2\pi a_{IB} \hbar^2}{m_r} = g_{IB} - \frac{g_{IB}^2}{(2\pi)^3} \int_{|q|<q_c} dq \frac{2m_r}{\hbar^2 q^2}. \qquad (34) $$

Using this expression, the $n_0 g_{IB}$ term in equation (14) can be replaced by:

$$ n_0 g_{IB} \rightarrow \frac{2\pi a_{IB} n_0 \hbar^2}{m_r} + E_{ren}(q_c), \qquad (35) $$

whereby we have defined $E_{ren}(q_c)$:

$$ E_{ren}(q_c) = \frac{n_0 g_{IB}^2}{(2\pi)^3} \int_{|q|<q_c} dq \frac{2m_r}{\hbar^2 q^2}. \qquad (36) $$

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

$$ E_p^{MC,F} = E_p^{MC,F}(q_c \rightarrow \infty) + E_{ren}(q_c \rightarrow \infty). \qquad (37) $$

This renormalization procedure preserves the absolute difference between the $E_p^{MC}(q_c \rightarrow \infty)$ and the $E_p^F(q_c \rightarrow \infty)$ results.

In order to illustrate the convergence of the equation (37) in both approaches, in figure 5 the energies $[E_p^{MC}(q_c) + E_{ren}(q_c)]$ and $[E_p^F(q_c) + E_{ren}(q_c)]$ are plotted as a function of $q_c$ for a representative value $\alpha_{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 \gtrsim 3000$. As a matter of fact, the $E_p^{MC,F}(q_c)$ and $E_{ren}(q_c)$ are of equal magnitude. This imposes severe constraints on the required accuracy of the used many-body technique. Small deviations in the computed $E_p^{MC,F}(q_c)$ can result in large changes in the obtained values of the renormalized polaron ground-state energy. This explains the relatively large differences (more than an order of magnitude for the results of figure 5) between the unbiased DiagMC results $[E_p^{MC}(q_c) + E_{ren}(q_c)]$ and the approximate $[E_p^F(q_c) + E_{ren}(q_c)]$ ones.

Figure 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 figure 4 for the non-renormalized energies, is no longer observed for the renormalized energies. Indeed, the latter are obtained with equation (37), which amounts to subtracting 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. Figure 7 illustrates that for small $\alpha_{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-
<figure><img src="image_8.png" /><figcaption>**Figure 5.** The renormalized BEC-polaron energies [$E_p^{MC}(q_c) + E_{ren}(q_c)$] at $\alpha_{IB} = 3$ are given as a function of the momentum cutoff $q_c$. The inset figure shows [$E_p^F(q_c) + E_{ren}(q_c)$] as a function of $q_c$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** The renormalized BEC-polaron energies [$E_p(q_c) + E_{ren}(q_c)$] as a function of $\alpha_{IB}$ for different values of the momentum cutoff $q_c$. Lines are the Feynman path-integral and symbols are the DiagMC results.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 7.** The renormalized BEC-polaron energies [$E_p(q_c) + E_{ren}(q_c)$] at small values of $\alpha_{IB}$ at $q_c = 2000$. The dot-dashed line is the Feynman path-integral result, symbols represent the DiagMC results, while the short dashed line is the prediction from second-order perturbation theory (PT).</figcaption></figure>

energy $\Sigma(\tau, \mu)$. At fixed diagram order, we have observed that the number of first-order subdiagrams—the definition of which is explained in the caption of figure 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 $\alpha_{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$.

## 4.2. Acoustic polaron

We now discuss the numerical results for the ground-state energy of the acoustic polaron introduced in section 2.2. In figures 9 and 10 we show a selection of the predictions $E_p^F$ from the Feynman upper-bound
<figure><img src="image_11.png" /><figcaption>**Figure 8.** A diagram of order five for the one-body self-energy. Line conventions as in figure 2. Imaginary time runs from left to right. A first-order subdiagram occurs whenever a first-order 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.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Non-renormalized ground-state energies $E_p^F$ and $E_p^{MC}$ for the acoustic polaron as a function of $\alpha_{AC}$ for $k_0 = 10$. The inset shows $\Delta E = \frac{E_p^{MC} - E_p^F}{E_p^{MC}}$ as a function of $\alpha_{AC}$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 10.** As in figure 9 but for $k_0 = 50$. The vertical dashed line denotes the coupling strength $\alpha_{AC} = 0.052$ corresponding with the transition as computed in [5].</figcaption></figure>

method of equation (17) together with the DiagMC results $E_p^{MC}$ which are computed with the aid of equation (30). For $k_0 = 10$ and $k_0 = 50$ an excellent agreement between $E_p^F$ and $E_p^{MC}$ is found. From the relative difference $\Delta E = \frac{E_p^{MC} - E_p^F}{E_p^{MC}}$, a value $\alpha_{AC}$ can be found where $\Delta E$ is largest in the considered region of $\alpha_{AC}$. For $k_0 = 10$ we find $\alpha_{AC}^{k_0=10} = 0.28 \pm 0.04$ and for $k_0 = 50$, $\alpha_{AC}^{k_0=50} = 0.52 \pm 0.01$. For $\alpha < \alpha_c$, $\Delta E$ increases with $\alpha_{AC}$ and for $\alpha > \alpha_c$, $\Delta E$ decreases with increasing $\alpha_{AC}$. We remark that $\alpha_c^{k_0=10}$ and $\alpha_c^{k_0=50}$ coincides with the coupling strength for the transition [30] as computed with the Feynman approach.

From a detailed analysis of the DiagMC results for $k_0 = 50$ we find that the class of diagrams of the type sketched in figure 8 plays a dominant role for $\alpha_{AC} < \alpha_c$. For $\alpha_{AC} > \alpha_c$ we observe a dramatic change in the importance of those diagrams, and we can no longer identify a class of diagrams that provides the major contribution to the self-energy $\Sigma(\tau, \mu)$.

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) propagators
$$ \Sigma^{(i-1)}(\mathbf{p}, \omega, \mu) = \int d\omega' \int \frac{d\mathbf{q}}{(2\pi)^3} G^{(i-1)}(\mathbf{p}-\mathbf{q}, \omega-\omega', \mu) D(\mathbf{q}, \omega') \\ \times G^{(i)}(\mathbf{p}, \omega, \mu) = \frac{1}{G^{0-1}(\mathbf{p}, \omega, \mu) - \Sigma^{(i-1)}(\mathbf{p}, \omega, \mu)} $$

with $\omega$ and $\omega'$ the imaginary frequencies. The self-energy $\Sigma^{(i-1)}$ and the dressed Green's function $G^{(i)}(\mathbf{p}, \omega, \mu)$ are calculated for subsequent values of $i$, starting from $i=1$, until $G^{(i)}(\mathbf{p}, \omega, \mu)$ is converged. In this way $\Sigma^{(i)}(\mathbf{p}, \omega, \mu)$ will contain all diagrams for which the lines of the phonon propagators do not cross.

## 5. 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, BEC-polaron energies from the non-renormalized ones, an additional procedure is required. Thereby, the ultraviolet behavior of the contact interaction is renormalized with the aid of the lowest-order correction obtained from the Lippmann-Schwinger equation (36). In the regime of strong coupling this results in significantly different Feynman and DiagMC BEC-polaron renormalized energies.

## Acknowledgments

This work is supported by the Flemish Research Foundation (FWO Vlaanderen) through project numbers G.0119.12N and G.0115.12N. Discussions with SN Klimin and LA Pena-Ardila are gratefully acknowledged. The computational resources (Stevin Supercomputer Infrastructure) and services used in this work were provided by Ghent University, the Hercules Foundation, and the Flemish Government.

## References

[1] Landau L D 1933 Phys. Z. Sowjetunion 3 664

[2] Wellein G, Röder H and Fehske H 1996 Phys. Rev. B 53 9666

[3] Marchand D J J and Berciu M 2013 Phys. Rev. B 88 6060301

[4] Fröhlich H 1954 Adv. Phys. 3 325

[5] Peeters F M and Devreese J T 1985 Phys. Rev. B 32 3515

[6] Cucchietti F M and Timmermans E 2006 Phys. Rev. Lett. 96 210401

[7] Sacha K and Timmermans E 2006 Phys. Rev. A 73 063604

[8] Tempere J, Casteels W, Oberthaler M K, Knoop S, Timmermans E and Devreese J T 2009 Phys. Rev. B 80 184504

[9] Shikin V and Monarkha Y 1974 J. Low Temp. Phys. 16 193

[10] Jackson S A and Platzman P M 1981 Phys. Rev. B 24 499

[11] Marques G E and Studart N 1989 Phys. Rev. B 39 4133

[12] Mahan G D and Hopfield J J 1964 Phys. Rev. Lett. 12 241

[13] Devreese J T and Alexandrov A S 2010 Advances In Polaron Physics vol 159 (Berlin: Springer)

[14] Devreese J T 2010 arXiv:1012.4576

[15] Lee T D, Low F E and Pines D 1953 Phys. Rev. 90 297

[16] Landau L D and Pekar S I 1946 Zh. Eksp. Teor. Fiz. 16 341

[17] Pekar S I 1963 Research in Electron Theory of Crystals (Washington, DC: US Atomic Energy Commission)

[18] Feynman R P 1955 Phys. Rev. 97 660

[19] Feynman R P 1990 Statistical Mechanics: A Set of Lectures (Reading, MA: Addison-Wesley)

[20] Prokof'ev N V and Svistunov B V 1998 Phys. Rev. Lett. 81 2514

[21] Mishchenko A S, Prokof'ev N V, Sakamoto A and Svistunov B V 2000 Phys. Rev. B 62 6317

[22] Massignan P, Zaccanti M and Bruun G M 2014 Rep. Prog. Phys. 77 034401

[23] Vlietinck J, Ryckebusch J and van Houcke K 2013 Phys. Rev. B 87 115133

[24] Vlietinck J, Ryckebusch J and van Houcke K 2014 Phys. Rev. B 89 085119

[25] Forbes M M, Gezerlis A, Hebeler K, Lesinski T and Schwenk A 2014 Phys. Rev. C 89 041301
[26] Thornber K K 1974 *Phys. Rev.* B **9** 3489

[27] Pitaevskii L and Stringari S 2003 *Bose-Einstein Condensation* 1st edn (Oxford: Oxford University Press)

[28] Bruderer M, Bao W and Jaksch D 2008 *Europhys. Lett.* **82** 30004

[29] Blinova A A, Boshier M G and Timmermans E 2013 *Phys. Rev. A* **88** 053610

[30] Sumi A and Toyozawa Y 1973 *J. Phys. Soc. Japan* **35** 137

[31] Fantoni R 2012 *Phys. Rev.* B **86** 144304