text:
# Gapless kinetic theory beyond the Popov approximation

## I Introduction

Finite-temperature theories of Bose-Einstein condensation (BEC) have been a very active field of study. The goal is a unified description of a dilute, atomic gas of bosons in a harmonic trap in terms of a condensate mean field interacting with a thermal cloud. The success of the zero-temperature Gross-Pitaevskii (GP) theory in describing BEC experiments spurred interest in effects not contained in this framework. Examples are damping of collective excitations and condensate growth through collisional redistributions of thermal atoms. This redistribution cannot be treated in theories that are still of first order in the binary interaction, such as, for example, Hartree-Fock-Bogoliubov (HFB) theories.

We present a kinetic theory of second order in the interaction, which is formulated in terms of Bogoliubov quasiparticles and contains collisional terms beyond the HFB approximation. The HFB interactions are expressed as many-body $`T`$ matrices to second order in the binary potential and thus include the correct renormalized scattering physics. This theory thus contains no ultraviolet divergences and has a gapless energy spectrum. We extend the papers of Walser et al. Walser1999a ; Walser2001a , which are the basis of the present results, by these two essential aspects.

We go beyond theoretical approaches that drop the anomalous pair matrix $`\stackrel{~}{m}=\stackrel{~}{a}\stackrel{~}{a}`$ in the Popov approximation Zaremba1999a ; AlKhawaja2000a ; Imamovic2000b ; Graham2000b . Monte Carlo simulations based on the semi-classical Zaremba-Nikuni-Griffin theory Zaremba1999a  show very good agreement for experimentally observed damping rates and response frequencies Jackson2001a ; Jackson2002c . However, recent experiments Donley2002a  and their theoretical explanations Kokkelmans2002b ; Koehler2002b  have shown that the pairing field plays an important role in Bose gases with resonance interactions.

Keeping the pair matrix $`\stackrel{~}{m}`$ in this theory would cause ultraviolet divergences if we replaced the bare interaction potentials with the $`s`$-wave scattering length $`a_s`$ in the contact-potential approximation. The vacuum part of the pairing field’s self interaction, for example, will diverge, because the delta-function potential contains unphysically high energy contributions. These divergences can be resolved by writing the interactions in terms of scattering $`T`$ matrices, which subsume the divergent sums over intermediate scattering states and correctly reduce to the scattering length $`a_s`$ in the zero-energy and -momentum limit Proukakis1998a ; Stoof1999a ; Hutchinson2000a ; Morgan2000a ; Koehler2002a . Using these scattering $`T`$ matrices, we can consistently eliminate all divergences to second order in the interaction.

We thus obtain a renormalized HFB operator, which has a gapless spectrum. The zero-energy eigenspace is spanned by the condensate, and excitations with non-vanishing energy are thus automatically orthogonal to the condensate. This is similar to the renormalized gapless HFB equations proposed in Ref. Hutchinson2000a . Other approaches have to explicitly project the excitations orthogonal to the condensate Morgan2000a ; Castin1997a . Using the condensate as the ground state of this adiabatic basis also simplifies the complicated representation of a condensate band in the Gardiner-Zoller master equation formulation Davis2000a .

The Monte Carlo simulations of Jackson and Zaremba Jackson2001a ; Jackson2002c  show that considering dynamic population-exchange between the condensate and the thermal cloud and within the thermal cloud leads to good agreement with the observed response spectra. We include these kinetic effects, thus going beyond collisionless descriptions Giorgini2000a ; Fedichev1998b .

Our presentation builds on the papers by Walser et al. Walser1999a ; Walser2001a  and we will begin by summarizing the kinetic equations derived in these papers in Sec. II. In Sec. III, we examine the first-order HFB propagator and rewrite it in terms of second-order $`T`$ matrices, by including second-order energy shifts and adiabatically eliminating the pairing field $`\stackrel{~}{m}`$. This shows that the theory is explicitly gapless and renormalized. We then make use of the gaplessness and write the kinetic equations in terms of Bogoliubov quasiparticles, which are orthogonal to the condensate by construction. The results in Sec. IV show that the complicated collision terms presented in the Walser et al. papers can be simplified dramatically by a basis transformation. Practical calculations of the quantum Boltzmann equation then will require only diagonal elements of the quasiparticle population matrix.

## II Single-Particle Kinetic Equations

We present the Walser et al. formulation of kinetic theory Walser2001a , which was originally derived using a statistical operator approach Walser1999a . The authors of this paper also recently showed a derivation of the same equations Wachter2001a  from the Kadanoff-Baym Kadanoff ; Kane1965a  theory of non-equilibrium Green functions. They used the gapless second-order Beliaev approximation Imamovic2000b ; Beliaev1958b  and showed equivalence to the work of Walser et al. under the Markov approximation. Recently, another independent derivation Proukakis2001a ; Proukakis2002a  connected this approach to kinetic theories by Morgan and Proukakis.

Neglecting three-body and higher interactions, we can describe the weakly interacting, dilute gas by the following Hamiltonian

$$\widehat{H}=H_{}^{(0)}{}_{}{}^{1^{}2^{}}\widehat{a}_1^{}^{}\widehat{a}_2^{}^{}+\varphi ^{1^{}2^{}3^{}4^{}}\widehat{a}_1^{}^{}\widehat{a}_2^{}^{}\widehat{a}_3^{}^{}\widehat{a}_4^{}^{},$$
(1)
where $`H_{}^{(0)}{}_{}{}^{1^{}2^{}}=1^{}\left|\widehat{H}^{(0)}\right|2^{}`$ denotes the matrix elements of the interaction-free single-particle Hamiltonian

$$\widehat{H}^{(0)}=\frac{\widehat{𝐩}^2}{2m}+V_{\mathrm{ext}}(\widehat{𝐱^{}})$$
(2)
with external harmonic potential $`V_{\mathrm{ext}}`$. The bosonic creation operator $`\widehat{a}_1^{}^{}`$ creates a particle in the state $`|1^{}`$, where $`1^{}`$ stands for a complete set of quantum numbers, which label a constant, single-particle energy basis, such as, for example, harmonic oscillator states or eigenstates of the GP equation. We use the summation convention for these abbreviated indices and indicate the single-particle basis with primes.

The two-particle matrix elements of the binary interaction potential $`V_{\mathrm{bin}}(𝐱_1^{},𝐱_2^{})`$ are defined by

$`\varphi ^{1^{}2^{}3^{}4^{}}={\displaystyle \frac{1}{4}}`$ $`{\displaystyle }`$ $`d𝐱_1^{}d𝐱_2^{}1^{}𝐱_1^{}2^{}𝐱_2^{}V_{\mathrm{bin}}(𝐱_1^{},𝐱_2^{})`$ (3)
$`\times `$ $`\left\{𝐱_1^{}3^{}𝐱_2^{}4^{}+𝐱_1^{}4^{}𝐱_2^{}3^{}\right\}.`$
These matrix elements are symmetric in the first and last two indices:

$$\varphi ^{1^{}2^{}3^{}4^{}}=\varphi ^{2^{}1^{}3^{}4^{}}=\varphi ^{1^{}2^{}4^{}3^{}}.$$
(4)
To determine the measurable quantities we want to calculate in this theory, we first define the condensate mean field $`\alpha `$ as the expectation value of the destruction operator

$$\alpha =\alpha _1^{}^{}|1^{}\widehat{a}_1^{}^{}|1^{}.$$
(5)
The total density matrix $`f`$ is defined by

$$f\widehat{a}_2^{}^{}\widehat{a}_1^{}^{}|1^{}2^{}|.$$
(6)
Subtracting the condensate density matrix

$$f^c\alpha \alpha ^{}\alpha _2^{}^{}\alpha _1^{}^{}|1^{}2^{}|,$$
(7)
we obtain the density matrix of thermal atoms $`\stackrel{~}{f}ff^c`$. The anomalous average $`m`$ is split analogously

$$m=\widehat{a}_1^{}^{}\widehat{a}_2^{}^{}|1^{}|2^{}=m^c+\stackrel{~}{m},$$
(8)
in the condensate part $`m^c\alpha \alpha `$ and fluctuations $`\stackrel{~}{m}`$.

This set of variables contains all possible combinations of up to two field operators. The course-grained statistical operator parametrized by these variables is thus Gaussian, and we can use Wick’s theorem to truncate the coupling to higher-order correlation functions, i. e., expectation values of more then two field operators. This approximation is valid because of the diluteness of the condensed gas. In a dilute gas the duration of a collision event is short compared to the essentially free evolution between collisions, which allows higher-order correlations to damp.

The procedure followed by Walser et al. Walser1999a  is then to write the Heisenberg equations of motion for the variables above and expand the expectation values using Wick’s theorem. Walser et al. thus obtain the equations of motion given in the following two sections Walser2001a .

### II.1 Mean-field equations

Since the anomalous fluctuations $`\stackrel{~}{m}`$ couple the mean field $`\alpha `$ to its conjugate $`\alpha ^{}=\alpha _1^{}^{}1^{}|`$, it is convenient to write the generalized Gross-Pitaevskii (GP) equation in a two-by-two matrix form

$$\frac{d}{dt}\chi =(i\mathrm{\Pi }+\mathrm{{\rm Y}}^<\mathrm{{\rm Y}}^>)\chi ,$$
(9)
where the two-component state vector

$$\chi =\left(\begin{array}{c}\alpha \\ \alpha ^{}\end{array}\right)$$
(10)
is defined in terms of $`\alpha =\alpha _1^{}^{}|1^{}`$ and also contains the time-reversed mean field $`\alpha ^{}`$.

The generalized GP propagator representing the reversible evolution of $`\chi `$ is defined as

$$\mathrm{\Pi }=\left(\begin{array}{cc}\mathrm{\Pi }_𝒩& \mathrm{\Pi }_𝒜\\ \mathrm{\Pi }_𝒜^{}& \mathrm{\Pi }_𝒩^{}\end{array}\right).$$
(11)
This symplectic propagator consists of the normal Hermitian Hamiltonian

$$\mathrm{\Pi }_𝒩=\widehat{H}^{(0)}+U_{f^c}+2U_{\stackrel{~}{f}}\mu ,$$
(12)
where $`\mu `$ removes rapid oscillations of the mean field, and the symmetric anomalous coupling

$$\mathrm{\Pi }_𝒜=V_{\stackrel{~}{m}}.$$
(13)
The energy shifts due to both the mean field and the normal fluctuations are given by the matrices

$$U_f=2\varphi ^{1^{}2^{}3^{}4^{}}f_{3^{}2^{}}^{}|1^{}4^{}|,$$
(14)
whereas the first-order anomalous coupling-strength is given by

$$V_m=2\varphi ^{1^{}2^{}3^{}4^{}}m_{3^{}4^{}}^{}|1^{}|2^{}.$$
(15)
The second-order irreversible evolution, consisting of damping rates and energy shifts, is given by the collision operator

$$\mathrm{{\rm Y}}^<=\left(\begin{array}{cc}\mathrm{{\rm Y}}_𝒩^<& \mathrm{{\rm Y}}_𝒜^<\\ \mathrm{{\rm Y}}_𝒜^>& \mathrm{{\rm Y}}_𝒩^>\end{array}\right)$$
(16)
and its time-reversed counterpart $`\mathrm{{\rm Y}}^>=\sigma _1\mathrm{{\rm Y}}^<\sigma _1`$, where $`\sigma _1`$ is the first Pauli matrix exchanging the positive- and negative-energy components of vectors. The matrix elements of the collision operator are given by

$`\mathrm{{\rm Y}}_𝒩^<`$ $`=`$ $`\mathrm{\Gamma }_{\stackrel{~}{f}\stackrel{~}{f}(1+\stackrel{~}{f})}+2\mathrm{\Gamma }_{\stackrel{~}{f}\stackrel{~}{m}\stackrel{~}{m}^{}},`$ (17a)
$`\mathrm{{\rm Y}}_𝒩^>`$ $`=`$ $`\mathrm{\Gamma }_{(1+\stackrel{~}{f})(1+\stackrel{~}{f})\stackrel{~}{f}}+2\mathrm{\Gamma }_{(1+\stackrel{~}{f})\stackrel{~}{m}\stackrel{~}{m}^{}},`$ (17b)
$`\mathrm{{\rm Y}}_𝒜^<`$ $`=`$ $`\mathrm{\Gamma }_{\stackrel{~}{m}\stackrel{~}{m}\stackrel{~}{m}^{}}+2\mathrm{\Gamma }_{\stackrel{~}{f}\stackrel{~}{m}(1+\stackrel{~}{f})},`$ (17c)
$`\mathrm{{\rm Y}}_𝒜^>`$ $`=`$ $`\mathrm{\Gamma }_{\stackrel{~}{m}\stackrel{~}{m}\stackrel{~}{m}^{}}+2\mathrm{\Gamma }_{(1+\stackrel{~}{f})\stackrel{~}{m}\stackrel{~}{f}},`$ (17d)
which is defined in terms of individual collisions $`\mathrm{\Gamma }`$. These elementary collision processes are defined explicitly as

