Transport and Phonon Damping in $^{\bf 4}$He

The dynamic structure function $S(k,\omega)$ informs about the dispersion and damping of excitations. We have recently (Phys. Rev. B {\bf 97}, 184520 (2018)) compared experimental results for $S(k,\omega)$ from high-precision neutron scattering experiment and theoretical results using the ``dynamic many-body theory'' (DMBT), showing excellent agreement over the whole experimentally accessible pressure regime. This paper focuses on the specific aspect of the propagation of low-energy phonons. We report calculations of the phonon mean-free path and phonon life time in liquid \he4 as a function of wave length and pressure. Historically, the question was of interest for experiments of quantum evaporation. More recently, there is interest in the potential use of $^4$He as a detector for low-energy dark matter (K. Schulz and Kathryn M. Zurek, Phys. Rev. Lett. {\bf 117}, 121302 (2016)). While the mean free path of long wave length phonons is large, phonons of intermediate energy can have a short mean free path of the order of $\mu$m. Comparison of different levels of theory indicate that reliable predictions of the phonon mean free path can be made only by using the most advanced many--body method available, namely, DMBT.


Introduction
It has been known for a long time [1]- [4] that low-energy phonons in liquid 4 He display an anomalous dispersion relation which allows these phonons to decay. Precise neutron scattering measurements of the phonon dispersion in liquid 4 He [5] provide the phonon dispersion relation for all experimentally accessible densites with unprecedented accuracy. They confirm the finding of earlier work that the phonon dispersion relation is anomalous up to densities of about ρ = 0.0245Å −3 . These data agree very well with our recent theoretical results for the dynamics of 4 He based on time-dependent multiparticle correlations [6]; our methods should therefore also be capable of quantitative microscopic predictions for the phonon lifetime. We apply our many-body theory of inelastic scattering, previously applied to scattering off 4 He droplets [7] and the surface of 4 He [8]. In the following sections, we briefly review both the theoretical methods used to calculate the relevant ground state properties (Section 2) and the dynamic features (Section 3) of the 4 He liquid, and the data analysis to obtain accurate values of the phonon dispersion coefficient.
The investigations are, among others, of interest because multiple scattering in superfluid 4 He has been proposed as a detection mechanism for low-mass dark matter [9]- [12], among other proposed detectors [13]. Such a detector design requires accurate knowledge of the propagation of low-energy phonons within the medium. In particular, when the phonon dispersion relation is anomalous, phonons are damped and have a finite mean free path.

Ground State Structure of 4 He
Microscopic calculation of properties of many-body systems begin with an accurate calculation of ground-state properties. For 4 He it is adequate to begin with a non-relativistic Hamiltonian where the pair-wise interaction v(r) is taken from Aziz et al. [14].
The most efficient evaluation of ground state properties is done by the variational Jastrow-Feenberg ansatz for the ground state: The correlation functions u i (r 1 , . . . , r i ) are obtained by minimizing the ground state energy E 0 δ δ u i (r 1 , . . ., r i ) The method is known as Jastrow-Feenberg-Euler-Lagrange (JF-EL) method. The key quantity obtained from such a ground state calculation and used for the calculation of the dynamics is the static structure function S(k). We show in Fig. 1 a comparison of different calculations and experiments.

