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

TL;DR

This paper investigates phonon propagation and damping in liquid helium-4, providing detailed calculations of phonon mean-free path and lifetime across different pressures using advanced many-body theory, relevant for quantum experiments and dark matter detection.

Contribution

It offers the first detailed calculation of phonon mean-free path and lifetime in liquid helium-4 using dynamic many-body theory, highlighting the importance of advanced methods for accurate predictions.

Findings

Long wavelength phonons have large mean free paths.

Intermediate energy phonons can have mean free paths of micrometers.

Advanced many-body methods are essential for reliable predictions.

Abstract

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.

Equations

$$H=-\frac{{\hslash }^2}{2m}\sum _i{\nabla }_i^2+\sum _{i<j}v(|{\mathbf{r}}_i-{\mathbf{r}}_j|)$$

$${\Psi }_0({\mathbf{r}}_1,\dots ,{\mathbf{r}}_N)=\exp \frac{1}{2}\left[\sum _{i<j}{u}_2({\mathbf{r}}_i,{\mathbf{r}}_j)+\sum _{i<j<k}{u}_3({\mathbf{r}}_i,{\mathbf{r}}_j,{\mathbf{r}}_k)\right].$$

$$\frac{\delta}{\delta u_{i}({\bf r}_{1},\dots,{\bf r}_{i})}\frac{\bigl{\langle}{\Psi_{0}}\bigr{|}H\bigl{|}{\Psi_{0}}\bigr{\rangle}}{\bigl{\langle}{\Psi_{0}}\big{|}{\Psi_{0}}\bigr{\rangle}}=0\,.$$

$$|\Psi (t)⟩=e^{-\mathrm{i}{E}_{0}t/\hslash }\frac{e^{\frac{1}{2}\delta U(t)}|{\Psi }_0⟩}{[⟨{\Psi }_0|e^{\frac{1}{2}\delta U^{†}(t)}{e}^{\frac{1}{2}\delta U(t)}|{\Psi }_0⟩{]}^{1/2}},$$

$$\delta U(t)=\sum _i\delta u_1({\mathbf{r}}_i;t)+\sum _{i<j}\delta u_2({\mathbf{r}}_i,{\mathbf{r}}_j;t)+\dots$$

$$\frac{\delta }{\delta u_i({\mathbf{r}}_1,\dots {\mathbf{r}}_i;t)}\int dt⟨\Psi (t)|H-\mathrm{i}\hslash {\partial }_t|\Psi (t)⟩=0.$$

$$\chi (q,\omega )=\frac{S(q)}{\hslash \omega -{\epsilon }_{\mathrm{F}}(q)-\Sigma (q,\hslash \omega )}+\frac{S(q)}{-\hslash \omega -{\epsilon }_{\mathrm{F}}(q)-\Sigma (q,-\hslash \omega )},$$

$$\Sigma (q,\hslash \omega )=\frac{1}{2}\int \frac{d^3{k}_1{d}^3{k}_2}{(2\pi {)}^3\rho }\frac{\delta (\mathbf{q}-{\mathbf{k}}_1-{\mathbf{k}}_2){\left|{\tilde{V}}_3(\mathbf{q};{\mathbf{k}}_1,{\mathbf{k}}_2)\right|}^2}{\hslash \omega -{\epsilon }_{\mathrm{F}}(k_1)-\Sigma (k_1,\hslash \omega -{\epsilon }_{\mathrm{F}}(k_2))-{\epsilon }_{\mathrm{F}}(k_2)-\Sigma (k_2,\hslash \omega -{\epsilon }_{\mathrm{F}}(k_1))}.$$

$${\tilde{V}}_3(\mathbf{q};{\mathbf{k}}_1,{\mathbf{k}}_2)=\frac{{\hslash }^2}{2m}\sqrt{\frac{S(k_1)S(k_2)}{S(q)}}\left[\mathbf{q}\cdot {\mathbf{k}}_1\tilde{X}(k_1)+\mathbf{q}\cdot {\mathbf{k}}_2\tilde{X}(k_2)-q^2{\tilde{X}}_3(\mathbf{q},{\mathbf{k}}_1,{\mathbf{k}}_2)\right],$$

$${\epsilon }_0(k)={\epsilon }_{\mathrm{F}}(k)+\Sigma (k,{\epsilon }_0(k)).$$

$${\epsilon }_0(k)\approx \hslash ck(1-\gamma k^2)$$

$${\epsilon }_0(k)=\hslash ck\left(1-\gamma k^2+{\alpha }_3{k}^3+{\alpha }_4{k}^4\right),$$