$`\mathrm{\Gamma }_{fff}`$ $`=`$ $`8\varphi ^{1^{}2^{}3^{}4^{}}\varphi _\eta ^{1^{\prime \prime }2^{\prime \prime }3^{\prime \prime }4^{\prime \prime }}f_{3^{}1^{\prime \prime }}^{}f_{4^{}2^{\prime \prime }}^{}f_{4^{\prime \prime }2^{}}^{}|1^{}3^{\prime \prime }|,`$ (18a)
$`\mathrm{\Gamma }_{fmf}`$ $`=`$ $`8\varphi ^{1^{}2^{}3^{}4^{}}\varphi _\eta ^{1^{\prime \prime }2^{\prime \prime }3^{\prime \prime }4^{\prime \prime }}f_{3^{}2^{\prime \prime }}^{}m_{4^{}3^{\prime \prime }}^{}f_{4^{\prime \prime }2^{}}^{}|1^{}|1^{\prime \prime },`$ (18b)
$`\mathrm{\Gamma }_{fmm^{}}`$ $`=`$ $`8\varphi ^{1^{}2^{}3^{}4^{}}\varphi _\eta ^{1^{\prime \prime }2^{\prime \prime }3^{\prime \prime }4^{\prime \prime }}f_{3^{}1^{\prime \prime }}^{}m_{4^{}4^{\prime \prime }}^{}m_{2^{\prime \prime }2^{}}^{}|1^{}3^{\prime \prime }|,`$ (18c)
$`\mathrm{\Gamma }_{mmm^{}}`$ $`=`$ $`8\varphi ^{1^{}2^{}3^{}4^{}}\varphi _\eta ^{1^{\prime \prime }2^{\prime \prime }3^{\prime \prime }4^{\prime \prime }}m_{3^{}4^{\prime \prime }}^{}m_{4^{}3^{\prime \prime }}^{}m_{2^{\prime \prime }2^{}}^{}|1^{}|1^{\prime \prime }.`$ (18d)
The in-rates of the collision operator $`\mathrm{{\rm Y}}`$ are depicted in Fig. 1.

In this theory, collisional interactions are considered to second order. The effect of higher order terms, which lead to a finite duration of a collision, can be modeled by introducing a parameter $`\eta `$, such that every second-order collision operator contains dispersive as well as dissipative parts from the complex-valued matrix element

$$\varphi _\eta ^{1^{\prime \prime }2^{\prime \prime }3^{\prime \prime }4^{\prime \prime }}=\varphi ^{1^{\prime \prime }2^{\prime \prime }3^{\prime \prime }4^{\prime \prime }}\frac{1}{\eta i\mathrm{\Delta }ϵ},$$
(19)
where the energy difference $`\mathrm{\Delta }ϵ=(ϵ_{1^{\prime \prime }}^0+ϵ_{2^{\prime \prime }}^0)+ϵ_{3^{\prime \prime }}^0+ϵ_{4^{\prime \prime }}^0`$ has to be smaller than the energy uncertainty $`\eta `$ to get a sizable contribution. Note that the papers Walser1999a ; Walser2001a  contain a sign error in the definition of $`\mathrm{\Delta }ϵ`$. For small $`\eta `$ we obtain

$$\frac{1}{\eta i\mathrm{\Delta }ϵ}\underset{\eta 0}{}\pi \delta _\eta (\mathrm{\Delta }ϵ)+i𝒫_\eta \frac{1}{\mathrm{\Delta }ϵ},$$
(20)
where $`𝒫`$ indicates that the Cauchy principal value has to be taken upon integration. The parameter $`\eta `$ thus represents off-the-energy-shell propagation after a collision. Most off-the-energy-shell coherences decay during subsequent propagation, but, due to the finite time between collisions $`\tau _c`$, energy cannot be conserved exactly, because $`\eta `$ has to be larger than the collision rate $`1/\tau _c`$.

### II.2 Equations for normal densities and anomalous fluctuations

The equations of motion for the fluctuation densities $`\stackrel{~}{f}`$ and $`\stackrel{~}{m}`$ are coupled and can also conveniently be written in terms of two-by-two matrices. To achieve this, we define the generalized single-time fluctuation-density matrix $`\stackrel{~}{G}^<`$ as

$$\stackrel{~}{G}^<=\left(\begin{array}{cc}\stackrel{~}{f}& \stackrel{~}{m}\\ \stackrel{~}{m}^{}& (1+\stackrel{~}{f})^{}\end{array}\right),$$
(21)
where $`\stackrel{~}{f}=\stackrel{~}{f}_{1^{}2^{}}|1^{}2^{}|`$ and $`\stackrel{~}{m}=\stackrel{~}{m}_{1^{}2^{}}|1^{}|2^{}`$ are the matrix representations of the master variables. In Ref. Wachter2001a , we use the property that this density matrix is the single-time limit of the corresponding time-ordered two-time Green’s function. This showed that the other time ordering is given by

$$\stackrel{~}{G}^>=\left(\begin{array}{cc}(1+\stackrel{~}{f})& \stackrel{~}{m}\\ \stackrel{~}{m}^{}& \stackrel{~}{f}^{}\end{array}\right)=\sigma _1\stackrel{~}{G}^<\sigma _1=\sigma _3+\stackrel{~}{G}^<.$$
(22)
Here, we use the third Pauli matrix $`\sigma _3=\mathrm{diag}(𝟙,𝟙)`$. Note that our naming of the fluctuation-density matrices $`\stackrel{~}{G}^<`$ and $`\stackrel{~}{G}^>`$ is consistent with the two-time formalism in Ref. Wachter2001a , but disagrees with Ref. Walser2001a .

The generalized Boltzmann equation of motion for this fluctuation-density matrix can be written as

$$\frac{d}{dt}\stackrel{~}{G}^<=i\mathrm{\Sigma }\stackrel{~}{G}^<+\mathrm{\Gamma }^<\stackrel{~}{G}^>\mathrm{\Gamma }^>\stackrel{~}{G}^<+\mathrm{H}.\mathrm{c}.$$
(23)
This equation has to be solved under the constraints

$$\alpha ^{}\stackrel{~}{f}=0\text{and}\alpha ^{}\stackrel{~}{m}=0,$$
(24)
which force the fluctuations to be orthogonal to the condensate.

Again, the equation of motion (23) has two parts: The reversible evolution is governed by the Hartree-Fock-Bogoliubov self-energy operator

$$\mathrm{\Sigma }=\left(\begin{array}{cc}\mathrm{\Sigma }_𝒩& \mathrm{\Sigma }_𝒜\\ \mathrm{\Sigma }_𝒜^{}& \mathrm{\Sigma }_𝒩^{}\end{array}\right),$$
(25)
which in turn consists of the Hermitian Hamiltonian

$$\mathrm{\Sigma }_𝒩=\widehat{H}^{(0)}+2U_{f^c}+2U_{\stackrel{~}{f}}\mu $$
(26)
and the symmetric anomalous coupling

$$\mathrm{\Sigma }_𝒜=V_{m^c}+V_{\stackrel{~}{m}}.$$
(27)
The irreversible evolution introduced by second-order collisional contributions now consists of the collisional operator

$$\mathrm{\Gamma }^<=\left(\begin{array}{cc}\mathrm{\Gamma }_𝒩^<& \mathrm{\Gamma }_𝒜^<\\ \mathrm{\Gamma }_𝒜^>& \mathrm{\Gamma }_𝒩^>\end{array}\right),$$
(28)
and its time-reversed counter part $`\mathrm{\Gamma }^>=\sigma _1\mathrm{\Gamma }^<\sigma _1`$. The diagonal components of the collision operator are defined as

$`\mathrm{\Gamma }_𝒩^<`$ $`=`$ $`\mathrm{\Gamma }_{(\stackrel{~}{f}+f^c)\stackrel{~}{f}(1+\stackrel{~}{f})}+\mathrm{\Gamma }_{\stackrel{~}{f}f^c(1+\stackrel{~}{f})}+\mathrm{\Gamma }_{\stackrel{~}{f}\stackrel{~}{f}f^c}`$
$`+2\left\{\mathrm{\Gamma }_{(\stackrel{~}{f}+f^c)\stackrel{~}{m}\stackrel{~}{m}^{}}+\mathrm{\Gamma }_{\stackrel{~}{f}m^c\stackrel{~}{m}^{}}+\mathrm{\Gamma }_{\stackrel{~}{f}\stackrel{~}{m}m^c}\right\},`$
$`\mathrm{\Gamma }_𝒩^>`$ $`=`$ $`\mathrm{\Gamma }_{(1+\stackrel{~}{f}+f^c)(1+\stackrel{~}{f})\stackrel{~}{f}}+\mathrm{\Gamma }_{(1+\stackrel{~}{f})f^c\stackrel{~}{f}}+\mathrm{\Gamma }_{(1+\stackrel{~}{f})(1+\stackrel{~}{f})f^c}`$
$`+2\left\{\mathrm{\Gamma }_{(1+\stackrel{~}{f}+f^c)\stackrel{~}{m}\stackrel{~}{m}^{}}+\mathrm{\Gamma }_{(1+\stackrel{~}{f})m^c\stackrel{~}{m}^{}}+\mathrm{\Gamma }_{(1+\stackrel{~}{f})\stackrel{~}{m}m^c}\right\},`$
and the off-diagonal, anomalous components as

$`\mathrm{\Gamma }_𝒜^<`$ $`=`$ $`\mathrm{\Gamma }_{(\stackrel{~}{m}+m^c)\stackrel{~}{m}\stackrel{~}{m}^{}}+\mathrm{\Gamma }_{\stackrel{~}{m}m^c\stackrel{~}{m}^{}}+\mathrm{\Gamma }_{\stackrel{~}{m}\stackrel{~}{m}m^c}`$
$`+2\left\{\mathrm{\Gamma }_{(\stackrel{~}{f}+f^c)\stackrel{~}{m}(1+\stackrel{~}{f})}+\mathrm{\Gamma }_{\stackrel{~}{f}m^c(1+\stackrel{~}{f})}+\mathrm{\Gamma }_{\stackrel{~}{f}\stackrel{~}{m}f^c}\right\},`$
$`\mathrm{\Gamma }_𝒜^>`$ $`=`$ $`\mathrm{\Gamma }_{(\stackrel{~}{m}+m^c)\stackrel{~}{m}\stackrel{~}{m}^{}}+\mathrm{\Gamma }_{\stackrel{~}{m}m^c\stackrel{~}{m}^{}}+\mathrm{\Gamma }_{\stackrel{~}{m}\stackrel{~}{m}m^c}`$
$`+2\left\{\mathrm{\Gamma }_{(1+\stackrel{~}{f}+f^c)\stackrel{~}{m}\stackrel{~}{f}}+\mathrm{\Gamma }_{(1+\stackrel{~}{f})m^c\stackrel{~}{f}}+\mathrm{\Gamma }_{(1+\stackrel{~}{f})\stackrel{~}{m}f^c}\right\}.`$
Both the diagonal and off-diagonal incoming rates are depicted diagrammatically in Fig. 2. For every term that appears in the collisional terms of the generalized GP equation in $`\mathrm{{\rm Y}}^<`$ \[Fig. (1)\], we here \[Fig. (2)\] have three additional terms, where in each of them, one of the three fluctuating contributions is replaced with the corresponding mean-field quantity. This replacement rule can be seen in the Beliaev collisional self energies presented in Ref. Wachter2001a  and is a consequence of the fact that the Boltzmann equation (23) can be generated from the GP equation (9) by functional differentiation Hohenberg1965a .

When the collision operator $`\mathrm{\Gamma }^<`$ is multiplied by $`\stackrel{~}{G}^>`$ as in Eq. (23), we get terms like

$$\mathrm{\Gamma }_{\stackrel{~}{f}\stackrel{~}{f}(1+\stackrel{~}{f})}(1+\stackrel{~}{f})\mathrm{\Gamma }_{(1+\stackrel{~}{f})(1+\stackrel{~}{f})\stackrel{~}{f}}\stackrel{~}{f},$$
(33)
where the second part comes from the time-reversed contribution $`\mathrm{\Gamma }^>\stackrel{~}{G}^<`$. The diagonal parts are exactly the in and out terms of the quantum Boltzmann equation for the single-particle distribution function $`\stackrel{~}{f}`$. The remaining second-order contributions couple to the anomalous fluctuations $`\stackrel{~}{m}`$ and do not have an analogue in the quantum Boltzmann equation. In Sec. IV, we will rewrite the kinetic equations presented here in terms of Bogoliubov quasiparticles. Then all collisional contribution take the form of Boltzmann terms.

## III Gaplessness — $`T`$ matrices