Many-Body Dynamics
Dynamics is treated at the same level as the ground state: We write the dynamic wave function as containing a small, time-dependent component where Ψ 0 is the ground state, and δU(t) is an excitation operator that is written, for the case of bosons, in exactly the same manner as the ground state wave function: PIMC [15] HNC−EL [6] Svensson et al. [17] Robkoff and Hallock [18] Cowley and Woods [16] Gibbs [19] Fig. 1 (color online) Static structure function S(k) for 4 He at equilibrium density. We compare JF-EL, Monte Carlo [15], and experiments [16,17,18,19]. The figure is from Ref. 6.
The amplitudes δ u i (r 1 , . . . r i ;t) are determined by the time-dependent generalization of the Ritz' variational principle: Assuming that δU(t) is a small perturbation of the ground state correlations, we can linearize the equations of motion for δ u i (r i , . . .;t), leading to the density-density response function χ(k, ω) from which we obtain the dynamic structure function S(k, ω) = ℑm χ(k, ω): where ε F (k) =h 2 k 2 /2mS(k) is the Feynman excitation spectrum, and the self-energy is given by an integral equation whereX(k) = 1 − 1/S(k).X 3 (q, k 1 , k 2 ) is the fully irreducible three-phonon coupling matrix element. In the simplest approximation,X 3 (q, k 1 , k 2 ) is replaced by the three-body correlationũ 3 (q, k 1 , k 2 ); this approximation ensures that long-wavelength properties of the excitation spectrum are preserved [20]. Improved calculations [21] sum a 3-point integral equation to ensure that exact properties ofX 3 (q, k 1 , k 2 ) as q → 0 and of the Fourier transform X 3 (r 1 , r 2 , r 3 ) for |r 1 − r 2 | → 0 and |r 1 − r 3 | → 0 are satisfied [21]. We include these corrections routinely, they have a visible effect only for wave vectors between the maxon and the roton. The CBF-Brillouin-Wigner (CBF-BW) approximation [22] is obtained by omitting the self-energy corrections in the energy denominator of Eq. (8).
The implementation of the method outlined only briefly here has led to an unprecedented agreement between theoretical predictions [6] and experimental results [23] which are still being explored [5].  (8), the lowest figure shows the CBF-BW approximation; the only difference to previous work [22] are more accurate input functions S(k) and X 3 (r 1 ,r 2 ,r 3 ).
In Fig. 2, we compare the experimental dynamic structure function S(k, ω) (top panel) [23] with results from DMBT calculations (middle panel) and CBF-BW calculation [22] (bottom panel). In the CBF-BW calculation we have, of course, used the best available values for S(k) in the three-phonon vertex (9) and included the irreducible partX 3 (q; k 1 , k 2 ) as described in Ref. 21. Thus, the numerical values are different from those of Ref. 22 but the physics described is basically the same.
The most significant difference between the CBF-BW approximation and the full solution of the integral equation (8) is the void regime above the phonon-maxon-roton dispersion relation. This is caused by the fact that, in the CBF-BW approximation, excitations decay into Feynman phonons ε F (k) and not into the physical phonons. Furthermore, CBF-BW overestimates the excitation energies in the maxon region, while DMBT agrees with the experiment.