$${\epsilon }_0(k)=\hslash ck\left(1-\gamma k^2\frac{1-k^2/Q_a^2}{1+k^2/Q_b^2}\right).$$

$${\alpha }_3=\frac{{\pi }^2}{24}\frac{\rho }{m{c}^2}{C}_6.$$

$${\epsilon }_{\mathrm{F}}(k)=\sqrt{t^2(k)+2t(k){\tilde{V}}_{\mathrm{p}-\mathrm{h}}(k)}$$

$${\tilde{V}}_{\mathrm{p}-\mathrm{h}}(k)\equiv \rho \int d^{3}r{V}_{\mathrm{p}-\mathrm{h}}(r)e^{\mathrm{i}\mathbf{k}\cdot \mathbf{r}}$$

$${\tilde{V}}_{\mathrm{p}-\mathrm{h}}(k)={\tilde{V}}_{\mathrm{p}-\mathrm{h}}(0+)+V_2{k}^2+\frac{{\pi }^2}{12}\rho C_6{k}^3$$

$$\Sigma (k,\hslash \omega )=-\frac{|{\tilde{V}}_3(\mathbf{k};-\mathbf{k}/2,-\mathbf{k}/2){|}^{2}k}{16\pi \rho {\epsilon }_0^{\prime }(k/2)}\sqrt{\frac{2{\epsilon }_0(k/2)-\hslash \omega }{{\epsilon }_0^{\prime \prime }(k/2)}}.$$

$$\tau (k)=\hslash \Im m[\Sigma (k,\hslash {\epsilon }_0(k),k){]}^{-1}$$

$$d(k)=v_g(k)\tau (k)$$

$$\mathbf{k}$$

$${\epsilon }_0(k)$$

$${\mathbf{j}}^{(2)}(\mathbf{r},t)$$

$$\frac{d}{dt}{j}^{(2)}(\mathbf{e})=\frac{1}{mS(k)}\int \frac{d^{3}q}{(2\pi {)}^3\rho }q\delta (\mathbf{e}-\frac{\mathbf{q}}{q})\Im m\sigma (\mathbf{k},\mathbf{q},\omega )$$

$$\sigma (\mathbf{k},\mathbf{q},\omega )=\frac{1}{2}\frac{{\left|\tilde{V}(\mathbf{k},\mathbf{k}-\mathbf{q},\mathbf{q})\right|}^2}{\hslash \omega -{\epsilon }_0(|\mathbf{k}-\mathbf{q}|)-{\epsilon }_0(q)+\mathrm{i}\eta }$$

$$U_{\mathrm{ext}}(\mathbf{r};t)=\left[U_{\mathrm{ext}}(\mathbf{r})e^{-\mathrm{i}\omega t}+U_{\mathrm{ext}}^{\ast }(\mathbf{r})e^{\mathrm{i}\omega t}\right]e^{\eta t}$$

$${\mathbf{j}}^{(2)}(\mathbf{r})=\frac{\hslash }{4m\mathrm{i}}\left[\delta {\rho }_1^{\ast }({\mathbf{r}}_1){\nabla }_1\delta u_1({\mathbf{r}}_1)+\int d^3{r}_2\delta {\rho }_2^{\ast }({\mathbf{r}}_1,{\mathbf{r}}_2){\nabla }_1\delta u_2({\mathbf{r}}_1,{\mathbf{r}}_2)\right].$$

$$\delta {\rho }_1({\mathbf{r}}_1)=\rho \left[S\ast \delta v_1\right]({\mathbf{r}}_1).$$

$$\left[S\ast f\right]({\mathbf{r}}_1)=f({\mathbf{r}}_1)+\rho \int d^3{r}_{2}h({\mathbf{r}}_1,{\mathbf{r}}_2)f({\mathbf{r}}_2).$$

$$\delta v_1(\mathbf{r})$$

$$\delta w_1(\mathbf{r})$$

Full text

††journal: Journal of Low Temperature Physics

1 Introduction

It has been known for a long time MarisRMP -Sridhar87 that low-energy phonons in liquid 4He display an anomalous dispersion relation which allows these phonons to decay. Precise neutron scattering measurements of the phonon dispersion in liquid 4He skwpress 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 $\rho=0.0245$Å*-3*. These data agree very well with our recent theoretical results for the dynamics of 4He based on time-dependent multiparticle correlations eomIII ; 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 4He droplets dropdyn and the surface of 4He hscatt . 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 4He 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 4He has been proposed as a detection mechanism for low-mass dark matter PRL117Zurek -Romani2019 , among other proposed detectors PhysRevLett.119.131802 . 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.