Our goal in this section is to explicitly show that the kinetic equations (23) and (9) are gapless. This should on one hand be obvious, because our previous paper Wachter2001a  showed them to be equivalent to the Kadanoff-Baym equations Kadanoff ; Kane1965a  in the gapless Beliaev approximation Imamovic2000b ; Beliaev1958b . The first-order HFB propagator $`\mathrm{\Sigma }`$, which appears in the kinetic equations, is, on the other hand, known to exhibit a non-physical energy gap in the long-wavelength, homogeneous limit Griffin1996a .

We resolve this discrepancy by including second-order collisional energy shifts $`𝒫\{\mathrm{\Gamma }\}`$ into the HFB propagator and adiabatically eliminating the anomalous average $`\stackrel{~}{m}`$ in the first-order anomalous potential $`V_{\stackrel{~}{m}}`$ (15). This upgrades the bare interaction potentials in the first-order propagators $`\mathrm{\Sigma }`$ and $`\mathrm{\Pi }`$ to the real parts of many-body $`T`$ matrices. We can then approximate these $`T`$ matrices with the $`s`$-wave scattering length without incurring ultra-violet divergences.

The upgraded HFB propagator $`\mathrm{\Sigma }`$, where all binary interactions are written as many-body $`T`$ matrices, is explicitly gapless and thus obeys the Hugenholtz-Pines theorem Hugenholtz1959a . The propagator thus has zero-energy modes, which are completely specified by the value of the condensate $`\alpha `$. If we use the non-zero energy Bogoliubov modes of $`\mathrm{\Sigma }`$ as a basis for the thermal excitations, the excitations will automatically be orthogonal to the condensate. We follow this idea in Sec. IV.

### III.1 Off-diagonal potentials

Here, we update the off-diagonal potentials $`V_{(m^c+\stackrel{~}{m})}`$ in the HFB propagator $`\mathrm{\Sigma }`$ by adiabatically eliminating the pairing field $`\stackrel{~}{m}`$. We integrate the first-order equation of motion for the anomalous average $`\stackrel{~}{m}`$

$$\frac{d}{dt}\stackrel{~}{m}=i\mathrm{\Sigma }_𝒩\stackrel{~}{m}i\stackrel{~}{m}\mathrm{\Sigma }_𝒩i\mathrm{\Sigma }_𝒜(1+\stackrel{~}{f})^{}i\stackrel{~}{f}\mathrm{\Sigma }_𝒜,$$
(34)
which is obtained by taking the $`\stackrel{~}{m}`$ component of the generalized Boltzmann equation (23) and dropping the second-order terms, because we want to substitute the result for $`\stackrel{~}{m}`$ into the anomalous potential $`V_{\stackrel{~}{m}}`$ and only keep terms up to second order. In stationarity, we solve for $`\stackrel{~}{m}`$ in the eigenbasis of $`\mathrm{\Sigma }_𝒩`$,

$$\mathrm{\Sigma }_𝒩|ϵ_1^{}^0=(ϵ_1^{}^0\mu )|ϵ_1^{}^0,$$
(35)
and obtain

$$\stackrel{~}{m}_{1^{}2^{}}^{}=𝒫\frac{\mathrm{\Sigma }_𝒜^{1^{}2^{\prime \prime }}(1+\stackrel{~}{f})_{2^{}2^{\prime \prime }}^{}+\stackrel{~}{f}_{1^{}2^{\prime \prime }}^{}\mathrm{\Sigma }_𝒜^{2^{\prime \prime }2^{}}}{2\mu (ϵ_1^{}^0+ϵ_2^{}^0)}$$
(36)
as an adiabatic solution. Adiabatic here means that this solution only includes time-variations with characteristic times long compared to the duration of a collision. We use the Cauchy principal value $`𝒫`$ to indicate omission of the divergent term in an energy integral or sum. This divergent $`\delta `$-function term gives rise to the imaginary part. We plug this result into the off-diagonal potential (15),

$$V_{\stackrel{~}{m}}^{1^{}2^{}}=4𝒫\frac{\varphi ^{1^{}2^{}3^{\prime \prime }4^{\prime \prime }}(1+2\stackrel{~}{f})_{4^{\prime \prime }2^{\prime \prime }}^{}\varphi ^{3^{\prime \prime }2^{\prime \prime }3^{}4^{}}}{2\mu (ϵ_{3^{\prime \prime }}^0+ϵ_{4^{\prime \prime }}^0)}m_{3^{}4^{}}^c,$$
(37)
where we dropped the recursive $`V_{\stackrel{~}{m}}`$ term in the anomalous coupling $`\mathrm{\Sigma }_𝒜`$ in order to keep Eq. (37) at second order. We will discuss the recursive term in Sec. III.4.

We then recognize that we can write the off-diagonal element of the HFB propagator $`\mathrm{\Sigma }_𝒜`$ as the real part of a many-body $`T`$ matrix

$`\mathrm{\Sigma }_𝒜`$ $`=`$ $`V_{m^c}+V_{\stackrel{~}{m}}=T_{m^c}(2\mu ),`$ (38)
which is defined by

$$T^{1^{}2^{}3^{}4^{}}(ϵ)=2\varphi ^{1^{}2^{}3^{}4^{}}+4𝒫\frac{\varphi ^{1^{}2^{}3^{\prime \prime }4^{\prime \prime }}(1+2\stackrel{~}{f})_{4^{\prime \prime }2^{\prime \prime }}^{}\varphi ^{3^{\prime \prime }2^{\prime \prime }3^{}4^{}}}{ϵ(ϵ_{3^{\prime \prime }}^0+ϵ_{4^{\prime \prime }}^0)}.$$
(39)
The energies $`ϵ^0`$ are dressed by the normal and mean-field shifts, but are not the full quasiparticle energies, because they do not include the effect of the pairing field. Contractions of this $`T`$ matrix with anomalous averages are defined by

$$T_m^{1^{}2^{}}(ϵ)=T^{1^{}2^{}3^{}4^{}}(ϵ)m_{3^{}4^{}}^{}.$$
(40)
The $`T`$ matrix defined in Eq. (39) is a function of energy through its last two indices in the sense that its argument $`ϵ=ϵ_3^{}^0+ϵ_4^{}^0`$.

### III.2 Diagonal Potentials

In this Section, we want to redefine the diagonal potentials $`U_{f^c}`$ and $`U_{\stackrel{~}{f}}`$ as the real parts of $`T`$ matrices by using the second-order energy shifts $`𝒫\{\mathrm{\Gamma }\}`$. With $`𝒫\{\mathrm{\Gamma }\}`$ we here denote the principal-value part of the collisional terms in Eqs. (28) and (16) according to Eq. (20). We will begin by considering the condensate potential.

$`i\mathrm{\hspace{0.33em}2}U_{f^c}+𝒫\{2\mathrm{\Gamma }_{f^c\stackrel{~}{f}(1+\stackrel{~}{f})}+\mathrm{\Gamma }_{\stackrel{~}{f}\stackrel{~}{f}f^c}`$
$`2\mathrm{\Gamma }_{f^c(1+\stackrel{~}{f})\stackrel{~}{f}}\mathrm{\Gamma }_{(1+\stackrel{~}{f})(1+\stackrel{~}{f})f^c}\}`$
$`=i\mathrm{\hspace{0.33em}2}U_{f^c}𝒫\{\mathrm{\Gamma }_{1(1+2\stackrel{~}{f})f^c}\}`$ (41)
We here assume real eigenfunctions for the single-particle energy basis and do not include any $`\mathrm{\Gamma }`$ terms involving the anomalous average $`\stackrel{~}{m}`$, because they are at least of order $`V_{\mathrm{bin}}^3/(\mathrm{\Delta }ϵ)^2`$ according to Eq. (36). The second-order terms that contain only normal fluctuations $`\stackrel{~}{f}`$ are used in Eq. (46) to rewrite the fluctuation potential $`U_{\stackrel{~}{f}}`$. The term in Eq. (41) can again be written in terms of a many-body $`T`$ matrix

$$U_{f^c}+\frac{1}{2i}𝒫\{\mathrm{\Gamma }_{1(1+2\stackrel{~}{f})f^c}\}=T_{f^c},$$
(42)
which is given by

$$T^{1^{}2^{}3^{}4^{}}(ϵ)=2\varphi ^{1^{}2^{}3^{}4^{}}+4𝒫\frac{\varphi ^{1^{}2^{}3^{\prime \prime }4^{\prime \prime }}(1+2\stackrel{~}{f})_{4^{\prime \prime }2^{\prime \prime }}^{}\varphi ^{3^{\prime \prime }2^{\prime \prime }3^{}4^{}}}{ϵ(ϵ_{3^{\prime \prime }}^0+ϵ_{2^{\prime \prime }}^0)}.$$
(43)
The slight difference compared to Eq. (39) is resolved when we assume diagonal quasiparticle populations $`P_{\overline{1}\overline{2}}=P_{\overline{1}}\delta _{\overline{1}\overline{1}}`$ as will be justified in Sec. IV:

$$\stackrel{~}{f}_{1^{}2^{}}^{}=U_1^{}^{\overline{1}}P_{\overline{1}}U_2^{}^{\overline{1}}=\stackrel{~}{f}_{2^{}1^{}}^{},$$
(44)
where $`U`$ is the transformation matrix to the quasiparticle basis. Alternatively, we note that the $`T`$ matrix is essentially constant for energy differences up to the duration of a collision. Contractions of the $`T`$ matrix with normal averages are performed according to

$$T_f^{1^{}4^{}}T^{1^{}2^{}3^{}4^{}}(ϵ_3^{}^0+ϵ_4^{}^0)f_{3^{}2^{}}.$$
(45)
We now consider the fluctuation potential $`U_{\stackrel{~}{f}}`$, again include the truly second-order energy shifts, and obtain

$`i\mathrm{\hspace{0.33em}2}U_{\stackrel{~}{f}}+𝒫\{\mathrm{\Gamma }_{\stackrel{~}{f}\stackrel{~}{f}(1+\stackrel{~}{f})}\mathrm{\Gamma }_{(1+\stackrel{~}{f})(1+\stackrel{~}{f})\stackrel{~}{f}}\}`$
$`=i\mathrm{\hspace{0.33em}2}U_{\stackrel{~}{f}}𝒫\{\mathrm{\Gamma }_{1(1+\stackrel{~}{f})\stackrel{~}{f}}\}=T_{\stackrel{~}{f}}^{},`$ (46)
where we get a different $`T`$ matrix defined by

$$T^{1^{}2^{}3^{}4^{}}(ϵ)=2\varphi ^{1^{}2^{}3^{}4^{}}+4𝒫\frac{\varphi ^{1^{}2^{}3^{\prime \prime }4^{\prime \prime }}(1+\stackrel{~}{f})_{4^{\prime \prime }2^{\prime \prime }}^{}\varphi ^{3^{\prime \prime }2^{\prime \prime }3^{}4^{}}}{ϵ(ϵ_{3^{\prime \prime }}^0+ϵ_{2^{\prime \prime }}^0)},$$
(47)
which does not have a factor of $`2`$ in the intermediate-population term $`(1+\stackrel{~}{f})`$. This difference is due to the fact that the mean field $`\alpha `$ is not bosonically enhanced. If we assume diagonal population $`\stackrel{~}{f}_{1^{}2^{}}=\delta _{1^{}2^{}}\stackrel{~}{f}_{1^{}1^{}}`$, which is not a good approximation in this basis, we reproduce the results of Ref. Proukakis1998b  for the GP equation to second order in the interaction potential. In particular, the factor of $`2`$ in their many-body $`T`$ matrix in the term corresponding to Eq. (46) gets canceled with a negative term from adiabatically eliminating their triple average.

### III.3 Renormalized Propagators

Using the $`T`$ matrices defined in the previous sections, we can now rewrite the generalized GP propagator $`\mathrm{\Pi }`$ given in Eq. (11) and the generalized Boltzmann propagator $`\mathrm{\Sigma }`$ given in Eq. (25).

The Hamiltonian of the GP equation is now

$$\mathrm{\Pi }_𝒩^{}=\widehat{H}^{(0)}+T_{f^c}^{}+2T_{\stackrel{~}{f}}^{}\mu ,$$
(48)
and the anomalous coupling $`\mathrm{\Pi }_𝒜^{}`$ vanishes, because of the identity

$$V_{\stackrel{~}{m}}\alpha ^{}=\left\{T_{m^c}(2\mu )V_{m^c}\right\}\alpha ^{}=\left\{T_{f^c}(2\mu )U_{f^c}\right\}\alpha .$$
(49)
This means that the coupling between $`\alpha `$ and $`\alpha ^{}`$ vanishes and they are no longer independent quantities, and that we can without loss of generality treat $`\alpha `$ as real.

When we write the Boltzmann propagator in terms of $`T`$ matrices, we can explicitly show that the energy spectrum is gapless; this theory thus fulfills the Hugenholtz-Pines theorem Hugenholtz1959a . The diagonal part of the propagator $`\mathrm{\Sigma }`$ is now