Phonon dispersion
The theoretical phonon dispersion relation is obtained by solving the implicit equation At long wave lengths, the phonon dispersion relation is where c is the speed of sound and γ is the phonon dispersion coefficient. If γ < 0, we speak of anomalous dispersion, long-wavelength phonons are damped.
To obtain the dispersion coefficient γ we have fitted both experimental and theoretical data by a polynomial of the form from which the phonon dispersion coefficient was extracted. The polynomial form turned out to be more flexible than the Padé approximation [1,2] In particular, the Padé approximation (13) does not contain the term proportional to k 3 which can be calculated analytically from the asymptotic form of the microscopic two-body interaction. Assuming the typical asymptotic form V (r) = C 6 r −6 , one arrives at [24,25] We should note that fitting procedure should not be understood as a rigorous low-k expansion in the sense of a Taylor expansion of the dispersion relation ε 0 (k) around k = 0, but rather as a fit to the data in the theoretically and experimentally relevant regime. The fit works well for k < 0.6Å −1 .
To explain the fact that the above fitting procedure should not be considered to be a rigorous Taylor expansion, we must review a little more of the theoretical background. The Feynman spectrum ε F (k) is derived from a Bogoliubov formula  (17) Feynman spectrum (15) Spectrum with moment expansion where t(k) =h 2 k 2 /2m is the kinetic energy, and is the "particle-hole" interaction or, in the language of Aldrich and Pines [26] the "pseudopotential". For what follows we only need the property that, for large distances, whereṼ p−h (0+) = ρ d 3 rV (r) = mc 2 is related to the speed of sound, and V 2 is the second Higher moments do not exist due to the van der Waals tail. Fig. 4 shows the fullṼ p−h (k) the way it is used in the Bogoliubov formula (15) and the expansion (17) as calculated from the potential moments. Evidently, the agreement is very good only in the regime 0 ≤ k ≤ 0.1Å −1 which is inaccessible to neutron scattering measurements. The figure also shows the Feynman spectrum ε F (k) as obtained from the full V p−h (k) and from the moment expansion.
Both the DMBT and the CBF-BW results for γ agree quite well with the experimental data, whereas the Feynman approximation predicts anomalous dispersion at all densities.
The calculated phase velocities c = ε 0 (k)/hk are given in Fig. 4 for densities covering the full range between the saturated vapor pressure and solidification. Also shown are experimental data [23,5] at four pressures corresponding to the same the same density range. Clearly the agreement is excellent. Only a very small shift in density is observed between the theoretical calculations and the experimental results, already discussed in previous publications [23,5].
The polynomial expression (12) provides an excellent fit of the theoretical curves in the small wave vector range. In practice, fits of the experimental curves were done in the range 0.18 < k < 0.6Å −1 , the lower bound being determined by the neutron detectors smallest angle. The highest bound was determined as the maximum wave vector where stable fits could be obtained using the form (12). This was also verified for the theoretical curves; for the latter, the fit range could be extended to k = 0, without affecting significantly the results of the fits, as seen in Fig. 5.
The experimental results have been analyzed using the polynomial expression (12) to determine the dispersion coefficient γ for several 4 He densities. The speed of sound c determined by the fit agrees well with the well know values measured by ultrasonic or thermodynamic techniques [27,28], as can be seen in the extrapolations to k = 0 in Fig. 5. The agreement is not perfect, however, and this affects the value of the next term in the expansion, the dispersion coefficient γ. We thus show in Fig. 5 the curves for γ determined using either the fitted values of c, or the ultrasonic ones.
For completeness, we have also extracted other data from both the calculations and the experiments. Table 1 gives the calculated coefficientshc, the speed of sound c, the Grüneisen-Constant u and the dispersion coefficient γ. Static quantities like c and u can also be obtained from Monte Carlo simulations, Evidently, the agreement between the DMBT results and the experiment is excellent at all densities. The dispersion coefficient γ turns positive above ρ ≈ 0.0245Å −3 , meaning that long-wavelength phonons can propagate freely only at high pressures. The Feynman approximation predicts, on the other hand, a negative dispersion coefficient at all densites.