2 Ground State Structure of 4He

Microscopic calculation of properties of many–body systems begin with an accurate calculation of ground–state properties. For 4He it is adequate to begin with a non-relativistic Hamiltonian

[TABLE]

where the pair-wise interaction $v(r)$ is taken from Aziz et al. AzizII .

The most efficient evaluation of ground state properties is done by the variational Jastrow-Feenberg ansatz for the ground state:

[TABLE]

The correlation functions $u_{i}({\bf r}_{1},\dots,{\bf r}_{i})$ are obtained by minimizing the ground state energy $E_{0}$

[TABLE]

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.

3 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

[TABLE]

where \bigl{|}{\Psi_{0}}\bigr{\rangle} is the ground state, and $\delta 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:

[TABLE]

The amplitudes $\delta u_{i}({\bf r}_{1},\ldots{\bf r}_{i};t)$ are determined by the time-dependent generalization of the Ritz’ variational principle:

[TABLE]

Assuming that $\delta U(t)$ is a small perturbation of the ground state correlations, we can linearize the equations of motion for $\delta u_{i}({\bf r}_{i},\dots;t)$, leading to the density–density response function $\chi(k,\omega)$ from which we obtain the dynamic structure function $S(k,\omega)=\Im m\,\chi(k,\omega)$:

[TABLE]

where $\varepsilon_{\rm F}(k)=\hbar^{2}k^{2}/2mS(k)$ is the Feynman excitation spectrum, and the self-energy is given by an integral equation

[TABLE]

$\tilde{V}_{3}({\bf k};{\bf p},{\bf q})$ is the three-phonon vertex

[TABLE]

where $\tilde{X}(k)=1-1/S(k)$. $\tilde{X}_{3}({\bf q},{\bf k}_{1},{\bf k}_{2})$ is the fully irreducible three-phonon coupling matrix element. In the simplest approximation, $\tilde{X}_{3}({\bf q},{\bf k}_{1},{\bf k}_{2})$ is replaced by the three–body correlation $\tilde{u}_{3}({\bf q},{\bf k}_{1},{\bf k}_{2})$; this approximation ensures that long–wavelength properties of the excitation spectrum are preserved Chuckphonon . Improved calculations eomI sum a 3-point integral equation to ensure that exact properties of $\tilde{X}_{3}({\bf q},{\bf k}_{1},{\bf k}_{2})$ as $q\rightarrow 0$ and of the Fourier transform $X_{3}({\bf r}_{1},{\bf r}_{2},{\bf r}_{3})$ for $|{\bf r}_{1}-{\bf r}_{2}|\rightarrow 0$ and $|{\bf r}_{1}-{\bf r}_{3}|\rightarrow 0$ are satisfied eomI . 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 JacksonSkw 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 eomIII and experimental results skw4lett which are still being explored skwpress .