$$\mathrm{\Sigma }_𝒩^{}=\widehat{H}^{(0)}+2T_{f^c}^{}+2T_{\stackrel{~}{f}}^{}\mu ,$$
(50)
and the anomalous coupling is

$$\mathrm{\Sigma }_𝒜^{}=T_{m^c}(2\mu ).$$
(51)
We now show that this renormalized Boltzmann propagator $`\mathrm{\Sigma }^{}`$ has a zero-energy eigenvector, i. e., its spectrum is gapless. We begin by writing the generalized GP equation for the condensate ground state

$$\left\{\widehat{H}^{(0)}+T_{f^c}^{}(2\mu )+2T_{\stackrel{~}{f}}^{}(2\mu )\right\}\alpha =\mu \alpha .$$
(52)
This equation can be written in terms of the renormalized GP Hamiltonian (48) as $`\mathrm{\Pi }_𝒩\alpha =0`$, where the energies are now measured relative to the adiabatic chemical potential $`\mu `$. An immediate consequence of Eq. (52) is that the quasiparticle ground state

$$P_\alpha =C\left(\begin{array}{c}\alpha \\ \alpha ^{}\end{array}\right),$$
(53)
with a normalization constant $`C`$, is a zero-energy eigenvector of the renormalized $`\mathrm{\Sigma }^{}`$, because of the identity

$$T_{m^c}(2\mu )\alpha ^{}=T_{f^c}(2\mu )\alpha .$$
(54)
These zero-energy eigenvectors of the HFB propagator are proportional to the condensate mean field $`\alpha `$. All non-zero-energy eigenvectors $`W_{E_{\overline{1}}0}`$ of $`\mathrm{\Sigma }`$ are thus automatically orthogonal to the condensate, and we can use the complete set $`\{\alpha ,W_{E_{\overline{1}}0}\}`$ as a basis to describe the condensate interacting with thermal excitations Fetter1972a . Other approaches to finite-temperature theories Morgan2000a ; Castin1997a  have to explicitly orthogonalize their bases using projection operators. We discuss this new basis in Sec. IV.

### III.4 Ladder approximation

We can extend the second-order $`T`$ matrices introduced in the previous sections to include ladder diagrams to all orders by using consistency arguments. We first note that by keeping the recursive $`V_{\stackrel{~}{m}}`$ term on the right side of Eq. (37), we can extend our definition of the off-diagonal $`T`$ matrix to

$$T^{1^{}2^{}3^{}4^{}}(ϵ)=2\varphi ^{1^{}2^{}3^{}4^{}}+4𝒫\frac{\varphi ^{1^{}2^{}3^{\prime \prime }4^{\prime \prime }}(1+2\stackrel{~}{f})_{4^{\prime \prime }2^{\prime \prime }}^{}T^{3^{\prime \prime }2^{\prime \prime }3^{}4^{}}(ϵ)}{ϵ(ϵ_{3^{\prime \prime }}^0+ϵ_{4^{\prime \prime }}^0)}.$$
(55)
This is the real part of the many-body $`T`$ matrix in the ladder approximation.

In the previous section, we showed that the HFB propagator $`\mathrm{\Sigma }`$ is gapless with the $`T`$ matrices defined to second order. If we used the ladder $`T`$ in Eq. (55) on the off-diagonal while keeping the second-order $`T`$ in Eq. (43) in the diagonal, we would find an energy gap of third order, because the cancellation in Eq. (54) only works if the two $`T`$s are identical. However, the Hugenholtz-Pines theorem tells us that the full theory should again be gapless. We thus conclude that we have to upgrade the $`T`$ matrix on the diagonal of $`\mathrm{\Sigma }`$ to the ladder approximation as well.

We would like to finish our discussion of the scattering matrices in terms of the single-particle energy basis with two remarks. First, as the Liouville-space formulation Mukamel  of density-matrix evolution shows, the scattering should really be formulated in terms of Liouville-space scattering $`𝒯`$ matrices Fano1963a , which can be expressed in terms of Hilbert-space $`T`$ matrices as

$$𝒯=T𝟙+𝟙𝕋^{}+𝕋𝕋^{}.$$
(56)
Since we only consider Hilbert-space $`T`$ matrices, we thus would miss higher-order terms of the type $`TT^{}`$, even if we included the full many-body Hilbert-space scattering matrix. The imaginary part of $`𝒯`$ gives rise to the inelastic rates in the kinetic equations in Sec. IV.

Second, since the asymptotic states in this scattering problem are typically trapped harmonic-oscillator states for dilute, trapped atomic gases, we are strictly speaking dealing with bound-state $`R`$ matrices Lane1958a  instead of $`T`$ matrices.

## IV Quasiparticle Kinetic Equations

### IV.1 Quasiparticle Basis

We now want to write the kinetic equations (23) in the Bogoliubov quasiparticle basis. In this basis, the complicated and non-linear evolution due to the HFB propagator $`\mathrm{\Sigma }^{}`$ is replaced with a simple commutator with the eigenenergies and a slow basis rotation. As this theory is gapless, the $`E0`$ quasiparticle states together with the condensate $`\alpha `$ form an orthogonal basis, i. e., the thermal fluctuations are by definition orthogonal to the condensate. Another motivation for transforming to a diagonal first-order Hamiltonian can be found in the numerical results of Walser et al.; they find that the reversible first-order evolution leaves the quasiparticle populations constant. Thus, in the quasiparticle basis, only the second-order collisional terms change the populations. A more detailed account of the transformation to the quasiparticle basis can be found in Ref. Wachter2000a .

Since the quasiparticles consist of the eigenvectors of the self-energy propagator $`\mathrm{\Sigma }^{}`$, we consider the propagator’s $`2n`$ by $`2n`$ ($`n`$ is the number of single-particle states considered) eigenvector matrix $`W`$ defined by

$$\mathrm{\Sigma }^{}W=WE$$
(57)
at each time with the diagonal quasiparticle eigenenergy matrix $`E`$, which is labeled with the quasiparticle indices $`\overline{1}E_{\overline{1}}`$. The eigenvalue equations (57) are the Bogoliubov-de-Gennes equations. We decompose their solution $`W`$ into two $`n`$ by $`2n`$ matrices $`U`$ and $`V^{}`$,

$$W=\left(\begin{array}{c}U\\ V^{}\end{array}\right)=\left(\begin{array}{cc}u_+& u_{}\\ v_+^{}& v_{}^{}\end{array}\right),$$
(58)
which are in turn split into $`n`$ by $`n`$ matrices for positive ($`u_+`$, $`v_+^{}`$) and negative quasiparticle energies ($`u_{}`$, $`v_{}^{}`$). These quasiparticle eigenenergies $`\overline{1}`$ are the column indices and the original single-particle energies $`1^{}`$ the row indices, such that $`U=\left[U_1^{}^{\overline{1}}\right]_{\{1^{},\overline{1}\}}`$.

Since the propagator $`\sigma _3\mathrm{\Sigma }^{}`$ is positive semi-definite, the eigenvalues $`E`$ in Eq. (57) are real and come in positive and negative pairs of equal magnitude Blaizot . We argued in the previous section that there exists a zero-energy eigen vector $`P_\alpha `$ defined in Eq. (53). This means that the propagator $`\mathrm{\Sigma }^{}`$ is degenerate and has a two-dimensional null space. The eigenvalue equation (57), however, does not yield the second linearly independent zero-energy eigenvector. Instead, the associated vector $`Q_\alpha `$ is given by Blaizot

$$\mathrm{\Sigma }^{}Q_\alpha =i\frac{P_\alpha }{M},$$
(59)
with a positive constant $`M`$. In our case, we find

$$Q_\alpha =iC\left(\begin{array}{c}\alpha \\ \alpha ^{}\end{array}\right),$$
(60)
up to a normalization constant $`C`$. In order to find a complete set of basis states, we now define the quadrature components of $`P_\alpha `$ and $`Q_\alpha `$:

$`W^{+0}`$ $``$ $`{\displaystyle \frac{1}{\sqrt{2}}}(P_\alpha +iQ_\alpha )=\sqrt{2}C\left(\begin{array}{c}\alpha \\ 0\end{array}\right),`$ (61a)
$`W^0`$ $``$ $`{\displaystyle \frac{1}{\sqrt{2}}}(P_\alpha iQ_\alpha )=\sqrt{2}C\left(\begin{array}{c}0\\ \alpha ^{}\end{array}\right).`$ (61b)
If we now substitute these two states $`W^{+0}`$ and $`W^0`$ into the zero-energy columns of $`W`$, we can normalize all eigenvectors using the symplectic norm

$$W^{}\sigma _3W=\sigma _3.$$
(62)
The negative-energy states thus have negative norm, and the positive-energy states positive norm. In particular, the zero-energy vector $`W^{+0}`$ has positive norm and thus belongs to the positive-energy part of $`W`$. This normalization fixes the constant $`C`$ for the zero-energy states to

$$C=\frac{1}{\sqrt{2\alpha ^{}\alpha }}\frac{1}{\sqrt{2N_c}},$$
(63)
where we defined the number of condensate atoms $`N_c`$. The completeness relation for the quasiparticle basis also has symplectic structure

$$W\sigma _3W^{}\sigma _3=𝟙_{\mathrm{𝟚}𝕟},$$
(64)
where $`𝟙_{\mathrm{𝟚}𝕟}`$ indicates the $`2n`$-dimensional unit-matrix. The symplectic structure of the orthonormalization and completeness relations Eqs. (62) and (64), i. e., the appearance of the matrix $`\sigma _3`$ in these expressions, guarantees that the transformation to the quasiparticle basis is canonical Blaizot . Canonical here means that the new quasiparticle operators obey the boson commutation relations.

To further examine the structure of the quasiparticle states $`W`$, we consider the following symmetry of the HFB propagator $`\mathrm{\Sigma }^{}`$

$$\mathrm{\Sigma }^{}=\sigma _1\mathrm{\Sigma }^{}\sigma _1,$$
(65)
which holds according to its definition in Eq. (25). This symmetry implies (Walser2001a, , IV.B.) the following relation for the quasiparticle states $`W`$:

$$W=\left(\begin{array}{cc}u_+& u_{}\\ v_+^{}& v_{}^{}\end{array}\right)=\left(\begin{array}{cc}u_+& v_+\\ v_+^{}& u_+^{}\end{array}\right).$$
(66)
This relation in Eq. (66) shows that the positive- and negative-energy eigenvectors of $`\mathrm{\Sigma }^{}`$ are not independent but related by

$$W^{\overline{1}}=\sigma _1W^{\overline{1}}.$$
(67)
Plugging $`W`$ from Eq. (66) into the full completeness relation Eq. (64), we find a completeness relation for the independent elements of $`W`$ alone:

$$\frac{1}{N_c}\alpha \alpha ^{}+\underset{\overline{1}>0}{}(u_+^{\overline{1}}u_+^{\overline{1}}v_+^{\overline{1}}v_+^{\overline{1}})=𝟙_𝕟.$$
(68)
This is the basis we want to use for the kinetic equations. The explicit split in a condensate mode $`\alpha `$ and the orthogonal fluctuation modes $`u^{\overline{1}>0}`$ and $`v^{\overline{1}>0}`$ is similar to the formalism of Ref. Fetter1972a . While the condensate mode $`\alpha `$ evolves according to Eq. (9) with the updated propagator $`\mathrm{\Pi }^{}`$, we have to find an evolution equation for the quasiparticle populations.

We thus write the generalized single-time fluctuation-density matrix $`\stackrel{~}{G}^<`$ defined in Eq. (21) in terms of quasiparticles as

$$\stackrel{~}{G}^<=WPW^{}=\left(\begin{array}{cc}UPU^{}& UPV^{}\\ V^{}PU^{}& V^{}PV^{}\end{array}\right),$$
(69)
where $`P`$ is the not necessarily diagonal $`2n`$ by $`2n`$ quasiparticle population matrix. The time-reversed density matrix $`\stackrel{~}{G}^>`$ transforms according to its definition Eq. (22). From the same equation, we can also deduce that the quasiparticle population matrix $`P`$ is Hermitian and fulfills a similar identity:

$$\sigma _3+P=\sigma _1P^{}\sigma _1=\sigma _1P^{}\sigma _1.$$
(70)
This identity implies that $`P`$ can be written as

$$P=\left(\begin{array}{cc}p& q\\ q^{}& (1+p)^{}\end{array}\right),$$
(71)
in a structure similar to $`\stackrel{~}{G}^<`$ in Eq. (21). Since the classical mean field $`\alpha `$ only undergoes stimulated emission, the ground-state factors $`P_{+0+0}=P_{00}`$ do not contain enhancement terms. Thus, the zero-energy components of $`\sigma _3`$ in Eq. (70) and similar enhancement terms associated with $`\alpha `$ actually have to be zero.

We also represent the condensate (10) by a population matrix $`P^c`$ in the quasiparticle basis