Phonon Mean Free Path
If γ < 0, a phonon of energy/momentum (hω, k) can decay into two phonons of lower energy and longer wave length. As a consequence, the self-energy becomes complex on the phonon dispersion relation. The onset of the imaginary part of an excitation of energy and momentum (hω, k) is at the critical energyhω crit = 2ε 0 (k/2), where ε 0 (k) is the dispersion relation calculated with DMBT, see eq. (10). At a corresponding critical wave number k crit determined by solving 2ε 0 (k/2) = ε 0 (k) for k, a phonon can decay into two phonons, parallel to the original phonon but with half the wave number. At this point the imaginary part of the self-energy is [29,6] Phonons with lower wave numbers, k < k crit , have more decay channels because the wave vectors of the produced phonons need not be parallel to the original phonon; this angular spread is discussed further below. Phonons with wave numbers k > k crit do not decay into pairs of phonons and have infinite life-time in the DMBT approximation. However, higher-order processes, ı.e. decay into three and more phonons, lead to a long, but finite life-time also for k > k crit .
The life time of phonons with k < k crit can be readily calculated from the imaginary part of the self-energy, The mean free path of a phonon, important for the suggested application of superfluid 4 He for dark matter detection mentioned above, can be obtained from where v g = dω(k) dk is the group velocity. Evidently, two things are needed for a reliable theoretical prediction of the phonon life time τ(k) and of the mean free path d(k). (1) An accurate dispersion coefficient is required for the low-momentum kinematics, which in particular determines the range of momenta where phonon damping occurs. (2) The three-phonon vertex V 3 is required for an accurate damping strength. The comparison of the dynamic structure function between experiment and the DMBT result in Fig. 2 shows that DMBT is accurate enough that eqns. (19) and (20) indeed provide a reliable theoretical prediction of phonon life time and mean free path. Figs. 6 show our results for the phonon mean free path d(k) in both CBF-BW approximation (left panel) and in DMBT (right panel) for several helium densities ρ. The density range is smaller in the latter case, because the dispersion relation is not anomalous anymore for ρ 0.0245Å −3 , see the dispersion coefficient γ shown in Fig. 5. In the CBF-BW calculation phonons decay by producing Feynman-phonons, which overestimates the anomalous dispersion and which have a negative γ for all densities considered here (see Fig.5). Therefore CBF-BW yields a finite mean free path even for ρ = 0.0260Å −3 . In fact d(k) is almost independent of ρ and thus of the pressure in the CBF-BW approximation, while the much improved DMBT result shows that the range of momenta where phonons can decay strongly depends on ρ, hence on the pressure. While CBF-BW predicts a minimal decay length of about 0.1µm for all densities, the DMBT results show that the minimal decay length is 1 − 10µm around the equilibrium density, and increases significantly with density and pressure, because γ crosses zero and becomes positive at high density, Fig.5, until the phonons do not decay anymore for densities higher than 0.0240Å −3 ,  Fig. 6 The phonon mean free path d(k) as function of the wave number for different pressures. The left panel shows the result in CBF-BW approximation and the right one shows the improved result obtained with the DMBT method which shows that the phonon mean free path depends strongly on the pressure.
While the DMBT calculation predicts a much smaller range of wave numbers at which phonons are dampened, the values for mean free path d(k) at different densities are almost identical for a given k, if there is damping at all. In fact, the d(k) values predicted by CBF-BW and DMBT are essentially the same. The reason is that the lifetime and thus the mean free path are mostly determined by the vertex V 3 , which is the same for CBF-BW and DMBT. The strength of damping does depend on the first derivative of the dispersion which, however, for the small k range relevant here is very similar in all approximations (Feynman, CBF-BW, and DMBT) -only the higher order derivatives, essential for the correct decay kinematics, are sensitive to the approximation.

Inelastic currents
Anomalous dispersion is a prerequisite for phonon decay, but the deviation of the dispersion relation from a linear dispersion is, in the wave number regime under consideration here, hardly visible, see Fig. 4. Therefore, a phonon decays into a pair of more or less collinear phonons. The angular distribution of these decay phonons can be obtained from the transport current, which is the second order expectation value of the current operator in the fluctuations δU(r 1 , . . . , r N ,t) The derivation of workable formulas is somewhat tedious, some essential steps will be presented in the Appendix. A detailed derivation for the case of impurity scattering off the 4 He surface may be found in Ref. 8, and for the case of impurity and 4 He scattering in Ref. 7. The decay of a phonon with wave vector k produces two phonons according to the kinematics (21) and (22). We are interested in a measure for the probability that a decay product is ejected in direction e (e is a unit vector with Euler angles θ and ϕ). For this purpose we calculate the rate of inelastic phonon current in direction e, for which we obtain The δ (e − q/q) obviously selects only phonons with wave vector q parallel to the direction of interest. Here σ (k, q, ω) is The integration over q yields the self energy d 3 q (2π) 3 ρ σ (q; k, ω) = Σ (k, ω), see Eq. (8). For the calculation of σ (k, q, ω), we have replaced the non-local self-energy appearing in Eq. (8) by the dispersion relation (10), this is legitimate in the low momentum regime under consideration. Fig. 6 shows, at the density ρ = 0.0215Å −3 , the DMBT prediction for the rate with which a phonon with wave number k decays into phonons with direction θ with respect to k (obviously, the rate does not depend on ϕ). The figure summarizes both the kinematics and the strength of phonon damping by decay into lower momentum phonons. In agreement with the decay length in the right panel of Fig. 6, the decay rate is highest for phonons with momentum k close to 0.5Å −1 . Fig.6 shows that these higher energy phonons decay into lower energy phonons which lie in a narrow cone (small θ ). Phonons with k = 0.2 − 0.3Å −1 decay at a smaller rate (have a larger decay length, see Fig.6), but generate phonons in a larger cone up to angles of θ = 20 • .

