{'fraction_non_alphanumeric': 0.04454439930354034, 'fraction_numerical': 0.015851712130005805, 'mean_word_length': 4.144429930957268, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Mean-field vs RPA calculation of the energy of an impurity immersed in a spin 1/2 superfluid

In this article we calculate the energy of an impurity weakly coupled to a spin 1/2 fermionic superfluid. We show that the divergences resulting from three-body physics can only be cured using a proper description of the excitations of the many-body background. We highlight the crucial role played by interactions between quasiparticles which are overlooked within BCS (Bardeen-Cooper-Schrieffer) mean-field theory of fermionic superfluidity. By contrast, we prove that their addition using the Random Phase Approximation (RPA) allows us to regularize the energy of the impurity. Finally, we show that these beyond mean-field corrections should be observable by the analysis of the frequency shift of the impurity center of mass oscillations in an external confining potential.

In this article we calculate the energy of an impurity weakly coupled to a spin 1/2 fermionic superfluid. We show that the divergences resulting from three-body physics can only be cured using a proper description of the excitations of the many-body background. We highlight the crucial role played by interactions between quasiparticles which are overlooked within BCS (Bardeen-Cooper-Schrieffer) mean-field theory of fermionic superfluidity. By contrast, we prove that their addition using the Random Phase Approximation (RPA) allows us to regularize the energy of the impurity. Finally, we show that these beyond mean-field corrections should be observable by the analysis of the frequency shift of the impurity center of mass oscillations in an external confining potential.

I. INTRODUCTION AND MAIN RESULTS
The concept of quasi-particle provides a simple and powerful framework for the study of the low-energy physics of complex many-body systems. In practice, a quasi-particle can be described as a free particle dressed by a cloud of excitations of the surrounding many-body medium leading to a renormalization of its physical properties, and most notably its mass.
This mechanism was first suggested by Landau and Pekar to describe the coupling of electrons to the vibration modes of a crystal lattice (the so-called polaron problem) [1] and since then has been generalized to a host of physical situations, from solid state physics where polarons play a crucial role in the study a photovoltaic materials [2], to nuclear physics [3].
In the recent years, mixtures of ultracold atoms have provided an ideal playground for the study of polaron physics [4]. In this context, two simple systems were studied. Firstly, an analogous of Landau-Pekar polaron was obtained by immersing an impurity in a weakly interacting Bose-Einstein condensate (BEC) -the so-called Bose-polaron [5][6][7]. In this case, the bogoliubov spectrum describing the low-lying excitations of the BEC possess a structure that is similar to the phonons of a crystal. Another simple situation was obtained by considering the case of a particle swimming in a sea of non-interacting spin-polarized fermions, the so-called Fermi polaron [8][9][10][11][12][13][14] that was also realized in exciton-polariton systems [15].
A generalization of these two polaronic model systems is provided by the study of an impurity immersed in an ensemble of attractive spin 1/2 fermions [16][17][18]. When the attraction between the particles of the medium is varied, the ground state of the many-body background evolves from an ideal gas of fermions to a Bose-Einstein condensate of strongly bound dimers, thus realizing the celebrated BEC-BCS crossover [19]. As a consequence, the state describing the impurity immersed in a fermionic superfluid interpolates between the Fermi and Bosepolarons. While in reported experiments, the impurity is weakly coupled to the medium and most can be captured quantitatively by treating the impurity-fermion interaction within a mean-field approximation, a theoretical study of beyond-mean field effects was initiated in [20][21][22] that highlighted the role of three-body interactions. In particular, [22] showed that the leading order corrections were related to the compressibility of the background medium after a regularisation of UV divergences made possible by the introduction of explicit 3body interactions [23].
Building on this work, we show here that a precise description of the properties of the background is required for the regularization procedure used in [22] to be effective. Indeed, describing the background superfluid within BCS approximation that considers only pairbreaking excitations is incompatible with this renormalization scheme. To obtain a finite beyond-mean field correction, we need to take into account collective modes of the system, and for this we work within the framework of the Random Phase Approximation (RPA) [24][25][26] that allows us to make quantitative predictions for the energy of the impurity.
More precisely, consider an impurity of mass m i immersed in a many-body ensemble of spin 1/2 fermions of mass m f . a and a are respectively the fermion-fermion and impurity-fermion scattering lengths (the latter being assumed to be spin independent). We note |α the eigenstates of the medium in the absence of impurity, and E α the corresponding eigenenergies. By convention, α = 0 corresponds to the ground state. For a total density n of fermions, the state of the many-body background is characterized by the dimensionless parameter 1/k F a, where k F = (3π 2 n) 1/3 is the Fermi wavector, the limits 1/k F a → −∞ (resp. +∞) corresponding respectively to the weakly (resp. strongly) attractive regimes.
In the quasi-particle picture, the low-lying energy of an impurity of momentum q takes the form of a free arXiv:2202.03222v1 [cond-mat.quant-gas] 7 Feb 2022 particle dispersion relation where ∆E is the interaction energy of the impurity with the many-body ensemble and m * is its effective mass. Using perturbation theory, we have at second order in impurity-fermion coupling [22] 1 Here, g = 2π 2 a /m r is the coupling constant of the impurity-medium two-body contact interaction, m r is the reduced mass of an impurity-fermion pair, ε (α=i,r) q = 2 q 2 /2m α and wheren q = k,σ c † k,σ c k+q,σ is the Fourier transform of the density operator and N is the total number of fermions. χ p plays a central role in the following and we note that it depends only on the properties of the excitation spectrum of the many-body background and can be related to the density-density response function (see Sec. II).
n q is an operator that transfers a momentum − q to the many-body system. For large q, we can assume that it couples the ground-state to free particle excitations for which E α − E 0 2 q 2 /2m f . We therefore have As a consequence, we readily conclude that the sum appearing in the expression of the effective mass (Eq. (2)) is convergent and does not need any regularization. By contrast, the case of the interaction energy is more involved. Indeed, although χ 1 and 1/ε (r) compensate at leading order in Eq. (3), it was shown in [22] that the sum is log-divergent as a consequence of the following large momentum behaviour of χ 1 where C 2 is Tan's contact parameter of the many-body background [27], η = m i /m f and Following [23], this divergence can be regularized using an effective three-body interaction leading to the following regularized expression for the polaron energy shift where R 3 is an effective length characterizing the threebody scattering amplitude (R 3 /a 1.5 for the Lithium 6/Lithium 7 mixture used in [16]) and To calculate χ p , a first approach is to use the extension of the BCS theory-a mean-field theory-to the whole BEC-BCS crossover [22]. In this case, we readily obtain where (u k , v k ) are the amplitudes of the Bogoliubov modes and E k their energy [28]. However, this type of approximation scheme is restricted to the pair-breaking excitations sector and does not account for collective modes (phonons) that are dominant on the BEC side of the crossover. A well known consequence of this peculiarity is for instance that the compressibility calculated using BCS theory does not satisfy the f-sum rule. In our case, this leads to an underestimated value of κ (the third term appearing in Eq. (7) is missing, see sec. III) which leads to a UV-divergence of F when Λ → ∞ and prevents the regularization of the energy of the impurity.
Here, we solve this issue using the Random Phase Approximation (RPA) that incorporates quasi-particles interactions and recovers the collectives modes of the system (knonwn as Bogoliubov-Anderson's modes in this context) [24,25,29]. In this framework, we will show that χ 1 satisfies Eq. (6), with C 2 being given by BCS mean-field value C 2 /N = m 2 f ∆ 2 /n (here ∆ is the superconducting gap, see Sec. III). This allows us to regularize the polaron energy, leading to a universal function F displayed in Fig. 1. The numerical values obtained using RPA coincide with the asymptotic behaviours predicted in [22] for the BEC and BCS limits (1/k F a → ±∞, orange solid lines): On the weakly attractive side, F converges towards the finite limit F (−∞) = 3/2π that can be calculated analytically for an ideal Fermi gas, while on the far BEC limit, we recover the asymptotic expansion imposed by the matching between the polaron energy with the mean-field interaction of an impurity with a Bose-Einstein condensate of deeply-bound dimers: where C ad comes from the analysis of the three-body scattering problem and depends only on the mass ratio η. Note that the resummation scheme used in [20] also yields the correct velue for κ and is therefore compatible with the regularization procedure described previously. However this work incorporates quasi-particle interactions using only a scalar response function while the RPA considers a more general 3 × 3 response matrix that describes the coupling to both density and order parameter -see section II below. Likewise, RPA can be used to calculate the secondorder correction to the effective mass. Using dimensional analysis, we have The plot of G 1 is displayed in Fig.1, both within RPA and mean-field approximations. Both approaches coincide in the BCS limit where the excitation spectrum is dominated by pair-breaking excitations and where they asymptotically connect with the ideal gas prediction (orange line). In the BEC limit, the impurity is mostly dressed by a cloud of phonons. The BCS predictions therefore strongly underestimate the effective mass of the quasi-particle and RPA has to be used to obtain a quantitatively correct result. We see that in this regime the RPA prediction coincides with a calculation where the superfluid is described by a weakly repulsive Bose-Einstein condensate of dimers whose low-lying excitation spectrum is dominated by Bogoliubov excitations (orange line, see also Sec. IV).

II. RESPONSE FUNCTIONS IN THE MEAN-FIELD AND IN THE RANDOM PHASE APPROXIMATIONS
We now turn to a more detailed derivation of the results sketched in the previous secction. χ p is the central quantity giving access to both the energy and the effective mass of the polaron. To calculate it, we first note that it can be related to the density-density response of the system χ R . More precisely, let's consider the response of the superfluid to a perturbating potential V given by In the linear regime, the density response can be written as n q = Ω uχ R e −iωt , where, using standard perturbation theory, and Ω is a quantization volume ensuring that χ R is an intensive quantity. We note that χ 1 is proportional to the second term of χ R , and is associated with poles located on the negative part of the real axis. From this remark, we deduce that for E > 0 The expression of χ p for arbitrary p is then obtained by diferentiation with respect to E that yields We now turn to actual calculation of χ R of an ensemble of fermions with zero range interactions characterized by a scattering length a. The Hamiltonian describing the system is given by Here µ is the chemical potential and g 0 is the bare coupling constant that we can relate to g = 4π 2 a/m f using Lippman-Schwinger's equation A first approach to calculate the response of the system is to consider a mean-field approximation where we replace H by a quadratic Hamiltonian The optimal values for the Hartree potential U and the order parameter ∆ are obtained variationally [28] and we have U = g 0 ψ † ↑ ψ ↑ (r) = g 0 ψ † ↓ ψ ↓ (r) (we consider here a spin-balanced system) and ∆(r) = g 0 ψ ↓ (r)ψ ↑ (r) . The spectrum of the system can then be calculated using a standard self-consistent Bogoliubov transformation of the field operators and we note ∆ 0 and U 0 the ground state values of ∆ and U .
We can now calculate χ R by considering the response of the fermionic superfluid to the following perturbation: (19) where u and φ are driving fields coupling respectively to the density and the order parameter of the system. We first consider the response of the fermionic superfluid by describing its dynamics using the Hamiltonian H 0 [U 0 , ∆ 0 ] + H drive , where we assume that the expressions of the Hartree potential and of the order parameter appearing in H are not affected by the perturbation. As a consequence, this first approach restricts the response of the system to the pair-breaking sector.
In the framework of linear response theory, the response of the system in terms of order parameter and density is related to the drive fields u, φ + = (φ) and φ − = (φ) by a 3-by-3 correlation matrix M (q, ω) where the components M ij can be calculated explicitly as a function of the Bogoliubov amplitudes (u k , v k ) [26] [45]. In particular, we note that M 33 corresponds to the density-density response of the superfluid. We have indeed χ R MF = 2M 33 and applying Eq. (14) to χ R MF yields the mean field Eq. (10) for χ p .
The perturbation of the order parameter and of the Hartree potential contradicts the assumption that U and ∆ are not modified in the expression of H 0 . In the Random Phase Approximation, we solve this contradiction by considering a self-consistent response of the system to H drive . In other words, we now describes the system by We note that, by construction, , from which we conclude that the response of the system can be described using the previous calculation, but considering now an effective self-consistent drive defined by Using the non self-consistent approach, we have hence, using the definition of the effective driving fields, with The compressibility of the system corresponds to χ R RPA = 2 (M RPA ) 33 . In the zero-range limit where g 0 → 0, the diagonal terms M ii with i = 1, 2 need to be renormalized to obtain finite values. Introducing M ii = M ii − 1/g 0 , the compressibility can be written as a sum of a mean-field term and an interaction-induced term proportional to ∆ 2 [25]: with We note that, all matrix elements being calculated in the BCS superfluid, their imaginary part is non-zero only if the frequency ω is greater than 2∆ ; however, the zeros of the denominator of the beyond-mean-field term yield a non-zero imaginary part for frequencies below 2∆. In other terms, in addition to the contribution of the pair-breaking continuum, the RPA gives us both the dispersion relation of the collective modes and their contribution to the density-density response function. These collective modes were shown to turn into the usual Bogoliubov condensate modes in the BEC limit [25], and they become negligible in the BCS limit due to the exponential decrease of the gap.
We calculate numerically χ R and χ p using the approach laid out in [26], At large q, both the mean-field and RPA results behave as 1/q 2 and we display in Fig. 2 the sub-leading contribution for a unitary Fermi gas. As discussed in the introduction, we observe that the meanfield prediction underestimates the value of χ 1 while the RPA approach yields the asymptotic behavior given in Eq. (6) and is therefore compatible with a renormalization using three-body collisions.
III. HIGH-q BEHAVIOR OF χ1(q, q 2 /(2mi)) Here, we prove analytically that, as revealed numerically in previous section, χ 1 follows the asymptotic behavior (6). In section III A, we first prove that the dominant O(q −2 ) term comes from χ R MF , the mean-field contribution to the density-density response function and  (7) ). We see that the numerics agree with the predicted asymptotic behaviours. We see that, contrary to the MF result, the RPA result satisfies the asymptotic behaviour required by renormalization of the polaron energy.
that the mean-field contribution gives a O(q −3 ) term in (6), with a constant κ MF given by the first two terms in (7). Finally, in section III B, we prove that χ R int , the non mean-field contribution to the density-density response function gives a O(q −3 ) contribution to χ 1 (q, q 2 /(2m i )) and the third term in equation (7) for κ.

A. Subdominant O(q −3 ) contribution: MF term
As mentionned in the first section, the 1/q 2 behaviour of χ 1 originates from the high energy response of the system that corresponds to free particle excitations. As a consequence, the leading UV behaviour of χ is fully contained in the BCS mean-field term χ MF (see Eq. (5)(6). In order to reveal the sub-leading q −3 dependence, we subtract the term of order q −2 to χ MF, 1 and we have where we have used the notations E ± = E k± with k ± = k ± q/2 and ξ ± = ξ k± and E k = ξ 2 k + ∆ 2 , ξ k = k 2 /(2m f ) − µ. To obtain Eq.(27), we have furthermore used that n/2 = d 3 k (2π) 3 In the q → ∞ limit, we define a cutoff Λ such that where k low is a low energy wavevector scale defined by k low = max √ 2m f ∆, 2m f |µ| . In Appendix A, we show that that the contribution from the domains {k ± < Λ} (27) are negligible at leading order. In the domain {k ± > Λ}, we can replace at lowest order ( and E ± by k 2 ± /(2m f ) in the denominators. The only remaining scale is q and we rescale k by q. It is easily seen that the integral converges in +∞ and near k ± = 0, due to the first term in the integrand of (27) [46]. The contribution of order q −3 is where 29) J MF can be evaluated analytically and we find If we add the contributions of order q −2 and q −3 , we find the following high-q expansion for χ MF, 1 (q, q 2 2mi ) where In this expression, m 2 f ∆ 2 is Tan's contact per unit volume C 2 /Ω in a mean-field BCS theory. The ratio m 2 f ∆ 2 n is therefore equal to C 2 /N , where N is the total (summed on spins) number of particles. We recover a result similar to the Eq.(6), but with κ MF instead of the full κ. The non-MF part will complete the value of κ, as we will see in the next section.

B. Subdominant O(q −3 ) contribution: non MF term
We find the high-q behaviour of the non-mean field contribution of the response function χ R int (q, ω) in (26) by taking the ratio ω = ω/(q 2 /(2m f )) fixed. This means that in the q → ∞ limit, the frequency ω also tends to infinity.
We find that in this regime, These scaling behaviors are easily found using the expressions of the Π ij 's of [26] and the correspondance between Π ij and M ij .
Let us consider for instance M 23 . From equations (36) and (37) of [26], we have We separate the k space into two parts: k < Λ and k > Λ, where the cutoff Λ fulfills k low Λ q. For k < Λ, at lowest order, the integrand is of order ∆/(q 2 /m f ) 2 . Indeed, we can neglect k compared to q and ∆ and µ compared to q 2 /m f . The integration gives a volume factor of the order Λ 3 and the contribution to M 23 is of order ∆Λ 3 /(q 2 /m f ) 2 . Compared to the q −1 contribution, of the k > Λ domain this is negligible. Indeed, we find a ratio of order Λ 3 /q 3 1. In the k > Λ domain, we rescale k by q and we can neglect ∆ and µ in ξ ± and E ± . The integration volume gives a factor q 3 and the integrand a factor ∆/(q 2 /(2m f )) 2 . This gives a q −3 dependance.
In the same manner, we find If we use the asymptotic expressions of the M ij 's in (26), we find χ R int ∼ . Finally, we use (14) to determine the large q behavior of χ 1 (q, q 2 /(2m i )). where 2 ) . In the last step, we use the integral expression for X + and write X 2 + as a double integral on wave vectors and exchange the order of integrations (we integrate on ω first). After integration on the frequency ω, the integrals on the wave vectors factorize. The integral we need to cal- It is found to be equal to 1 2π arctan η η+2 . As a result we find for the high q behavior of the non mean field part of χ 1 As mentionned before, (m f ∆) 2 is Tan's contact per unit volume, and we see from (6) that the non mean-field contribution gives the third contribution to κ in (7).

IV. BCS AND BEC LIMITS FOR THE EFFECTIVE MASS
The expressions of χ R found using RPA and the meanfield approximations can also be used to obtain the effective mass of the polaron. As explained in the end of Sec.I, we can write the second-order perturbation to the effective mass m * in units of the impurity mass m i : it takes the form of a mean-field energy in units of the Fermi energy g n/E F squared times a dimensionless function G η (1/k F a) (Eq. (12)). Along the crossover, Eq.(2) gives us a relation between the response function χ 3 (q, ε (i) q ) and the correction to the effective mass G η : As mentionned in the introductory part, we see in Fig.  1 that the mean-field and RPA results are consistent in the BCS limit (that is, when 1/k F a tends to −∞). In this regime, the bath is indeed made of two non-interacting Fermi seas of opposite spins and its elementary excitations are particle-hole pairs which are correctly captured within the BCS approximation, as well as its RPA extension. χ 3 is then given by In particular, in the equal-mass case, an analytical calculation yields G 1 (−∞) = 3/8 that coincides with the numerical results displayed in Fig. (1).
In the opposite limit, 1/(k F a) → +∞, the mean-field approximation strongly underestimates the correction to the effective mass. Indeed, in this regime the fermions become so tightly bound that Cooper pairs lose their internal degrees of freedom, and we can treat the system as a condensate of composite bosons [30,31]. Pair breaking excitations are therefore suppressd and We can apply Eq.(2) to an ensemble of bosons, provided we replace χ 3 by the response function χ b 3 of the condensate, calculated within the mean-field Bogoliubov approximation. The fermionic physical parameters n, m f , a and µ are replaced by their bosonic equivalents: the boson density n b = n/2, the boson mass m b = 2m f and the boson-boson scattering length (in the framework of BCS theory) a bb = 2a which yields us the bosonic chemical potential µ b = 4πa bb m b n b = gn/2. We also need to replace the impurity-fermion coupling constant g by its impurity-boson counterpart given by Born's approximation g b = 2g . Eq.(2) transforms into with u b q and v b q being the amplitudes of the Bogoliubov modes of the condensate and E b q their energy. In this integral, energies and wave-vectors can be rescaled by the chemical potential and the healing length of the bosonic dimers. We then readily find that G 1 (1/k F a) scales like the square root of 1/k F a with where with I(1) 0.265. As shown in Fig. (1), this asymptotic behaviour coincides with the numerical results based on the RPA approximation.

V. CONCLUSION AND OUTLOOK
A possible way to test this prediction for the energy shift ∆E is to use radio-frequency spectroscopy which provided a very sensitive probe for the study of Fermi [12,13,32,33] and Bose [5,6,34] polarons, However, as described below, signatures of beyond mean-field effects might be already observable in the measurement of the oscillation frequency of an impurity of 7 Li inside a superfluid of 6 Li reported in [16].
To model this experiment, we assume that the impurity and the background superfuid are trapped by the same potential V . We can then describe the semiclassical dynamics of the polaron using the Hamiltonian h(r, p) = p 2 2m * + V (r) + ∆E(n(r)).
Within LDA, the density profile of the fermionic background can be obtained from its equation of state n(µ) using the prescription µ(r) = µ 0 − V (r), where µ 0 is the global chemical potential of the fermions.
For vanishingly small amplitude oscillations of the impurity, we have V (r) µ 0 and the energy of the quasiparticle can be approximated by h(r, p) = ∆E(n 0 )+ p 2 2m * + 1 + d∆E dn dn dµ 0 V (r), (40) where the subscript 0 means that the quantity is evaluated at the center of the fermionic cloud. If V is a harmonic potential characterized by a frequency ω, the effective oscillation frequency of the impurity is now If we treat the impurity-bath interaction in a meanfield approximation (in other terms, if we restrain ourselves to the first-order term in the perturbative expansions in g ), the relative frequency correction (ω * (1) − ω)/ω is given by g Consequently, if the beyond-mean-field relative frequency correction β is small, we can write it as the sum of the contributions of the second-order term in g in the ground-state energy and in the effective mass respectively: At unitarity, we get and β (2) We will here compute the value of these second-order frequency correction in the conditions reported in [16], namely k F a = 10 −2 , R 3 = 1.50a and m i = 7 6 m f . Knowing that at unitarity the equation of state of the superfluid is the same as the one of an ideal Fermi gas up to a multiplying constant ξ = 0.376 [19,[35][36][37][38][39][40], we obtain that the second term in the above expression (which we may call the three-body term) gives a relative correction about −4.7%, whereas the two-body term (the one involving F ) can be estimated to be around 1.1%.
Along the crossover, the correction due to the effective mass β which would give us a value of about −0.02%.
All these terms being second-order in g , their relative values at unitarity are fixed: the effective mass will always (in the RPA and as long as g is small) have an effect on the frequency correction that is two orders of magnitude smaller than the two-body term of the ground-state energy and about 250 times smaller than the three-body term.
Interestingly, [16] reported a ∼ 14 % upshift of the mean-field prediction with respect to the measured oscillation frequency at unitarity. This discrepancy is slightly beyond the 10% experimental error bar and the beyondmean-field correction calculated in the present paper lowers this difference to only ∼ 10 % that is now within the experimental uncertainty.
In this estimate of the beyond-meand field contribution to the frequency shift, we have neglected the drag between the two isotopes due to Andreev-Bashkin's effect [41] that was predicted first in nuclear physics and has recently been the focus of several theoterical studies in the context of ultracold gases [42][43][44]. The study of the interplay between these two effects is beyond the scope of this article and will be addressed in future work.
We now briefly justify why the integrations in domains {k − < Λ} and {k + < Λ} in (27) are negligible at this order. In the domain {k − < Λ}, we take the limit q → +∞ and keep k − finite. We easily find E+−ξ+ After angular integration, the first term of order q −3 vanishes. The term of order q −4 behaves as k −2 − if k − k low . After integration on k − (k − < Λ), it gives a contribution of order Λ/q 4 . Compared to the term of order q −3 , we find a ratio Λ/q. Therefore it is negligible. In the same way, we find the integration in domain 2 is negligible.