$$\chi \chi ^{}=WP^cW^{}.$$
(72)
In order to maintain orthogonality between the condensate and the thermal excitations, we have to demand an analogue to the constraint (24) in the quasiparticle basis. We here have to distinguish two cases.

In the first case, the system evolves slowly enough that the ground-state of the adiabatic basis given in Eq. (68) is the true condensate state $`\alpha (t)=\alpha `$ (real). We can then explicitly give the condensate matrix as

$$P^c=N_ce_0^{}e_0^{},$$
(73)
with a vector $`e_0^{\overline{1}}=\delta _{\overline{1},\pm 0}`$, such that all non-zero-energy components of $`P_c`$ vanish. To implement the orthogonality constraint, we then explicitly set the $`E=\pm 0`$ elements of the quasiparticle matrix $`P`$ to zero:

$$P_{0\overline{1}}=P_{\overline{1}0}=0,\text{for all }\overline{1}.$$
(74)
In the second case, the system evolves too fast for the quasiparticle basis to follow. Then, the condensate matrix contains components of non-zero energy. In this case, the more general orthogonality constraint

$$\alpha (t)|P=\alpha (t)\overline{1}P_{\overline{1}\overline{2}}\overline{2}|=0$$
(75)
has to be fulfilled. This constraint is particularly important when the condensate is coherently excited in linear response. For adiabatic evolution, Eq. (75) reduces to Eq. (74). In the next section, we will find the evolution equation for $`P`$.

### IV.2 Kinetic Equations

With these ingredients, we can now try to obtain an equation of motion for the occupation number matrix $`P`$ from the kinetic equation (23). We use Eq. (69) to substitute the fluctuation-density matrices $`\stackrel{~}{G}`$ and obtain

$`{\displaystyle \frac{d}{dt}}P`$ $`=`$ $`i[E,P]_{}+\{W^1{\displaystyle \frac{dW}{dt}}P+\mathrm{H}.\mathrm{c}.\}`$
$`+`$ $`\{W^1\mathrm{\Gamma }^<W(\sigma _3+P)W^1\mathrm{\Gamma }^>WP+\mathrm{H}.\mathrm{c}.\}.`$
We are now left with the task of transforming the collisional contributions to the quasiparticle basis. To this end, we define new two-particle matrix elements in the quasiparticle energy basis by

$$\mathrm{\Phi }^{\overline{1}\overline{2}\overline{3}\overline{4}}\varphi ^{1^{}2^{}3^{}4^{}}V_1^{}^{\overline{1}}V_2^{}^{\overline{2}}U_3^{}^{\overline{3}}U_4^{}^{\overline{4}}.$$
(77)
These quasiparticle matrix elements have the same symmetries as the original ones given in Eq. (4). Furthermore, they fulfill

$$\mathrm{\Phi }^{\overline{1}\overline{2}\overline{3}\overline{4}}\left(\mathrm{\Phi }^{\overline{1}\overline{2}\overline{3}\overline{4}}\right)^{}=\mathrm{\Phi }^{\overline{3}\overline{4}\overline{1}\overline{2}}.$$
(78)
The careful examination of symmetries in Ref. Wachter2000a  shows that the collision operator Eq. (28) reduces to terms containing the completely symmetric matrix elements

$`\psi ^{\overline{1}\overline{2}\overline{3}\overline{4}}{\displaystyle \frac{1}{\sqrt{3}}}{\displaystyle \underset{\pi }{}}\mathrm{\Phi }^{\pi (\overline{1}\overline{2}\overline{3}\overline{4})},`$ (79)
which are defined as a sum over all index permutations $`\pi `$.

The resulting kinetic equations for the quasiparticle populations $`P`$ are

$`{\displaystyle \frac{d}{dt}}P=iEP+B_WP+C_{PP}+C_{\alpha P}+\mathrm{H}.\mathrm{c}.,`$ (80)
with the non-adiabatic basis-rotation term

$$B_W=W^1\frac{dW}{dt}=\sigma _3W^{}\sigma _3\frac{dW}{dt},$$
(81)
and the Boltzmann collision rates between fluctuations

$`C_{PP}`$ $`=`$ $`C_{PP}^<(\sigma _3+P)C_{PP}^>P`$
$`=`$ $`\mathrm{\Gamma }_{PP(\sigma _3+P)}(\sigma _3+P)\mathrm{\Gamma }_{(\sigma _3+P)(\sigma _3+P)P}P`$
and between fluctuations and the condensate

$$C_{\alpha P}=3\mathrm{\Gamma }_{P^cP(\sigma _3+P)}(\sigma _3+P)3\mathrm{\Gamma }_{P^c(\sigma _3+P)P}P.$$
(83)
These collision rates are defined in terms of the quasiparticle collision operator

$$\sigma _3\mathrm{\Gamma }_{PPP}=\frac{1}{2}\psi ^{\overline{1}\overline{2}\overline{3}\overline{4}}\psi _\eta ^{\overline{1}^{}\overline{2}^{}\overline{3}^{}\overline{4}^{}}P_{\overline{1}\overline{1}^{}}^{}P_{\overline{2}\overline{2}^{}}^{}P_{\overline{4}^{}\overline{4}}^{}|\overline{3}\overline{3}^{}|.$$
(84)
In this operator, $`P`$ can stand for any one of the three possibilities $`P`$, $`(\sigma _3+P)`$, or $`P^c`$, as needed in Eqs. (IV.2) and (83). The approximately energy conserving matrix element $`\psi _\eta `$ is explicitly given by

$$\psi _\eta ^{\overline{1}^{}\overline{2}^{}\overline{3}^{}\overline{4}^{}}=\psi ^{\overline{1}^{}\overline{2}^{}\overline{3}^{}\overline{4}^{}}\pi \delta _\eta (E_{\overline{1}^{}}+E_{\overline{2}^{}}+E_{\overline{3}^{}}+E_{\overline{4}^{}}),$$
(85)
where the quasiparticle energies can be positive or negative depending on their index. In an on-shell scattering event, this delta function forces two of the indices to be positive and the remaining two to be negative. This has to be considered in interpreting the collision terms $`C_{PP}`$ and $`C_{\alpha P}`$, because there all the sums run over positive and negative indices. The principal-value part in Eq. (20), which appears in the single-particle kinetic equations, is absorbed in the $`T`$ matrices in Sec. III and thus does not appear anymore in the quasiparticle equations.

To complete the presentation of the quasiparticle kinetic equations, we write the generalized Gross-Pitaevskii equation (9) with the updated GP propagator $`\mathrm{\Pi }^{}`$ and the second-order collisional contributions expressed in terms of quasiparticle populations $`P`$ and obtain

$$\frac{d}{dt}\chi =i\mathrm{\Pi }^{}\chi +W\sigma _3\left(\mathrm{\Gamma }_{PP(\sigma _3+P)}\mathrm{\Gamma }_{(\sigma _3+P)(\sigma _3+P)P}\right)W^1\chi $$
(86)
In Sec. V, we write this equation for $`\alpha (t)`$ alone. This equation together with Eq. (80) gives a complete description in terms of Bogoliubov quasiparticles of a condensate coupled to a thermal cloud at finite temperatures.

### IV.3 Orders of magnitude

We now want to discuss the orders of magnitude of several of the quantities of this theory. This will suggest some approximations to the full quasiparticle kinetic equations (80).

We first consider the basis transformation $`W`$. The completeness relation for the quasiparticle basis Eq. (64) tells us that

$$W=𝒪(1).$$
(87)
For example, for high-energy eigenfunctions, the effect of the condensate becomes small, and the quasiparticle transformation reduces to

$$u1\text{and}v0,\text{as}E\mathrm{}.$$
(88)
Now, considering the basis-rotation term defined in Eq. (81) we find

$$B_W=W^1\frac{dW}{dt}=\frac{𝒪(1)}{dt}=𝒪(\mathrm{\Gamma })<𝒪\left(\frac{1}{\tau _c}\right),$$
(89)
because the time-scale for population changes, which change the quasiparticle basis $`W`$, is limited by the time between collisions $`dt>\tau _c`$. In equilibrium, the populations are constant due to detailed balance of the in and out rates. Thus, the net collision rate $`\mathrm{\Gamma }C_{PP}+C_{\alpha P}+\mathrm{H}.\mathrm{c}.`$ gives a better estimate for population changes $`dt\mathrm{\Gamma }^1`$. This also confirms that $`B_W=0`$ in equilibrium, since $`dW/dt=0`$.

We now show that the stationary solutions $`P`$ of the Boltzmann equation (80) are diagonal. Considering the stationary solution $`\frac{d}{dt}P_{\overline{1}\overline{2}}=0`$ of Eq. (80) for an off-diagonal element with $`\overline{1}\overline{2}`$, we obtain

$$P_{\overline{1}\overline{2}}=i\frac{B_WP+C_{PP}+C_{\alpha P}+\mathrm{H}.\mathrm{c}.}{E_{\overline{2}}E_{\overline{1}}}=𝒪\left(\frac{\mathrm{\Gamma }}{\mathrm{\Delta }ϵ}\right).$$
(90)
This shows that the off-diagonal elements of the quasiparticle population are small compared to the diagonal ones, which are of the order of the number of particles $`N`$, and vanish at equilibrium.

### IV.4 Quasiparticle $`T`$ matrices

The $`T`$ matrix in Eq. (39) (or the ladder extension Eq. (55)) has been obtained appropriately to second order. The second-order term with the factor of 1 in Eq. (39) is the divergent part, which is renormalized in Eq. (55) when replaced by the $`T`$-matrix. The term containing $`2\stackrel{~}{f}`$ of Eq. (39) is a convergent many-body second-order term. It is thus reasonable to replace the $`T`$ of Eq. (55) with

$$T^{1^{}2^{}3^{}4^{}}(E)=T_{2B}^{1^{}2^{}3^{}4^{}}(E+2\mu )+8𝒫\frac{\varphi ^{1^{}2^{}3^{\prime \prime }4^{\prime \prime }}\stackrel{~}{f}_{4^{\prime \prime }2^{\prime \prime }}\varphi ^{3^{\prime \prime }2^{\prime \prime }3^{}4^{}}}{E(E_{3^{\prime \prime }}+E_{4^{\prime \prime }})},$$
(91)
where we replaced the single-particle energies in the denominator by quasiparticle energies, because the difference is of higher order. We here use the following two-body $`T`$ matrix

$$T_{\text{2B}}^{1^{}2^{}3^{}4^{}}(ϵ)=2\varphi ^{1^{}2^{}3^{}4^{}}+4𝒫\frac{\varphi ^{1^{}2^{}3^{\prime \prime }4^{\prime \prime }}\varphi ^{3^{\prime \prime }4^{\prime \prime }3^{}4^{}}}{ϵ(ϵ_{3^{\prime \prime }}^{\text{sp}}+ϵ_{4^{\prime \prime }}^{\text{sp}})},$$
(92)
which is given in terms of single-particle energies $`ϵ^{\text{sp}}`$ defined by

$$\widehat{H}^{(0)}|ϵ_1^{}^{\text{sp}}=\left(\frac{\widehat{𝐩}^2}{2m}+V_{\mathrm{ext}}(\widehat{𝐱^{}})\right)|ϵ_1^{}^{\text{sp}}=ϵ_1^{}^{\text{sp}}|ϵ_1^{}^{\text{sp}}.$$
(93)
The collisional terms of the kinetic equations (80) correspond to the imaginary part of a Liouville-space $`𝒯`$-matrix (56). Examining the argument of the $`\delta _\eta `$-function in Eq. (85), which defines $`\mathrm{\Psi }_\eta `$, in comparison with the ladder $`T`$ in Eq. (55) shows that (when the energies are correlated with the appropriate elements of W, see Wachter2000a ) two of the energies in the denominator must be positive and two negative, which, together with the fact that $`P`$ is diagonal close to equilibrium, leads to terms in the kinetic equation of the form $`PP(1+P)(1+P)`$, etc. The kinetic equation (80) then becomes the quantum Boltzmann equation for quasiparticle populations. The equilibrium solution is therefore the expected Bose-Einstein distribution for the quasiparticles, as the steady-state solution of Eqs. (80) and (86) shows (Walser2001a, , Sec. V). The interaction matrix elements $`\varphi `$ in Eq. (77) can also be upgraded to $`T`$s. These equations in terms of $`T`$ are now consistent with an impact approximation treatment (with elastic scattering not contributing when $`TT^{}`$ terms are considered, cf Eq. (56)).

### IV.5 Conservation laws

Because of the non-vanishing pairing field, the trace of the quasiparticle density matrix $`P`$ is in general not equal to the number of excited particles in the system. Hence, the number of quasi-particles is not conserved. Our single-particle kinetic equations (9) and (23) do, however, conserve the mean total number

$$N=\alpha ^{}\alpha +\mathrm{Tr}\{\stackrel{~}{f}\}=N_c+\mathrm{Tr}\{\stackrel{~}{f}\}=\mathrm{const};$$
(94)
this can be proven explicitly by plugging these kinetic equations into $`\frac{d}{dt}N`$ and canceling terms. We thus adopt a self-consistent procedure for the quasiparticle kinetic equations Eqs. (80) and (86), by requiring