Conclusion
The dispersion and damping of low momemtum phonons in superfluid helium-4 is of practical interest for cosmological applications in weakly interacting particle detectors to test hypotheses for dark matter. We presented a comparison of the low momentum dispersion in 4 He between experimental data and the results from dynamic many-body theory (DMBT) and find excellent agreement. The dispersion obtained with DMBT is a significant improvement over the CBF prediction in the full momentum and energy range of interest. This allows us to make quantitative predictions for the decay of phonons due to the anomalous dispersion without phenomenological or experimental input parameters. The phonon mean free path is similar to CBF-BW results, but the CBF-BW approximation fails to predict the pressure dependence, leading to decay lengths as short as 0.1µm. This shortcoming is remedied by the DMBT which shows that the minimal decay length is strongly pressure dependent.

Transport current
We present in this appendix the essential steps of the derivation of the transport currents; the details of the derivation are rather similar to the ones for the impurity currents derived in Refs. 8 and 7. We derive the inelastic current at the level of pair fluctuations, i.e. we truncate the expansion (5) at the two-body level. To derive the integral equation (8) one needs to include fluctuations to all orders [6]; our result will be plausible enough such that we can avoid these complications.
As in the derivation of the linearized equations of motion for the correlation operator (5), leading to the density-density response function χ(k,ω), Eq. (7), it is convenient to apply an weak harmonic perturbation (e.g. induced by a neutron beam) As usual, the perturbation is switched on adiabatically with an infinitesimal rate η > 0. After linearization, the positive and negative frequency terms leads to a corresponding response of the correlation operator δU(t) split into (complex) positive and negative frequency terms; similarly for the response of one-body, two-body densities, δ ρ n . Inserting the correlation operator (5) at the two-body level in the transport currents (23) gives (we omit the time dependence for clarity) j (2) (r) =h 4mi δ ρ * 1 (r 1 )∇ 1 δ u 1 (r 1 ) + d 3 r 2 δ ρ * 2 (r 1 ,r 2 )∇ 1 δ u 2 (r 1 ,r 2 ) .
where we have used that we are only interested in the time-averaged currents, since only those are observable by a detector. While looking very simple, the expression (27) for the current contains both the fluctuations of the correlation and of the densities. In accordance with the derivation of the self energy Σ (k,hω), summarized in the text, we express δ u 1 and δ ρ 2 in terms of δ ρ 1 and δ u 2 . It turns out advantageous to introduce δ v 1 (r;t) such that δ ρ 1 (r 1 ) = ρ [S * δ v 1 ] (r 1 ) .
In the limit of adiabatically switching on the perturbation, η → 0, we need to calculate the rate with which the inelastic current is generated, i.e. the time derivative. Using d dt e 2ηt (ε 0 (q) − ε 0 (k 1 ) − ε 0 (k 2 )) 2 +h 2 η 2 → 2π h δ (ε 0 (q) − ε 0 (k 1 ) − ε 0 (k 2 )) (40) we obtain the rate of the total inelastic current d dt j (2) apart from the arbitrary strength factor λ 2 . Selecting only the inelastic phonon current (the decay products) in a specific direction e relative to the incoming phonon with wave vector q, we obtain the rate Eq. (24). Note that the two mixed term contributions to j (2) (r), Eqs. (35) and (36), do not have a finite contribution to the rate in the adiabatic limit.