In Fig. 2, we compare the experimental dynamic structure function $S(k,\omega)$ (top panel) skw4lett with results from DMBT calculations (middle panel) and CBF-BW calculation JacksonSkw (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 part $\tilde{X}_{3}({\bf q};{\bf k}_{1},{\bf 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 $\varepsilon_{\rm 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.

4 Phonon dispersion

The theoretical phonon dispersion relation is obtained by solving the implicit equation

[TABLE]

At long wave lengths, the phonon dispersion relation is

[TABLE]

where $c$ is the speed of sound and $\gamma$ is the phonon dispersion coefficient. If $\gamma<0$, we speak of anomalous dispersion, long–wavelength phonons are damped.

To obtain the dispersion coefficient $\gamma$ we have fitted both experimental and theoretical data by a polynomial of the form

[TABLE]

from which the phonon dispersion coefficient was extracted. The polynomial form turned out to be more flexible than the Padé approximation MarisRMP ; PhysRevA.8.1980

[TABLE]

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 Pitaevskii1970 ; Feenberg1971

[TABLE]

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 $\epsilon_{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 $\varepsilon_{\rm F}(k)$ is derived from a Bogoliubov formula

[TABLE]

where $t(k)=\hbar^{2}k^{2}/2m$ is the kinetic energy, and

[TABLE]

is the “particle-hole” interaction or, in the language of Aldrich and Pines Aldrich the “pseudopotential”. For what follows we only need the property that, for large distances, $V_{\rm p-h}(r)$ falls off like the bare interaction, i.e. $V_{\rm p-h}(r)\sim C_{6}r^{-6}$ for $r\rightarrow\infty$. Then, $\tilde{V}_{\rm p-h}(k)$ has the expansion

[TABLE]

where $\tilde{V}_{\rm p-h}(0+)=\rho\int d^{3}rV(r)=mc^{2}$ is related to the speed of sound, and $V_{2}$ is the second moment $V_{2}=-\frac{\rho}{6}\int d^{3}rV_{\rm p-h}(r)r^{2}$. Higher moments do not exist due to the van der Waals tail.

Fig. 3 shows the full $\tilde{V}_{\rm 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\leq k\leq 0.1\,$Å*-1* which is inaccessible to neutron scattering measurements. The figure also shows the Feynman spectrum $\varepsilon_{\rm F}(k)$ as obtained from the full $\tilde{V}_{\rm p-h}(k)$ and from the moment expansion.

Both the DMBT and the CBF-BW results for $\gamma$ agree quite well with the experimental data, whereas the Feynman approximation predicts anomalous dispersion at all densities.

The calculated phase velocities $c=\varepsilon_{0}(k)/\hbar k$ are given in Fig. 4 for densities covering the full range between the saturated vapor pressure and solidification. Also shown are experimental data skw4lett ; skwpress 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 skw4lett ; skwpress .

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 $\gamma$ for several 4He densities. The speed of sound $c$ determined by the fit agrees well with the well know values measured by ultrasonic or thermodynamic techniques Abraham ; CaupinJordi , 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 $\gamma$. We thus show in Fig. 5 the curves for $\gamma$ 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 coefficients $\hbar c$, the speed of sound $c$, the Grüneisen-Constant $u$ and the dispersion coefficient $\gamma$. 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 $\gamma$ turns positive above $\rho\approx 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.

5 Phonon Mean Free Path

If $\gamma<0$, a phonon of energy/momentum $(\hbar\omega,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 $(\hbar\omega,k)$ is at the critical energy $\hbar\omega_{\rm crit}=2\varepsilon_{0}(k/2)$, where $\varepsilon_{0}(k)$ is the dispersion relation calculated with DMBT, see eq. (10). At a corresponding critical wave number $k_{\rm crit}$ determined by solving $2\varepsilon_{0}(k/2)=\varepsilon_{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 bosegas ; eomIII

[TABLE]

Phonons with lower wave numbers, $k<k_{\rm 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_{\rm 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_{\rm crit}$.

The life time of phonons with $k<k_{\rm crit}$ can be readily calculated from the imaginary part of the self-energy,

[TABLE]

The mean free path of a phonon, important for the suggested application of superfluid 4He for dark matter detection mentioned above, can be obtained from

[TABLE]

where $v_{g}={d\omega(k)\over dk}$ is the group velocity.

Evidently, two things are needed for a reliable theoretical prediction of the phonon life time $\tau(k)$ and of the mean free path $d(k)$. (1) An accurate dispersion coefficient is required for the low-momentum kinematics,

[TABLE]

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 $\rho$. The density range is smaller in the latter case, because the dispersion relation is not anomalous anymore for $\rho\gtrsim 0.0245$Å*-3*, see the dispersion coefficient $\gamma$ 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 $\gamma$ for all densities considered here (see Fig.5). Therefore CBF-BW yields a finite mean free path even for $\rho=0.0260$Å*-3*. In fact $d(k)$ is almost independent of $\rho$ 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 $\rho$, hence on the pressure. While CBF-BW predicts a minimal decay length of about $0.1\mu$m for all densities, the DMBT results show that the minimal decay length is $1-10\mu$m around the equilibrium density, and increases significantly with density and pressure, because $\gamma$ crosses zero and becomes positive at high density, Fig.5, until the phonons do not decay anymore for densities higher than $0.0240$Å*-3*,

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.

6 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 $\delta U({\bf r}_{1},\ldots,{\bf r}_{N},t)$

[TABLE]

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 4He surface may be found in Ref. 8, and for the case of impurity and 4He scattering in Ref. 7.

The decay of a phonon with wave vector ${\bf 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 ${\bf e}$ (${\bf e}$ is a unit vector with Euler angles $\theta$ and $\varphi$). For this purpose we calculate the rate of inelastic phonon current in direction ${\bf e}$, for which we obtain

[TABLE]

The $\delta({\bf e}-{{\bf q}/q})$ obviously selects only phonons with wave vector ${\bf q}$ parallel to the direction of interest. Here $\sigma({\bf k},{\bf q},\omega)$ is

[TABLE]

The integration over ${\bf q}$ yields the self energy $\int\frac{d^{3}q}{(2\pi)^{3}\rho}\sigma({\bf q};{\bf k},\omega)=\Sigma({\bf k},\omega)$, see Eq. (8). For the calculation of $\sigma({\bf k},{\bf q},\omega)$, 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. 7 shows, at the density $\rho=0.0215$Å*-3*, the DMBT prediction for the rate with which a phonon with wave number ${\bf k}$ decays into phonons with direction $\theta$ with respect to ${\bf k}$ (obviously, the rate does not depend on $\varphi$). 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.7 shows that these higher energy phonons decay into lower energy phonons which lie in a narrow cone (small $\theta$). 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 $\theta=20^{\circ}$.

7 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 4He 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$\mu$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 eomIII ; 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 $\chi(k,\omega)$, Eq. (7), it is convenient to apply an weak harmonic perturbation (e.g. induced by a neutron beam)

[TABLE]

As usual, the perturbation is switched on adiabatically with an infinitesimal rate $\eta>0$. After linearization, the positive and negative frequency terms leads to a corresponding response of the correlation operator $\delta U(t)$ split into (complex) positive and negative frequency terms; similarly for the response of one-body, two-body densities, $\delta\rho_{n}$.

Inserting the correlation operator (5) at the two-body level in the transport currents (23) gives (we omit the time dependence for clarity)

[TABLE]

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 $\Sigma(k,\hbar\omega)$, summarized in the text, we express $\delta u_{1}$ and $\delta\rho_{2}$ in terms of $\delta\rho_{1}$ and $\delta u_{2}$. It turns out advantageous to introduce $\delta v_{1}({\bf r};t)$ such that

[TABLE]

where we define the convolution product with the static structure function

[TABLE]

In linear order, we can then write

[TABLE]

where the $Y({\bf r}_{1};{\bf r}_{2},{\bf r}_{3})$ is the set of those contributions to the three-body distribution function that is non-nodal in the point ${\bf r}_{1}$; details can be found in Ref. eomIII .

In all further manipulations, we use the “convolution approximation” which is the simplest approximation for high-order distribution functions that maintains the exact long-wavelength properties. In convolution approximation, $Y({\bf r}_{1};{\bf r}_{2},{\bf r}_{3})$ is simply

[TABLE]

Using Eqs. (30) and (31) and expressing $\delta\rho_{2}({\bf r}_{1},{\bf r}_{2})$ in terms of $\delta\rho_{1}({\bf r})$ and the fluctuation of the pair distribution function $\delta g({\bf r}_{1},{\bf r}_{2})$,

[TABLE]

we obtain

[TABLE]

Finally we eliminate $\delta g({\bf r}_{1},{\bf r}_{2})$ using the convolution approximation:

[TABLE]

which yields the length expression for the current

[TABLE]

The first term (34) is the current in Feynman approximation which is the elastic channel. The second and third terms, (35) and (36), are correlations between the density fluctuation $\delta\rho_{1}$ and the pair correlation fluctuation $\delta u_{2}$. We will show below that only the fourth term (37) contains the decay rate of an elementary excitation (in the above derivation given in Feynman approximation).

We now need explicit expressions for the one- and two-body fluctuations. These come from solving the equations of motion, formally given by Eq. (6), in detail given in Refs. 8; 6. We can assume that the perturbation $U_{\rm ext}({\bf r})$ is a plane wave with wave vector ${\bf q}$ which induces a plane wave density fluctuation, propagating with phase velocity given by the dispersion relation, $\omega_{0}(q)=\varepsilon_{0}(q)/\hbar$:

[TABLE]

where $\lambda\ll 1$ is an arbitrary strength factor.

The two-body correlation fluctuations are then coupled to the density fluctuation according to the equations of motion. Their spatial Fourier transform is given by

[TABLE]

We finally use Eq. (38) for $\delta\rho_{1}({\bf r})$ and Eq. (39) for $\delta u_{2}({\bf r}_{1},{\bf r}_{2})$ in the current ${\bf j}^{(2)}({\bf r})$ above. The fourth term (37) which becomes

[TABLE]

In the limit of adiabatically switching on the perturbation, $\eta\to 0$, we need to calculate the rate with which the inelastic current is generated, i.e. the time derivative. Using

[TABLE]

we obtain the rate of the total inelastic current $\frac{d}{dt}{\bf j}^{(2)}$ apart from the arbitrary strength factor $\lambda^{2}$. Selecting only the inelastic phonon current (the decay products) in a specific direction ${\bf e}$ relative to the incoming phonon with wave vector ${\bf q}$, we obtain the rate Eq. (24). Note that the two mixed term contributions to ${\bf j}^{(2)}({\bf r})$, Eqs. (35) and (36), do not have a finite contribution to the rate in the adiabatic limit.