$$N=N_c+\mathrm{Tr}\{UPU^{}\}=\mathrm{const}.$$
(95)
This equation self-consistently constrains the number of condensate atoms $`N_c`$: in equilibrium at temperature $`\beta ^1`$, the quasiparticle matrix $`P`$ consists of $`q=0`$ and

$$p(E)=\frac{1}{e^{\beta (E\mu )}1}$$
(96)
according to Eq. (71) with the chemical potential $`\mu `$ given by the GP equation (52). As temperature tends to zero, we obtain the usual corrections for the anomalous average and condensate depletion Morgan2000a . Note that the number of excited atoms is not given by the trace of $`p`$, but has to include the basis transformation $`U`$ as discussed above. If we drop the basis-rotation term $`B_W`$ in numerical simulations, we will incur number-non-conservation on the order of $`\mathrm{\Gamma }`$, while away from equilibrium.

Since we use a Markovian collision integral with a damping function of finite width $`\eta `$ in order to include off-the-energy-shell propagation, this theory is not exactly energy conserving. Markovian theories fail to track the decay of initial correlations and thus do not account for the decay of the correlation energy Bonitz . In our case, with a self-consistent $`\eta `$, energy is conserved to order $`\eta `$, which is consistent with the order of approximation. For a detailed discussion of these memory effects and how they affect the conservation laws see Bhongale2002a .

## V Summary

We summarize the main results of this paper given in Secs. III and IV. We formulate a kinetic theory in terms of Bogoliubov quasiparticle modes $`W`$, which are defined by the eigenvalue equation for the renormalized Hartree-Fock-Bogoliubov propagator $`\mathrm{\Sigma }^{}`$

$$\left(\begin{array}{cc}\widehat{H}^{(0)}+2T_{f^c}^{}+2T_{\stackrel{~}{f}}^{}\mu & T_{m^c}\\ \mathrm{\Sigma }_𝒜^{}& \mathrm{\Sigma }_𝒩^{}\end{array}\right)W=WE.$$
(57)
The $`T`$ matrices are defined to second order in Eqs. (39) and (47) and Eq. (55) gives an extension to the ladder approximation. The Gross-Pitaevskii equation

$$\left\{\widehat{H}^{(0)}+T_{f^c}^{}(2\mu )+2T_{\stackrel{~}{f}}^{}(2\mu )\right\}\alpha =\mu \alpha $$
(52)
for the ground state $`\alpha `$ is contained in Eq. (57), because the renormalized $`\mathrm{\Sigma }^{}`$ is gapless. The GP equation defines the value of the chemical potential $`\mu `$. To find a complete basis $`W`$, we have to find the second zero-energy mode and form quadrature components as discussed in Sec. IV.1.

We find the following Boltzmann equation for the thermal excitations in terms of Bogoliubov quasiparticles:

$`{\displaystyle \frac{d}{dt}}P`$ $`=`$ $`i[E,P]+\{\sigma _3W^{}\sigma _3{\displaystyle \frac{dW}{dt}}P+\mathrm{H}.\mathrm{c}.\}`$
$`+`$ $`\{\mathrm{\Gamma }_{PP(\sigma _3+P)}(\sigma _3+P)\mathrm{\Gamma }_{(\sigma _3+P)(\sigma _3+P)P}P+\mathrm{H}.\mathrm{c}.\}`$
$`+`$ $`3\{\mathrm{\Gamma }_{P^cP(\sigma _3+P)}(\sigma _3+P)\mathrm{\Gamma }_{P^c(\sigma _3+P)P}P+\mathrm{H}.\mathrm{c}.\}.`$
The basis-rotation term containing $`\frac{d}{dt}W`$ can be dropped for adiabatic evolution. The quasiparticle density matrix $`P`$ is diagonal close to equilibrium, and its elements obey an equilibrium Bose-Einstein distribution as the second and third lines show. The collision terms containing the general condensate matrix $`P^c`$ defined in Eq. (72) represent population exchange between the thermal cloud and the condensate. They are balanced in the following equation for the condensate

$`{\displaystyle \frac{d}{dt}}\alpha (t)`$ $`=`$ $`\left\{\widehat{H}^{(0)}+T_{f^c}^{}+2T_{\stackrel{~}{f}}^{}\mu \right\}\alpha (t)`$
$`+`$ $`U\sigma _3\left\{\mathrm{\Gamma }_{PP(\sigma _3+P)}\mathrm{\Gamma }_{(\sigma _3+P)(\sigma _3+P)P}\right\}\sigma _3`$
$`\times \left\{U^{}\alpha (t)V^{}\alpha ^{}(t)\right\}.`$
In general, $`\alpha (t)`$ can be different from the $`\alpha `$ used as the ground state of the adiabatic basis. A non-adiabatic change in a driving force, for example, would cause $`\alpha (t)`$ to change quickly.

The coupled Eqs. (80) and (86) have to be solved under the orthogonality constraints (75) or (74) depending on whether the evolution is adiabatic or not. Furthermore, the total particle number has to fulfill the self-consistent constraint Eq. (95).

## VI Conclusions

We have extended the non-equilibrium kinetic theory of Walser et al. Walser1999a ; Walser2001a  in two important respects. First, we write the binary interactions as many-body $`T`$ matrices to second order in the interaction by subsuming ultra-violet divergent terms. This procedure removes divergences caused by the anomalous average $`\stackrel{~}{m}`$ and contained in the second-order terms. We can then replace the low-energy limit of the $`T`$ matrix with the $`s`$-wave scattering length. This is to our knowledge the only consistent treatment of these difficulties associated with the anomalous average in a theory that includes second-order collisional terms. The updated Hartree-Fock-Bogoliubov propagator $`\mathrm{\Sigma }^{}`$ is then gapless, which greatly facilitates a consistent treatment of the condensate coupled to thermal fluctuations.

The second extension of the Walser et al. theory makes use of the gapless HFB propagator by using its quasiparticle eigenmodes to parameterize the thermal fluctuations, which are then automatically orthogonal to the condensate mean field. We find that this basis greatly simplifies the second-order collision terms of the Walser et al. theory. Another important result reported here is that, in equilibrium, the Bogoliubov modes diagonalize the quasiparticle population matrix $`P`$. This means that, close to equilibrium, the second-order terms can be evaluated in $`n^4`$ operations, where $`n`$ is the number of energy levels considered. This is a vast improvement compared to $`n^8`$ operations for a basis that does not diagonalize the population matrix, such as, for example, the single-particle basis used by Walser et al.

## Acknowledgments

J.W. acknowledges financial support from the National Science Foundation. R.W. gratefully acknowledges support from the Austrian Academy of Sciences through an APART grant. M.H. acknowledges support from the U.S. Department of Energy, Office of Basic Energy Sciences via the Chemical Sciences, Geosciences, and Biosciences Division. J.C. is retired and acknowledges support from TIAA-CREF.


format: text

text:
# Ultrasonic Attenuation in the Vortex State of d-wave Superconductors

## Abstract

We calculate the low temperature quasi-particle contribution to the ultrasonic attenuation rate in the mixed state of d-wave superconductors. Our calculation is performed within the semi-classical approximation using quasi-particle energies that are Doppler shifted, with respect to their values in the Meissner phase, by the supercurrent associated with the vortices. We find that the attenuation at low temperatures and at fields $`H_{c1}HH_{c2}`$ has a temperature independent contribution which is proportional to $`\sqrt{H}`$ where $`H`$ is the applied magnetic field. We indicate how our result in combination with the zero-field result for ultrasonic attenuation can be used to calculate one of the parameters $`v_F`$, $`H_{c2}`$ or $`\xi `$ given the values for any two of them.

The discovery of high temperature superconductivity Bednorz  in the cuprates has led to an intense theoretical effort at understanding the origins of the novel phenomena seen in these materials. The normal state of the cuprates is highly anomalous PWA  and remains a puzzle that is still not understood. However, there is evidence Mohit  that well defined quasi-particles exist in the superconducting state. The superconducting state is now known Leggett ; Kirtley  to be a d-wave state whose order parameter has $`d_{x^2y^2}`$ symmetry. The order parameter thus has nodes along the lines $`k_x=\pm k_y`$ in the two-dimensional Brillouin zone.

While the origin of the pairing interaction that leads to the occurence of d-wave superconductivity is not known at the present time, a great deal of progress in understanding the superconducting state can be made by focussing on the low energy nodal quasi-particles near the points where the Fermi surface intersects the $`k_x=\pm k_y`$ lines. These quasi-particles, which have an energy dispersion corresponding to Dirac fermions Lee  in the Meissner phase, are the dominant low-energy excitations which determine the low temperature properties in the superconducting state.

It is possible to derive MPLB  the d-wave gap from a model Hamiltonian but the connection of that Hamiltonian with the underlying physics of strongly interacting fermions is unclear. Thus for the purposes of this paper, we assume on phenomenological grounds, the existence of BCS like quasi-particles with an energy gap that has d-wave symmetry. These quasi-particles might not be the real strongly interacting electrons, to whom their connection is unclear at the present time.

Ultrasonic attenuation has historically been a very useful tool in investigating the energies of quasi-particles in conventional superconductors. Early verification Morse  of the BCS prediction for the temperature dependence of the energy gap was done using this technique. Even afterwards Bishop ; Batlogg  this has been a useful tool in the investigation of of superconductivity in heavy fermion superconductors. On the theoretical front calculations of the ultrasonic attenuation in d-wave superconductors have been carried out both in the clean MPLB ; Klemm ; Vekh ; Carb  and dirty Coleman ; New  limit in the Meissner phase. The former limit which corresponds to $`Ql1`$ , Q being the ultrasound wave-vector and l being the electronic mean free path, is the one considered by us in this paper.

Assuming 100 MHz to be a typical frequency Bishop ; Batlogg  (although ultrasound experiments can be carried out over a large frequency range from kHz to GHz) at which ultrasound experiments are done and taking the sound velocity to beShobo  $`4\times 10^5`$m/s we find that this translates into the requirement $`l4\mu m`$. This restricts the applicability of our work to the cleanest samplesClean .

While ultrasonic attenuation in the vortex state has recieved attention in the past Pankert ; Faraday1 ; Faraday2 ; Blatter  the emphasis of earlier workers has been to study effects arising from the coupling of ultrasonic waves to vortex motion. This coupling comes from two different effects: a) The pinning potentials from impurity ions tend to drag the vortices along with them , b) The ionic current due to the ionic displacement will exert a Lorentz force on the vortices. In the clean limit($`l4\mu m`$) in which we work, the effects of pinning are expected to be weak. The Lorentz force coupling is in any case weak because it is suppressed by a factor of $`v/c`$ where $`v`$ is the ionic velocity and $`c`$ is the velocity of light. We therefore ignore these processes and focus entirely on the phonon damping from the electronic quasi-particles.

In conventional s-wave superconductors, low energy electronic excitations in the vortex state are the bound states localised in the core region of the vortices Caroli . However in d-wave superconductors, because of the presence of nodes in the gap function, the dominant low energy excitations are those in the far region away from the vortex cores. Their contribution is expected to be overwhelmingly larger Volovik  than that coming from the the cores. Further, experimental studies from STM measurements STM1 ; STM2  reveal that that there are just a few bound states in the vortex cores. Thus we focus exclusively on the excitations in the far region.

To describe the low energy excitations in the mixed phase we make the use of the semi-classical approximation, first discussed DeGennes  in the case of s-wave superconductors. This approximation has been employed for d-wave superconductors in recent times Volovik ; Hirschfeld  and its regime of applicability is the range $`H_{c1}HH_{c2}`$, where H is the applied magnetic field and $`H_{c1}`$ and $`H_{c2}`$ are the lower critical and upper critical magnetic fields repectively. A striking success of this method has been the prediction Volovik  of a term linear in the temperature T whose co-efficient scales as $`\sqrt{H}`$ (H being the applied magnetic field). This prediction has been verified spht  and has given greater credibility to the semi-classical description. Recently this approximation has been systematized and put on a firm footing by Ramakrishnan and RajagopalTVR  who have derived it microscopically.

Using the semi-classical approximation to describe the electronic Green functions in the superconducting state we calculate the imaginary part of the electron density-density correlation function which is proportional to the inverse phonon lifetime. As the sound velocity is only weakly dependent on temperatureShobo , the temperature dependence of the attenuation comes almost entirely from the the inverse lifetime. We find that at low tempeartures, and for parameters appropiate to the cuprates, the ultrasonic attenuation co-efficient in the mixed state has a temperature independent contribution which scales as $`\sqrt{H}`$. The co-efficient of this term and its dependence on the ultrasound wave-vector as well as the leading temperature corrections to it have been explicitly evaluated. These are our main results.

We now present the details of our calculation. The semi-classical approximation can be understood DeGennes as ignoring, in the first instance, the spatial variation of the supercurrents around a vortex outside its core region. This can be justified on the grounds that the spatial variation of the electronic wave-function has a characteristic length scale $`k_F^1`$ which is much smaller than the smallest length scale associated with vortex currents $`\xi `$ (the core radius). In the cuprates the parameter $`(k_F\xi )^110`$ and thus the semi-classical approximation is expected to provide a reasonable description. We then proceed to evaluate the inverse phonon lifetime to first order in $`\stackrel{}{q}`$ where a current has been introduced by taking the energy gap function to be $`\mathrm{\Delta }_\stackrel{}{q}(\stackrel{}{R_i},\stackrel{}{R_j})=\mathrm{\Delta }_{i,j}e^{i\stackrel{}{q}(\stackrel{}{R_i}+\stackrel{}{R_j})}`$. Here $`\stackrel{}{R_i}`$ and $`\stackrel{}{R_j}`$ are neighbouring sites on a square lattice whose fermions are paired up in a singlet state and $`\mathrm{\Delta }_{i,j}=\mathrm{\Delta }(\mathrm{\Delta })`$ for $`\stackrel{}{R_i}\stackrel{}{R_j}`$ being in the $`\widehat{x}(\widehat{y})`$ direction. The spatial dependence of $`\stackrel{}{q}=\widehat{\varphi }/2r`$ is now restored and the the inverse phonon lifetime is averaged over a unit-cell of the vortex lattice to get our final result for the attenuation co-efficient.

Solving the Bogulibov-deGennes equation for the mean-field d-wave superconductor to linear order in $`\stackrel{}{q}`$ we obtain the results:

$$G_\stackrel{}{q}^\sigma (\stackrel{}{R_i}\stackrel{}{R_j},i\omega _n)=\frac{1}{N}\underset{\stackrel{}{k}}{}e^{i\stackrel{}{k}(\stackrel{}{R_i}\stackrel{}{R_j})}[\frac{u_{\stackrel{}{k}\stackrel{}{q}}^2}{i\omega _nE_{\stackrel{}{k}\stackrel{}{q}}}+\frac{v_{\stackrel{}{q}\stackrel{}{k}}^2}{i\omega _n+E_{\stackrel{}{q}\stackrel{}{k}}}]$$
(1)
where $`u_{\stackrel{}{k}}^{}{}_{}{}^{2}=(1+ϵ_\stackrel{}{k}/E_\stackrel{}{k}^0)/2,v_{\stackrel{}{k}}^{}{}_{}{}^{2}=(1ϵ_\stackrel{}{k}/E_\stackrel{}{k}^0)/2,ϵ_\stackrel{}{k}=\xi _\stackrel{}{k}\mu ,\xi _\stackrel{}{k}`$ being the band energy and $`\mu `$ the chemical potential, $`E_\stackrel{}{k}^0=\sqrt{ϵ_\stackrel{}{k}^2+\mathrm{\Delta }_\stackrel{}{k}^2},\mathrm{\Delta }_\stackrel{}{k}=\mathrm{\Delta }(\mathrm{cos}k_xa\mathrm{cos}k_ya)`$ being the d-wave gap and $`E_\stackrel{}{k}=E_\stackrel{}{k}^0+\stackrel{}{q}_k\xi _\stackrel{}{k}`$ being the Doppler shifted quasi-particle energy. Here

$$G_\stackrel{}{q}^\sigma (\stackrel{}{R_i}\stackrel{}{R_j},i\omega _n)=_0^\beta 𝑑\tau <c_{i,\sigma }(\tau )c_{j,\sigma }^{}(0)>e^{i\omega _n\tau }$$
(2)
is the ”normal” Green function. Similarly the anomalous Green functions are given by:

$$F_\stackrel{}{q}(\stackrel{}{R_i},\stackrel{}{R_j},i\omega _n)=\frac{e^{i\stackrel{}{q}(\stackrel{}{R_i}+\stackrel{}{R_j})}}{N}\underset{\stackrel{}{k}}{}u_\stackrel{}{k}v_\stackrel{}{k}(\frac{e^{i\stackrel{}{k}(\stackrel{}{R_i}\stackrel{}{R_j})}}{i\omega _nE_\stackrel{}{k}}\frac{e^{i\stackrel{}{k}(\stackrel{}{R_i}\stackrel{}{R_j})}}{i\omega _n+E_\stackrel{}{k}})$$
(3)
and

$$F_\stackrel{}{q}^+(\stackrel{}{R_i},\stackrel{}{R_j},i\omega _n)=\frac{e^{i\stackrel{}{q}(\stackrel{}{R_i}+\stackrel{}{R_j})}}{N}\underset{\stackrel{}{k}}{}u_\stackrel{}{k}v_\stackrel{}{k}(\frac{e^{i\stackrel{}{k}(\stackrel{}{R_i}\stackrel{}{R_j})}}{i\omega _nE_\stackrel{}{k}}\frac{e^{i\stackrel{}{k}(\stackrel{}{R_i}\stackrel{}{R_j})}}{i\omega _n+E_\stackrel{}{k}})$$
(4)
Here $`u_\stackrel{}{k}v_\stackrel{}{k}=\mathrm{\Delta }_\stackrel{}{k}/2E_\stackrel{}{k}^0`$ and the Green functions $`F_\stackrel{}{q}`$ and $`F_\stackrel{}{q}^+`$ are defined as

$$F_\stackrel{}{q}(\stackrel{}{R_i},\stackrel{}{R_j},i\omega _n)=_0^\beta 𝑑\tau <c_{i,}(\tau )c_{j,}(0)>e^{i\omega _n\tau }$$
(5)
$$F_\stackrel{}{q}^+(\stackrel{}{R_i},\stackrel{}{R_j},i\omega _n)=_0^\beta 𝑑\tau <c_{i,}^{}(\tau )c_{j,}^{}(0)>e^{i\omega _n\tau }$$
(6)
Using the Green functions in equations (1),(3) and (4) we compute the imaginary part of the density-density correlation function which is given by

$$\chi _\stackrel{}{q}^\mathrm{"}(\stackrel{}{Q},\omega )=X_1+X_2+X_3+X_4$$
(7)
where

$$X1=\frac{2\pi }{N}\underset{\stackrel{}{k}}{}[n(E_\stackrel{}{k})n(E_{\stackrel{}{k}+\stackrel{}{Q}})](u_{\stackrel{}{k}+\stackrel{}{Q}}^2u_\stackrel{}{k}^2\frac{\mathrm{\Delta }_\stackrel{}{k}\mathrm{\Delta }_{\stackrel{}{k}+\stackrel{}{Q}}}{4E_\stackrel{}{k}^0E_{\stackrel{}{k}+\stackrel{}{Q}}^0})\delta (\omega +E_\stackrel{}{k}E_{\stackrel{}{k}+\stackrel{}{Q}})$$
(8)
$$X2=\frac{2\pi }{N}\underset{\stackrel{}{k}}{}[n(E_\stackrel{}{k})n(E_{\stackrel{}{k}\stackrel{}{Q}})](v_{\stackrel{}{k}+\stackrel{}{Q}}^2v_\stackrel{}{k}^2\frac{\mathrm{\Delta }_\stackrel{}{k}\mathrm{\Delta }_{\stackrel{}{k}+\stackrel{}{Q}}}{4E_\stackrel{}{k}^0E_{\stackrel{}{k}+\stackrel{}{Q}}^0})\delta (\omega E_\stackrel{}{k}+E_{\stackrel{}{k}\stackrel{}{Q}})$$
(9)
$$X3=\frac{2\pi }{N}\underset{\stackrel{}{k}}{}[1n(E_\stackrel{}{k})n(E_{\stackrel{}{k}+\stackrel{}{Q}})](u_{\stackrel{}{k}+\stackrel{}{Q}}^2v_\stackrel{}{k}^2+\frac{\mathrm{\Delta }_\stackrel{}{k}\mathrm{\Delta }_{\stackrel{}{k}+\stackrel{}{Q}}}{4E_\stackrel{}{k}^0E_{\stackrel{}{k}+\stackrel{}{Q}}^0})\delta (\omega E_\stackrel{}{k}E_{\stackrel{}{k}+\stackrel{}{Q}}),$$
(10)
and

$$X4=\frac{2\pi }{N}\underset{\stackrel{}{k}}{}[1n(E_\stackrel{}{k})n(E_{\stackrel{}{k}\stackrel{}{Q}})](v_{\stackrel{}{k}+\stackrel{}{Q}}^2u_\stackrel{}{k}^2+\frac{\mathrm{\Delta }_\stackrel{}{k}\mathrm{\Delta }_{\stackrel{}{k}+\stackrel{}{Q}}}{4E_\stackrel{}{k}^0E_{\stackrel{}{k}+\stackrel{}{Q}}^0})\delta (\omega +E_\stackrel{}{k}+E_{\stackrel{}{k}\stackrel{}{Q}}).$$
(11)
In s-wave superconductors the contribution of $`X_3`$ and $`X_4`$ are zero as the ultrasound frequency $`\omega 2\mathrm{\Delta }`$ and so the $`\delta `$-fn condition in Eqs.(10) and (11) can never be satisfied. In the d-wave case the presence of nodes in $`\mathrm{\Delta }_\stackrel{}{k}`$ means that $`X_3`$ and $`X_4`$ will be finite. However, as the phase space for them is limited, their contribution is expected to be small and so we focus exclusively on $`X_1`$ and $`X_2`$ which make the dominant contribution to $`\chi _\stackrel{}{q}^\mathrm{"}(\stackrel{}{Q},\omega )`$.

Now expanding to the leading order in $`\omega `$ we find

$`\chi _\stackrel{}{q}^\mathrm{"}(\stackrel{}{Q},\omega )={\displaystyle \frac{2\pi \omega }{N}}{\displaystyle \underset{\stackrel{}{k}}{}}[n^{}(E_\stackrel{}{k})\delta (E_\stackrel{}{k}E_{\stackrel{}{k}+\stackrel{}{Q}})(u_{\stackrel{}{k}+\stackrel{}{Q}}^2u_\stackrel{}{k}^2{\displaystyle \frac{\mathrm{\Delta }_\stackrel{}{k}\mathrm{\Delta }_{\stackrel{}{k}+\stackrel{}{Q}}}{4E_\stackrel{}{k}^0E_{\stackrel{}{k}+\stackrel{}{Q}}^0}})+n^{}(E_\stackrel{}{k})`$
$`\delta (E_\stackrel{}{k}E_{\stackrel{}{k}\stackrel{}{Q}})(v_{\stackrel{}{k}+\stackrel{}{Q}}^2v_\stackrel{}{k}^2{\displaystyle \frac{\mathrm{\Delta }_\stackrel{}{k}\mathrm{\Delta }_{\stackrel{}{k}+\stackrel{}{Q}}}{4E_\stackrel{}{k}^0E_{\stackrel{}{k}+\stackrel{}{Q}}^0}})]`$ (12)
We next expand $`\chi ^\mathrm{"}`$ to leading order in the ultrasound wave-vector $`\stackrel{}{Q}`$ and on making the substitution $`\stackrel{}{k}\stackrel{}{k}`$ in the second term in Eq. (12) we arrive at the result

$$\chi _\stackrel{}{q}^\mathrm{"}(\stackrel{}{Q},\omega )=\frac{2\pi \omega }{N}\mathrm{\Sigma }_\stackrel{}{k}n^{}(E_\stackrel{}{k}^0+\stackrel{}{q}\stackrel{}{_k}\xi _\stackrel{}{k})\frac{ϵ_\stackrel{}{k}^2}{E_\stackrel{}{k}^{0^2}}\delta (\stackrel{}{Q}\stackrel{}{_k}E_\stackrel{}{k}^0+q_\alpha Q_\beta \frac{^2\xi _\stackrel{}{k}}{k_\alpha k_\beta })$$
(13)
where a sum over repeated indices is implied in the arguement of the $`\delta `$ function.

To do the $`\stackrel{}{k}`$ summation in Eq. (13) we only consider Lee  the fermions near the nodes. This approximation is reasonable at low temperatures $`(k_BT\mathrm{\Delta })`$. For the node near $`\stackrel{}{P_a}=(\pi /2a,\pi /2a)`$ upon linearising the band energies and gap function, we have: $`\xi _\stackrel{}{k}v_Fk_1`$ and $`\mathrm{\Delta }_\stackrel{}{k}v_\mathrm{\Delta }k_2`$ where $`k_1`$ and $`k_2`$ are co-ordinates normal and tangential respectively to the Fermi surface at the node. In terms of the co-ordinates $`k_1`$ and $`k_2`$ we find that the contribution from this node to $`\chi _\stackrel{}{q}^{\mathrm{"}(a)}`$ can be written as

$`\chi _\stackrel{}{q}^{\mathrm{"}(a)}(\stackrel{}{Q},\omega ){\displaystyle \frac{a^2\omega }{2\pi }}{\displaystyle 𝑑k_1𝑑k_2n^{}[\sqrt{v_F^2k_1^2+v_\mathrm{\Delta }^2k_2^2}+\rho _0]\frac{v_F^2k_1^2}{v_F^2k_1^2+v_\mathrm{\Delta }^2k_2^2}}`$
$`\delta [{\displaystyle \frac{\alpha v_Fk_1+\beta v_\mathrm{\Delta }k_2}{\sqrt{v_F^2k_1^2+v_\mathrm{\Delta }^2k_2^2}}}\alpha _1v_Fk_1\alpha _2v_Fk_2]`$ (14)
Here $`\alpha =v_FQ\mathrm{cos}(\pi /4\theta )`$, $`\beta =v_\mathrm{\Delta }Q\mathrm{sin}(\pi /4\theta )`$, $`\alpha _1=(Qqa^2/2)\mathrm{cos}(\psi \theta )`$, $`\alpha _2=(Qqa^2/2)\mathrm{cos}(\psi +\theta )`$, $`\rho _0=qv_F\mathrm{cos}(\psi \pi /4)`$ and $`\theta `$ and $`\psi `$ are the angles made by $`\stackrel{}{Q}`$ and $`\stackrel{}{q}`$ respectively with respect to the $`k_1`$ axis. We now introduce the polar co-ordinates $`v_Fk_1`$ = $`\rho \mathrm{cos}\varphi `$ and $`v_\mathrm{\Delta }k_2`$ = $`\rho \mathrm{sin}\varphi `$. Then upon performing the $`\varphi `$ integral we find $`\chi _\stackrel{}{q}^{\mathrm{"}(a)}`$ reduces to

$$\chi _\stackrel{}{q}^{\mathrm{"}(a)}(\stackrel{}{Q},\omega )=\frac{a^2\omega }{\pi v_Fv_\mathrm{\Delta }}_0^{\rho _c}𝑑\rho \rho n^{}[\rho +\rho _0]f_a[\rho ]$$
(15)
where $`f_a[\rho ]=\frac{B^2}{(A^2+B^2)^{3/2}}`$ and $`A=\alpha \alpha _1\rho `$ , $`B=\beta (v_F/v_\mathrm{\Delta })\alpha _2\rho `$ and $`\rho _c=\sqrt{\pi v_Fv_\mathrm{\Delta }}/a`$ is a cutoff introduced to preserve the volume of the Brillouin zone while doing the $`\stackrel{}{k}`$-integration.

Introducing the variable $`\rho ^{}=\rho +\rho _0`$ we have $`\chi ^{\mathrm{"}a}=\chi ^{\mathrm{"}(a1)}+\chi ^{\mathrm{"}(a2)}`$ where

$$\chi _\stackrel{}{q}^{\mathrm{"}(a1)}(\stackrel{}{Q},\omega )=\frac{a^2\omega }{\pi v_Fv_\mathrm{\Delta }}_{\rho _0}^{\rho _c}𝑑\rho ^{}\rho ^{}n^{}[\rho ^{}]f_a[\rho ^{}\rho _0]$$
(16)
and

$$\chi _\stackrel{}{q}^{\mathrm{"}(a2)}(\stackrel{}{Q},\omega )=\frac{a^2\omega }{\pi v_Fv_\mathrm{\Delta }}\rho _0_{\rho _0}^{\rho _c}𝑑\rho ^{}n^{}[\rho ^{}]f_a[\rho ^{}\rho _0]$$
(17)
In writing Eqs. (16) and (17) we have made use of $`\rho _c\rho _0`$ which follows from using $`q_{max}=1/2\xi `$ and $`\xi 30\AA `$, $`v_F/a192meV`$ and $`v_\mathrm{\Delta }/a28meV`$ for parameters appropriate to the cupratesLee .

In an analogous fashion we evaluate the corresponding contributions to $`\chi _\stackrel{}{q}^\mathrm{"}`$ from the nodes near $`\stackrel{}{P_b}=(\pi /2a,\pi /2a)`$, $`\stackrel{}{P_c}=(\pi /2a,\pi /2a)`$ and $`\stackrel{}{P_d}=(\pi /2a,\pi /2a)`$. We first focus on the type of terms in Eq. (16). In this case we find, on expanding to linear order in $`q`$, that the contribution linear in $`q`$ exactly vanishes due to cancellations from the terms coming from different nodes. We thus find $`\chi _\stackrel{}{q}^{\mathrm{"}(1)}=\chi _\stackrel{}{q}^{\mathrm{"}(a1)}+\chi _\stackrel{}{q}^{\mathrm{"}(b1)}+\chi _\stackrel{}{q}^{\mathrm{"}(c1)}+\chi _\stackrel{}{q}^{\mathrm{"}(d1)}`$ to be given by

$$\chi _\stackrel{}{q}^{\mathrm{"}(1)}=\frac{2ln2\omega a^2k_BT}{\pi v_Fv_\mathrm{\Delta }}[\frac{\beta ^2}{(\alpha ^2+\beta ^2)^{3/2}}+\frac{\eta ^2}{(\lambda ^2+\eta ^2)^{3/2}}]$$
(18)
to linear order in $`\stackrel{}{q}`$. Once again we have used $`k_BTqv_F\rho _c`$. Here $`\alpha `$ and $`\beta `$ are previously defined and we have introduced the parameters $`\lambda =v_FQ\mathrm{sin}(\pi /4\theta )`$, $`\eta =v_\mathrm{\Delta }Q\mathrm{cos}(\pi /4\theta )`$. The result obtained for this term is identical to the zero field result for the attenuation previously calculated in Ref. (12). We next turn our attention to the terms of the type written in Eq.(17). Here on adding the contributions from the nodes near $`\stackrel{}{P}_a`$ and $`\stackrel{}{P}_b`$, we find to leading order in $`q`$,

$$\chi _\stackrel{}{q}^{\mathrm{"}(a2)}+\chi _\stackrel{}{q}^{\mathrm{"}(b2)}=\frac{\omega a^2q}{\pi v_\mathrm{\Delta }}\mathrm{cos}(\psi \pi /4)\frac{\beta ^2}{(\alpha ^2+\beta ^2)^{3/2}}\mathrm{tanh}(\frac{qv_F\mathrm{cos}(\psi \pi /4)}{2k_BT})$$
(19)
Similarly the nodes near $`\stackrel{}{P}_c`$ and $`\stackrel{}{P}_d`$ yield the contribution

$$\chi _\stackrel{}{q}^{\mathrm{"}(c2)}+\chi _\stackrel{}{q}^{\mathrm{"}(d2)}=(\frac{\omega a^2q}{\pi v_\mathrm{\Delta }})\mathrm{cos}(\psi +\pi /4)\frac{\eta ^2}{(\lambda ^2+\eta ^2)^{3/2}}\mathrm{tanh}(\frac{qv_F\mathrm{cos}(\psi +\pi /4)}{2k_BT})$$
(20)
It is straightforward to see that for $`T0`$, Eqs. (19) and (20) reduce to

$$\chi _\stackrel{}{q}^{\mathrm{"}(a2)}+\chi _\stackrel{}{q}^{\mathrm{"}(b2)}=(\frac{\omega a^2q}{\pi v_\mathrm{\Delta }})\mathrm{cos}(\psi \pi /4)\frac{\beta ^2}{(\alpha ^2+\beta ^2)^{3/2}}[\theta (\mathrm{cos}(\psi \pi /4))\theta (\mathrm{cos}(\psi \pi /4))]$$
(21)
and

$$\chi _\stackrel{}{q}^{\mathrm{"}(c2)}+\chi _\stackrel{}{q}^{\mathrm{"}(d2)}=(\frac{\omega a^2q}{\pi v_\mathrm{\Delta }})\mathrm{cos}(\psi +\pi /4)\frac{\eta ^2}{(\lambda ^2+\eta ^2)^{3/2}}[\theta (\mathrm{cos}(\psi +\pi /4))\theta (\mathrm{cos}(\psi +\pi /4))]$$
(22)
Eqs.(21) and (22) are also obtained directly by taking the T = 0 limit for n’ in Eq.(15) and its analogues for the other 3 nodes.

To proceed further we now restore $`\stackrel{}{q}=\widehat{\varphi }/2r`$ corresponding to the supercurrents around a vortex where the vector potential has been ignored as the Ginzburg-Landau parameter $`\kappa 1`$. Assuming a vortex lattice with circular unit-cells we find that the radius of the cells $`R_c\xi \sqrt{H_{c2}/H}`$. We then average $`\chi _\stackrel{}{q}^\mathrm{"}`$ over one unit cell to obtain

$$\chi ^\mathrm{"}(H,T)=\frac{_\xi ^{R_c}r𝑑r_0^{2\pi }𝑑\psi \chi _\stackrel{}{q}^\mathrm{"}}{_\xi ^{R_c}r𝑑r_0^{2\pi }𝑑\psi }$$
(23)
to obtain our final result for $`\chi ^\mathrm{"}`$.

In order to perform the averaging in Eq. (23) we make the approximation $`tanh(x)x`$ for $`|x|<1`$, $`tanh(x)1`$ for $`x>1`$ and $`tanh(x)1`$ for $`x<1`$. This approximation interpolates between the asymptotically exact behaviour at small and large $`x`$. Further it reproduces the exact result that can be directly calculated for $`T=0`$. In that case the integral in Eq. (23) is elementary and can be done by substituting the expressions in Eqs. (21) and (22). We then obtain to leading order in the small parameter $`(\xi k_BT/v_F)^2(H_{c2}/H)`$:

$$\chi ^{\mathrm{"}(2)}(H,T)=\frac{2\omega a^2}{\pi ^2v_\mathrm{\Delta }\xi }(\frac{H}{H_{c2}})^{1/2}[1\frac{8\xi ^2}{9}(\frac{k_BT}{v_F})^2\frac{H_{c2}}{H}][\frac{\beta ^2}{(\alpha ^2+\beta ^2)^{3/2}}+\frac{\eta ^2}{(\lambda ^2+\eta ^2)^{3/2}}]$$
(24)
The assumption about the small parameter implies that the regime of validity of our result is for $`k_BT\frac{v_F}{\xi }(\frac{H}{H_{c2}})^{1/2}`$ together with the condition $`H_{c1}HH_{c2}`$. For parameters relevant to the cuprates $`H_{c2}/H_{c1}100`$ at low temperatures. This restricts the temperature window to $`T1K`$ at $`H=H_{c1}`$. Our result for a temperature independent contribution to the attenuation that scales as $`\sqrt{H}`$ is entirely understandable as it has its origin in a finite density of states at the Fermi-energy Volovik  whose size is proportional to $`\sqrt{H}`$ and whose signature is seen in the specific heat measurementsspht .

The actual ultrasonic attenuation coefficient is related to $`\chi ^\mathrm{"}`$ by

$$\alpha _S(T,H)=M(\stackrel{}{Q})\chi ^\mathrm{"}(\stackrel{}{Q},T,H)$$
(25)
where $`M(\stackrel{}{Q})`$ is a constant which depends on the sound velocity and the electron-phonon matrix element and the ultrasound frequency. Then on combining Eqs. (18) and (24) and considering the fact that Eq. (18) contains the result for the Meissner phase ($`H=0`$) we obtain the result:

$$\frac{\alpha _S(T,H)\alpha _S(T,H=0)}{\alpha _S(T,H=0)}=\frac{v_F}{\pi \xi }(\frac{H}{H_{c2}})^{1/2}\frac{1}{ln2k_BT}[1\frac{8\xi ^2}{9}(\frac{k_BT}{v_F})^2\frac{H_{c2}}{H}]$$
(26)
whose window of validity has been described above. This result is remarkable because it is independent of the ultrasound wave-vector $`\stackrel{}{Q}`$ as well as the gap increase parameter $`v_\mathrm{\Delta }`$. Thus a measurement of the ultrasonic attenuation at low temperatures and at fields above $`H_{c1}`$ would enable a determination of any one of the parameters $`v_F`$,$`\xi `$ and $`H_{c2}`$ given a knowledge of the other two.

We now discuss some shortcomings of our work. We have assumed perfectly well-defined quasi-particles thus ignoring the incoherent spectral weight seen in photoemission experiments Mohit . The semi-classical approximation used by us may not accurately describe all the physical effects due to the scattering of quasi-particles from the supercurrents in the vortex state. Our results are restricted to the clean limit and are applicable only to very clean samples.

Finally we conclude by recapitulating the main points of this paper. We have employed the semi-classical approximation to evaluate the phonon damping due to the electronic quasi-particles in a d-wave superconductor. We find that for parameters appropiate to the cuprates, in a temperature window $`k_BT1K`$, there is a temperature independent contribution to the ultrasonic attenuation whose magnitude scales as $`\sqrt{H}`$.

###### Acknowledgements.

One of us (DMG) thanks T.V.Ramakrishnan for useful comments.