text:
# Nonequilibrium relaxation in neutral BCS superconductors: Ginzburg-Landau approach with Landau damping in real time

## I Introduction

Nonequilibrium phenomena in superconductors continue to be the focus of attention. The dynamics of Josephson junctions, phase slip phenomena in the dynamics of vortices and relaxation of the order parameter and supercurrents are few examples of the experimental effort that probe nonequilibrium aspects of superconductivity.

Since the original work of Abrahams and Tsuneto tsuneto  there has been an ongoing effort in trying to obtain the effective time-dependent description of nonequilibrium phenomena from a microscopic Bardeen-Cooper-Schrieffer (BCS) BCS ; deGennes  theory. Whereas the effective Ginzburg-Landau (GL) GL  description in the *static* limit was derived by Gor’kov gorkov59 ; book:fetter , the effective time dependent Ginzburg-Landau description is still the focus of a substantial theoretical effort. There is a large body of work that established the validity of a time-dependent nonlinear Schroedinger equation that describes the dynamics of the order parameter at zero temperature fraser ; wilczek ; schakel ; aitchison95 ; stone ; palo .

At *finite* temperature the dynamical description is complicated by the presence of Landau damping which prevents a local description in time because the spectral densities feature branch cuts that prevent a derivative expansion. This problem was originally pointed out by Abrahams and Tsuneto tsuneto . At finite temperature Landau damping cuts are unavoidable and result from processes that involve scatterings of quasiparticles in the thermal bath. In derivations of the effective Lagrangian for dynamical phenomena from a microscopic theory the Landau damping contribution had often been ignored stoof .

The absorptive contributions to the effective action of long wavelength phase fluctuations at finite temperature have been studied by Aitchison et al. aitchison97  for a neutral BCS superconductor. These authors studied in detail the Landau damping contributions to the effective action of phase fluctuations and concluded that for temperatures $`0<T<0.6T_c`$ the effective propagator for the phase fluctuation (Bogoliubov-Anderson-Goldstone mode) can be well approximated by simple quasiparticle poles at complex energy and describe damped excitations with a linear and temperature dependent dispersion relation and narrow widths.

An alternative approach to study nonequilibrium aspects of superconductors is based on kinetic theory. Kopnin kopnin  studied the nonequilibrium dynamics of flux flow in clean superconductors but did not address the validity of the time dependent Landau-Ginzburg description near the critical point. Watts-Tobin et al. watts  studied the validity of the Landau-Ginzburg description near the critical point for *dirty* superconductors where relaxational processes are dominated by (elastic) collisions.

However, to the best of our knowledge, the description of the relaxational dynamics for the amplitude and the phase of the order parameter, as well as the validity of the Landau-Ginzburg description *near the critical point* in *clean* neutral superconductors had been the subject of several recent studies but has not been completely understood. The region near the critical temperature, where $`|\mathrm{\Delta }_0(T)|T`$ with $`\mathrm{\Delta }_0(T)`$ the finite-temperature gap (order parameter), is the region of validity of the Ginzburg-Landau theory.

The interest on a deeper understanding of the time dependent effective action for long-wavelength phase fluctuations has been rekindled by several recent developments. Recently there has been a substantial effort to obtain the time dependent effective action of long-wavelength collective excitations associated with phase fluctuations in $`d`$-wave superconductors sharapov ; takada ; para ; benfatto .

In particular these studies focused on the novel Carlson-Goldman modes carlson , which are Goldstone-like modes in *charged* superconductors that emerge near the critical temperature. While at $`T=0`$ the Anderson-Higgs mechanism combines the Goldstone and gauge fields into a gapped plasma mode, near the critical temperature a novel quasiparticle Goldstone-like excitation, the Carlson-Goldman mode, is present in charged superconductors. This mode is a superposition of the Bogoliubov-Anderson-Goldstone mode, present in *neutral* superconductors and the long-range gauge field which is screened at finite temperature. In Ref. takada,  it was pointed out that the existence of this mode is associated both with screening and Landau damping of phase fluctuations. The importance of the nonequilibrium dynamics of long-wavelength *phase* fluctuations has also been highlighted recently within the context of high-temperature superconductivity HTSC .

Furthermore, recent experiments in ultracold alkali atoms have demonstrated the trapping and cooling of *fermionic* alkalis, in particular $`{}_{}{}^{40}\mathrm{K}`$ and $`{}_{}{}^{6}\mathrm{Li}`$ schreck . One goal of this present experimental effort is to observe a transition to a *neutral fermi superfluid* for fermi systems with an *attractive* interaction between atoms in two different hyperfine states stoof2 . Recently the spectrum of low energy collective excitations in the *collisionless* regime has been studied, in particular focusing on the emergence of Goldstone or phase fluctuations in these *neutral Fermi superfluids* mottelson . A proposal for the detection of the phase transition to a neutral Fermi superfluid in $`{}_{}{}^{6}\mathrm{Li}`$ alkalis relies on the spectrum of long-wavelength collective (Goldstone) excitations zambelli .

The interest on neutral BCS Fermi superfluids is interdisciplinary, from the current experimental efforts in Fermi alkalis and in mixtures with Bose alkalis schreck , to neutron superfluidity in nuclear matter and neutron stars. For a recent discussion on neutral Fermi superfluids and their interest in a wide variety of fields see Ref. pethick, . Hence the study of the dynamics of phase fluctuations in neutral Fermi superfluids (or neutral BCS superconductors) continues to be of timely interest and of experimental relevance.

The goals of this article: In this article we focus on the nonequilibrium *real-time* dynamics of phase and amplitude fluctuations in neutral BCS superconductors in the Ginzburg-Landau regime near the critical temperature. In particular we obtain the effective *dynamical* Ginzburg-Landau description of nonequilibrium relaxation of long-wavelength, low-frequency fluctuations of the order parameter *near the critical point*.

While previous efforts, notably by Aitchison et al. aitchison95 ; aitchison97 , focused on the long-wavelength, low-frequency effective action well below the critical temperature for $`0<T<0.6T_c`$ our goal is to study the critical region $`|\mathrm{\Delta }_0(T)|T_c`$ with $`\mathrm{\Delta }_0(T)`$ the finite-temperature gap. Our study is different from previous attempts in several respects: (i) We implement the Schwinger-Keldysh formulation of nonequilibrium field theory schwinger  along with the recently introduced *tadpole method* tadpole  to obtain the equations of motion for small amplitude fluctuations of the order parameter in *real time*. (ii) The equations of motion obtained with these methods are retarded, lead to the Ward identities and allow to establish the retarded effective action at once. Furthermore, the equations of motion describe an initial value problem that allows a real-time study of relaxation and damping. (iii) We then focus on the Ginzburg-Landau regime $`|\mathrm{\Delta }_0(T)|/T1`$ and establish the *dynamical* Ginzburg-Landau effective action to one-loop order for long-wavelength, low-frequency fluctuations. This effective action is retarded and nonlocal because of Landau damping. (iv) We study the time evolution of small phase and amplitude fluctuations from the equilibrium configuration and reveal directly in real time the effect of Landau damping.

Main results: Implementing the Schwinger-Keldysh formulation of nonequilibrium field theory and the tadpole method we obtain the retarded equations of motion for small fluctuations, which in turn lead to Ward identities associated with (global) gauge invariance both *in* and *out of* equilibrium. From these equation of motion we obtain the retarded one-loop effective action which is *nonlocal* as a consequence of Landau damping.

We then focus on the Ginzburg-Landau ($`\mathrm{\Delta }_0T`$) region and study the real-time relaxation of small phase and amplitude fluctuations. While the spectral density for phase fluctuations features a peak that suggests a Goldstone-like dispersion relation, the relaxational dynamics is completely *overdamped* as a consequence of Landau damping.

Far away from the Ginzburg-Landau regime at low temperatures, the spectral densities for both phase and amplitude fluctuations feature narrow quasiparticle peaks confirming previous results aitchison95 ; aitchison97 . In particular, the real-time relaxation of long-wavelength phase fluctuations is weakly *underdamped* by Landau damping.

The article is organized as follows. In Sec. II we introduce the model and the linear response formulation to obtain the equations of motion. In Sec. III we introduce the Schwinger-Keldysh formulation in the Nambu-Gor’kov formalism to study the nonequilibrium aspects of Bogoliubov quasiparticles. In Sec. IV we introduce the tadpole method, obtain the equations of motion directly in real time and cast them in terms of an initial value problem. We obtain explicitly the retarded self-energies to one-loop order and their spectral representations and obtain the one-loop retarded effective action. In Sec. V we obtain the Ward identities and discuss the static limit of the self-energies. In Sec. VI we obtain the effective time dependent Ginzburg-Landau description focusing on the Ginzburg-Landau regime and the long-wavelength, low-frequency limit. In this section we provide a thorough numerical analysis of the real-time evolution of the relaxation of phase and amplitude fluctuations. Section VII presents our conclusions and poses new directions. An appendix is devoted to an alternative derivation of the Bogoliubov transformation, which facilitates the Schwinger-Keldysh nonequilibrium formulation.

## II Preliminaries

### II.1 Neutral BCS model

The BCS Hamiltonian of a neutral electron gas is given by

$`H`$ $`=`$ $`{\displaystyle \underset{\sigma =,}{}}{\displaystyle d^3x\psi _\sigma ^{}(𝐱,t)\left(\frac{^2}{2m}\right)\psi _\sigma (𝐱,t)}g{\displaystyle d^3x\psi _{}^{}(𝐱,t)\psi _{}^{}(𝐱,t)\psi _{}(𝐱,t)\psi _{}(𝐱,t)},`$ (1)
where $`\psi _\sigma (𝐱,t)`$ are the Heisenberg complex fields representing electrons of mass $`m`$ and spin $`\sigma `$, and $`g>0`$ is the strength of the $`s`$-wave pairing interaction between spin-up and spin-down electrons close to the Fermi surface. In this article, we set $`\mathrm{}=k_\mathrm{B}=1`$. The field $`\psi _\sigma (𝐱,t)`$ and its Hermitian conjugate satisfy the equal-time *anticommutation* relations

$$\{\psi _\sigma (𝐱,t),\psi _\sigma ^{}^{}(𝐱{}_{}{}^{},t)\}=\delta _{\sigma \sigma ^{}}\delta ^{(3)}(𝐱𝐱{}_{}{}^{}),$$

$$\{\psi _\sigma (𝐱,t),\psi _\sigma ^{}(𝐱{}_{}{}^{},t)\}=\{\psi _\sigma ^{}(𝐱,t),\psi _\sigma ^{}^{}(𝐱{}_{}{}^{},t)\}=0.$$
(2)
The Hamiltonian $`H`$ is invariant under the $`U(1)`$ gauge transformation

$$\psi _\sigma (𝐱,t)e^{i\theta }\psi _\sigma (𝐱,t),$$

$$\psi _\sigma ^{}(𝐱,t)e^{i\theta }\psi _\sigma ^{}(𝐱,t),$$
(3)
where $`\theta `$ is a constant phase. A consequence of this $`U(1)`$ gauge symmetry is conservation of the number of electrons. Indeed, the number operator of electrons

$$N=\underset{\sigma =,}{}d^3x\psi _\sigma ^{}(𝐱,t)\psi _\sigma (𝐱,t)$$
(4)
commutes with $`H`$ and hence is a constant of motion. However, it is convenient to work in the grand-canonical ensemble in which the grand-canonical Hamiltonian is given by

$`K`$ $``$ $`H\mu N`$ (5)
$`=`$ $`{\displaystyle \underset{\sigma =,}{}}{\displaystyle d^3x\psi _\sigma ^{}(𝐱,t)\left(\frac{^2}{2m}\mu \right)\psi _\sigma (𝐱,t)}g{\displaystyle d^3x\psi _{}^{}(𝐱,t)\psi _{}^{}(𝐱,t)\psi _{}(𝐱,t)\psi _{}(𝐱,t)},`$
where the chemical potential $`\mu `$ is the Lagrange multiplier associated with conservation of number of electrons. The chemical potential $`\mu `$ is determined by fixing the number of electrons and in general is a function of the temperature. However, for the situation under study in which the temperature is much lower than the Fermi temperature, $`\mu `$ can be approximated by its zero-temperature value, i.e., the Fermi energy. The Lagrangian (density) corresponding to $`K`$ is given by

$$[\psi ^{},\psi ]=\underset{\sigma =,}{}\psi _\sigma ^{}\left(i\frac{}{t}+\frac{^2}{2m}+\mu \right)\psi _\sigma +g\psi _{}^{}\psi _{}^{}\psi _{}\psi _{}.$$
(6)
Introducing the auxiliary complex scalar *pair field* and its Hermitian conjugate defined as

$$\mathrm{\Delta }(𝐱,t)=g\psi _{}(𝐱,t)\psi _{}(𝐱,t),$$

$$\mathrm{\Delta }^{}(𝐱,t)=g\psi _{}^{}(𝐱,t)\psi _{}^{}(𝐱,t),$$
(7)
and performing the Hubbard-Stratonovich transformationnegele , the Lagrangian can be written as

$$[\psi ^{},\psi ,\mathrm{\Delta }^{},\mathrm{\Delta }]=\underset{\sigma =,}{}\psi _\sigma ^{}\left(i\frac{}{t}+\frac{^2}{2m}+\mu \right)\psi _\sigma +\mathrm{\Delta }^{}\psi _{}\psi _{}+\psi _{}^{}\psi _{}^{}\mathrm{\Delta }\frac{1}{g}\mathrm{\Delta }^{}\mathrm{\Delta }.$$
(8)
We note that the pair field $`\mathrm{\Delta }(𝐱,t)`$ is *not* a dynamical field as there is no corresponding kinetic term in the Lagrangian (8).

In the superconducting phase, we decompose the paring field into the condensate and noncondensate parts

$$\mathrm{\Delta }(𝐱,t)=\mathrm{\Delta }(𝐱,t)+\chi (𝐱,t),\chi (𝐱,t)=0,$$
(9)
where $`𝒪(𝐱,t)=\mathrm{Tr}[\rho 𝒪(𝐱,t)]/\mathrm{Tr}\rho `$denotes the *expectation value* of the Heisenberg operator $`𝒪(𝐱,t)`$ in the *initial density matrix* $`\rho `$, $`\mathrm{\Delta }(𝐱,t)`$ is the superconducting order parameter, and $`\chi (𝐱,t)`$ describes the noncondensate operator. The presence of the condensate $`\mathrm{\Delta }(𝐱,t)0`$ leads to spontaneous breaking of the $`U(1)`$ gauge symmetry. In the absence of explicit symmetry breaking external sources, the condensate is homogeneous (i.e., space-time independent) $`\mathrm{\Delta }(𝐱,t)=\mathrm{\Delta }_0`$, which is the situation under consideration in this article.

### II.2 Real-time relaxation in linear response

The goal of this article is to obtain *directly in real time* the equations of motion for small amplitude perturbations of the homogeneous superconducting condensate in an initial value problem formulation. Our strategy to study the relaxation of the condensate perturbation as an initial value problem begins with preparing a superconducting state slightly perturbed away from equilibrium by applying an external source coupled to the pair field. Once the external source is switched off, the perturbed condensate must relax towards equilibrium. It is precisely this *real-time evolution* of the nonequilibrium fluctuations around the condensate the focus of this article.

Let $`\eta (𝐱,t)`$ be an external $`c`$-number source coupled to the pair field $`\mathrm{\Delta }(𝐱,t)`$, then the Lagrangian given by (8) becomes

$$[\psi ^{},\psi ,\mathrm{\Delta }^{},\mathrm{\Delta }][\psi ^{},\psi ,\mathrm{\Delta }^{},\mathrm{\Delta }]+\mathrm{\Delta }^{}\eta +\eta ^{}\mathrm{\Delta }.$$
(10)
The presence of the external source $`\eta `$ will induce a (linear) response of the system in the form of an induced expectation value

$$\mathrm{\Delta }(𝐱,t)_\eta =\mathrm{\Delta }_0+\delta (𝐱,t).$$
(11)
Here, $`\mathrm{\Delta }(𝐱,t)_\eta `$ denotes the expectation value of the paring field $`\mathrm{\Delta }(𝐱,t)`$ in the presence of the external source, $`\mathrm{\Delta }_0=\mathrm{\Delta }(𝐱,t)`$ is the homogeneous order parameter in the absence the external source, and $`\delta (𝐱,t)`$ is the space-time dependent perturbation of the homogeneous condensate $`\mathrm{\Delta }_0`$ *induced by the external source*. The linear response perturbation $`\delta (𝐱,t)`$ vanishes when the external source $`\eta (𝐱,t)`$ vanishes at all times. This is tantamount to decomposing the field into the homogeneous condensate ($`\mathrm{\Delta }_0`$), a small amplitude perturbation induced by the external source \[$`\delta (𝐱,t)`$\], and the noncondensate part \[$`\chi (𝐱,t)`$\] as

$$\mathrm{\Delta }(𝐱,t)=\mathrm{\Delta }_0+\delta (𝐱,t)+\chi (𝐱,t),\chi (𝐱,t)_\eta =0.$$
(12)
In linear response theory $`\delta (𝐱,t)`$ can be expressed in terms of the *exact* retarded Green’s function of the pair fields in the absence of external source book:fetter ; ivp . An experimentally relevant initial value problem formulation for the real-time relaxation of the condensate perturbation can be obtained by considering that the external source is adiabatically switched on at $`t=\mathrm{}`$ and switched off at $`t=0`$, i.e.,

$$\eta (𝐱,t)=\eta (𝐱)e^{ϵt}\mathrm{\Theta }(t),ϵ0^+.$$
(13)
The adiabatic switching-on of the external source induces a space-time dependent condensate perturbation $`\delta (𝐱,t)`$, which is prepared adiabatically by the external source with a given value $`\delta (𝐱,0)`$ at $`t=0`$ determined by $`\eta (𝐱)`$. For $`t>0`$ after the external source has been switched off, the perturbed condensate will evolve in the absence of any external source relaxing towards equilibrium. Thus, the external source $`\eta (𝐱,t)`$ is necessary for preparing an initial state at $`t=0`$ setting up an initial value problem. This method has been applied to study a wide variety of relaxation phenomena in different settings ivp ; mikheev , including the relaxation of condensate fluctuations in homogeneous Bose-Einstein condensates boyanbec .

Using the decomposition (12) we expand the Lagrangian density, and consistently with linear response, keep only the linear terms in $`\delta `$ and $`\delta ^{}`$, which are the small amplitude perturbations from the homogeneous condensate induced by the external source $`\eta `$. The Lagrangian $``$ becomes (in the presence of the external source $`\eta `$)

$$[\psi ^{},\psi ,\chi ^{},\chi ]=_0[\psi ^{},\psi ,\chi ^{},\chi ]+_{\mathrm{int}}[\psi ^{},\psi ,\chi ^{},\chi ],$$
(14)
with

$`_0[\psi ^{},\psi ,\chi ^{},\chi ]`$ $`=`$ $`{\displaystyle \underset{\sigma =,}{}}\psi _\sigma ^{}\left(i{\displaystyle \frac{}{t}}+{\displaystyle \frac{^2}{2m}}+\mu \right)\psi _\sigma +\mathrm{\Delta }_0^{}\psi _{}\psi _{}+\psi _{}^{}\psi _{}^{}\mathrm{\Delta }_0{\displaystyle \frac{1}{g}}\chi ^{}\chi ,`$
$`_{\mathrm{int}}[\psi ^{},\psi ,\chi ^{},\chi ]`$ $`=`$ $`(\delta ^{}+\chi ^{})\psi _{}\psi _{}{\displaystyle \frac{1}{g}}(\delta ^{}+\mathrm{\Delta }_0^{})\chi +\eta ^{}\chi +\mathrm{H}.\mathrm{c}.,`$ (15)
where we have discarded the $`c`$-number (field operators independent) terms. We note that the Lagrangian $`[\psi ^{},\psi ,\chi ^{},\chi ]`$ is obviously invariant under the gauge transformations

$`\psi _\sigma ,\chi ,\mathrm{\Delta }_0,\delta ,\eta e^{i\theta /2}\psi _\sigma ,e^{i\theta }\chi ,e^{i\theta }\mathrm{\Delta }_0,e^{i\theta }\delta ,e^{i\theta }\eta ,`$
$`\psi _\sigma ^{},\chi ^{},\mathrm{\Delta }_0^{},\delta ^{},\eta ^{}e^{i\theta /2}\psi _\sigma ^{},e^{i\theta }\chi ^{},e^{i\theta }\mathrm{\Delta }_0,e^{i\theta }\delta ,e^{i\theta }\eta ,`$ (16)
which, as will be seen below, is at the heart of the Ward identity.

Whereas in general a gauge transformation is invoked to fix the condensate $`\mathrm{\Delta }_0`$ to be *real* for convenience, this choice corresponds to *fixing a particular gauge*, which in turn hides the underlying gauge symmetry. In order to obtain the Ward identity associated with this symmetry we will keep a complex condensate $`\mathrm{\Delta }_0`$ and analyze in detail the transformation laws of the various contributions to the equations of motion.

The study of the static and dynamical properties of the BCS theory is simplified by introducing the Nambu-Gor’kov formulation. Let us introduce the two-component Nambu-Gor’kov fields nambu ; book:fetter

$$\mathrm{\Psi }(𝐱,t)=\left[\begin{array}{c}\psi _{}(𝐱,t)\\ \psi _{}^{}(𝐱,t)\end{array}\right],\mathrm{\Psi }^{}(𝐱,t)=[\psi _{}^{}(𝐱,t),\psi _{}(𝐱,t)],$$
(17)
and the $`2\times 2`$ Pauli matrices

$$\sigma _+=\left[\begin{array}{cc}0& 1\\ 0& 0\end{array}\right],\sigma _{}=\left[\begin{array}{cc}0& 0\\ 1& 0\end{array}\right],\sigma _3=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right],$$
(18)
in terms of which the Lagrangian can be written as

$$[\mathrm{\Psi }^{},\mathrm{\Psi },\chi ^{},\chi ]=_0[\mathrm{\Psi }^{},\mathrm{\Psi },\chi ^{},\chi ]+_{\mathrm{int}}[\mathrm{\Psi }^{},\mathrm{\Psi },\chi ^{},\chi ],$$
(19)
with

$`_0[\mathrm{\Psi }^{},\mathrm{\Psi },\chi ^{},\chi ]`$ $`=`$ $`\mathrm{\Psi }^{}\left[i{\displaystyle \frac{}{t}}+\sigma _3\left({\displaystyle \frac{^2}{2m}}+\mu \right)+\sigma _+\mathrm{\Delta }_0+\sigma _{}\mathrm{\Delta }_0^{}\right]\mathrm{\Psi }{\displaystyle \frac{1}{g}}\chi ^{}\chi ,`$
$`_{\mathrm{int}}[\mathrm{\Psi }^{},\mathrm{\Psi },\chi ^{},\chi ]`$ $`=`$ $`\mathrm{\Psi }^{}\sigma _+(\delta +\chi )\mathrm{\Psi }+{\displaystyle \frac{1}{g}}\chi ^{}(\mathrm{\Delta }_0+\delta )+\chi ^{}\eta +\mathrm{H}.\mathrm{c}..`$ (20)
## III Nonequilibrium formulation

### III.1 Generating functional

The general framework to study of nonequilibrium phenomena is the Schwinger-Keldysh formulation schwinger , which we briefly review here in a manner that leads immediately to a path integral formulation.

Consider that the system is described by an initial density matrix $`\rho `$ and a perturbation is switched on at a time $`t_0`$, so that the total Hamiltonian for $`t>t_0`$, $`H(t)`$, does not commute with the initial density matrix. The expectation value of a Heisenberg operator $`𝒪(t)=U^1(t,t_0)𝒪(t_0)U(t,t_0)`$ is given by

$$𝒪(t)=\frac{\mathrm{Tr}\rho U^1(t,t_0)𝒪U(t,t_0)}{\mathrm{Tr}\rho },$$
(21)
where $`U(t,t_0)`$ is the unitary time evolution operator in the Heisenberg picture

$$U(t,t_0)=T\mathrm{exp}\left[i_{t_0}^t𝑑t^{}H(t^{})\right],$$
(22)
with $`T`$ the time-ordering symbol. If the initial density matrix $`\rho `$ describes a state in thermal equilibrium at inverse temperature $`\beta `$ with the unperturbed Hamiltonian $`H(t<t_0)=H`$, i.e.,

$$\rho =e^{\beta H}=U(t_0i\beta ,t_0),$$
(23)
then the expectation value (21) can be written in the form

$$𝒪(t)=\frac{\mathrm{Tr}U(t_0i\beta ,t_0)U^1(t_0,t)𝒪U(t,t_0)}{\mathrm{Tr}U(t_0i\beta ,t_0)}.$$
(24)
The numerator of this expression has the following interpretation: evolve in time from $`t_0`$ up to $`t`$, insert the operator $`𝒪`$, evolve back from $`t`$ to the initial time $`t_0`$ and down the imaginary axis in time from $`t_0`$ to $`t_0i\beta `$. The denominator describes the evolution in imaginary time which is the familiar description of a thermal density matrix. We note that unlike the $`S`$-matrix elements or transition amplitudes, expectation values of Heisenberg operators require evolution *forward* and *backward* in time (corresponding to the $`U`$ and $`U^1`$ on each side of the operator $`𝒪`$).

The time evolution operators have a path-integral representation negele  in terms of the Lagrangian, and the insertion of operators can be systematically handled by introducing sources coupled linearly to the fields, i.e.,

$$[\mathrm{\Psi }^{},\mathrm{\Psi },\chi ^{},\chi ][\mathrm{\Psi }^{},\mathrm{\Psi },\chi ^{},\chi ]+\mathrm{\Psi }^{}J+J^{}\mathrm{\Psi }+\chi ^{}j+j^{}\chi ,$$
(25)
where $`J`$ and $`J^{}`$ are Grassmann-valued variables. The introduction of sources $`J`$, $`J^{}`$, $`j`$, and $`j^{}`$ also allows a systematic perturbative expansion. In such an expansion, powers of operators are obtained by functional derivatives with respect to these sources, which are set to zero after functional differentiation. We note that the sources $`j`$, $`j^{}`$ introduced in (25) to generate the perturbative expansion for the pair fields in terms of functional derivatives with respect to these, are *different* from the external sources $`\eta `$, $`\eta ^{}`$ introduced in (10) to generate an initial value problem in linear response and to displace the condensate from equilibrium.

Since there are *three* different time evolution operators, the forward, backward and imaginary, we introduce *three* different sources for each one of these time evolution operators, respectively. Taking $`t_0\mathrm{}`$, we are led to considering the generating functional comment

$$Z[J^+,J^{},J^\beta ]=\mathrm{Tr}U(\mathrm{}i\beta ,\mathrm{},J^\beta )U(\mathrm{},+\mathrm{},J^{})U(+\mathrm{},\mathrm{};J^+),$$
(26)
where $`U(t_f,t_i;J)`$ is the time evolution operator \[see (22)\] in the presence of the source $`J`$ and for simplicity of notation we have suppressed the spin index $`\sigma `$ and the complex conjugate of the sources $`J^{}`$. The denominator in (24) is given by $`\mathrm{Tr}\rho =Z[0,0,0]`$. The generating functional $`Z[J^+,J^{},J^\beta ]`$ can be written as a path integral along the contour in complex time plane (see Fig. 1)

$$Z[J^+,J^{},J^\beta ]=𝒟_𝒞\mathrm{\Psi }^{}𝒟_𝒞\mathrm{\Psi }\mathrm{exp}\left[i_𝒞d^4x_𝒞[\mathrm{\Psi }^{},\mathrm{\Psi };J]\right],$$
(27)
where $`𝒟_𝒞\mathrm{\Psi }^{}𝒟_𝒞\mathrm{\Psi }`$ denotes the functional integration measure along the contour $`𝒞`$ and

$`{\displaystyle _𝒞}d^4x_𝒞[\mathrm{\Psi }^{},\mathrm{\Psi };J]`$ $`=`$ $`{\displaystyle _{\mathrm{}}^+\mathrm{}}d^4x[\mathrm{\Psi }^+,\mathrm{\Psi }^+;J^+]{\displaystyle _{\mathrm{}}^+\mathrm{}}d^4x[\mathrm{\Psi }^{},\mathrm{\Psi }^{};J^{}]`$ (28)
$`+{\displaystyle _{\mathrm{}}^{\mathrm{}i\beta }}d^4x[\mathrm{\Psi }^\beta ,\mathrm{\Psi }^\beta ;J^\beta ],`$
with $`_{\mathrm{}}^+\mathrm{}d^4xd^3x_{\mathrm{}}^+\mathrm{}𝑑t`$, etc. Because of the trace and the fermionic nature of the operators, the path integral along the contour $`𝒞`$ requires *antiperiodic boundary conditions* on the fields. The superscripts $`+`$ and $``$ refer to fields defined in the upper and lower branches, respectively, corresponding to forward ($`+`$) and backward ($``$) time evolution, while the superscript $`\beta `$ refers to the field defined in the vertical branch running down parallel to the imaginary axis. The negative sign in front of the action along the backward branch is a result of the fact that backward time evolution is determined by $`U^1(+\mathrm{},\mathrm{})`$ with $`U`$ the time evolution operator. The contour source $`J`$ that enters in the contour Lagrangian $`_𝒞`$ in (27) takes the values of the sources $`J^\pm `$ and $`J^\beta `$ in the respective branches as displayed in Fig. 1.

Functional derivatives with respect to the sources in the forward branch give time-ordered Green’s functions, those with respect to the sources in the backward branch give the anti-time-ordered Green’s functions, and those with respect to the sources in the imaginary branch give the usual imaginary-time (Matsubara) Green’s functions. While the sources $`J^+`$, $`J^{}`$, and $`J^\beta `$ introduced to obtain the Green’s functions via functional differentiation are *different* in the different branches, as they generate the time-ordered, anti-time-ordered and Matsubara Green’s functions, respectively; the external source $`\eta `$, the homogeneous condensate $`\mathrm{\Delta }_0`$, and the departure from equilibrium $`\delta `$ are $`c`$-numbers and hence are the *same* in all branches.

Writing the Lagrangian as a free and an interaction part as $`=_0+_{\mathrm{int}}`$, the generating functional can be written as a power series expansion in the interaction part, which in turn can be generated by taking functional derivatives with respect to the sources $`J`$, $`J^{}`$ by identifying

$$\mathrm{\Psi }^\pm i\frac{\delta }{\delta J^\pm },\mathrm{\Psi }^\pm \pm i\frac{\delta }{\delta J^\pm },$$

$$\mathrm{\Psi }^\beta i\frac{\delta }{\delta J^\beta },\mathrm{\Psi }^\beta i\frac{\delta }{\delta J^\beta }.$$
(29)
As a result, the full generating functional along the contour $`𝒞`$ can be written as

$$Z[J]=\mathrm{exp}\left\{i_𝒞d^4x_{\mathrm{int},𝒞}[i\frac{\delta }{\delta J^{}},i\frac{\delta }{\delta J}]\right\}Z_0[J],$$
(30)
where free field generating functional $`Z_0[J]`$ is given by (27) and (28) but with $`[\mathrm{\Psi }^{},\mathrm{\Psi };J]`$ replaced by $`_0[\mathrm{\Psi }^{},\mathrm{\Psi };J]`$.

### III.2 Green’s functions

The equation of motion for the *free* Nambu-Gor’kov field $`\mathrm{\Psi }`$ in the presence of the source $`J`$ reads

$$\left[i\frac{}{t}+\sigma _3\left(\frac{^2}{2m}+\mu \right)+\sigma _+\mathrm{\Delta }_0+\sigma _{}\mathrm{\Delta }_0^{}\right]\mathrm{\Psi }(𝐱,t)=J(𝐱,t).$$
(31)
The solution of this equation of motion is given by

$$\mathrm{\Psi }_J(𝐱,t)=_𝒞d^4x^{}S(𝐱𝐱{}_{}{}^{},tt^{})J(𝐱{}_{}{}^{},t^{}),$$
(32)
where $`S(𝐱𝐱{}_{}{}^{},tt^{})`$ is the Green’s function matrix along the contour $`𝒞`$ and satisfies

$$[i\frac{}{t}+\sigma _3(\frac{^2}{2m}+\mu )+\sigma _+\mathrm{\Delta }_0+\sigma _{}\mathrm{\Delta }_0^{}]S(𝐱𝐱{}_{}{}^{},tt^{})=\delta _𝒞(tt^{})\delta ^{(3)}(𝐱𝐱{}_{}{}^{}),$$
(33)
with $`\delta _𝒞(tt^{})`$ the Dirac delta function along the contour $`_𝒞𝑑t^{}\delta _𝒞(tt^{})=1`$. The Green’s function $`S(𝐱𝐱{}_{}{}^{},tt^{})`$ has the form

$$S(𝐱𝐱{}_{}{}^{},tt^{})=S^>(𝐱𝐱{}_{}{}^{},tt^{})\mathrm{\Theta }_𝒞(tt^{})+S^<(𝐱𝐱{}_{}{}^{},tt^{})\mathrm{\Theta }_𝒞(t^{}t),$$
(34)
where $`\mathrm{\Theta }_𝒞(tt^{})`$ is the step function along the contour and $`S^{}(𝐱𝐱{}_{}{}^{},tt^{})`$ obey the homogeneous equations of motion.

The antiperiodic boundary conditions on the fields in the path integral, a result of the trace over fermionic fields in (26), lead to the following boundary condition on the Green’s function

$$\underset{t_0\mathrm{}}{lim}S(𝐱𝐱{}_{}{}^{},t_0t^{})=\underset{t_0\mathrm{}}{lim}S(𝐱𝐱{}_{}{}^{},t_0i\beta t^{}).$$
(35)
Since along the contour $`t_0\mathrm{}`$ is the *earliest* time and $`t_0i\beta `$ is therefore the *latest* time, (35) entails

$$\underset{t_0\mathrm{}}{lim}S^<(𝐱𝐱{}_{}{}^{},t_0t^{})=\underset{t_0\mathrm{}}{lim}S^>(𝐱𝐱{}_{}{}^{},t_0i\beta t^{}),t^{},$$
(36)
which is the Kubo-Martin-Schwinger (KMS) condition for equilibrium correlation functions book:kadanoff .

The free field generating functional $`Z_0[J]`$ is now obtained by writing

$$\mathrm{\Psi }(𝐱,t)=\stackrel{~}{\mathrm{\Psi }}(𝐱,t)+\mathrm{\Psi }_J(𝐱,t),$$

$$\mathrm{\Psi }^{}(𝐱,t)=\stackrel{~}{\mathrm{\Psi }}^{}(𝐱,t)+\mathrm{\Psi }_J^{}(𝐱,t),$$
(37)
which leads to the result

$$Z_0[J]=Z_0[0]\mathrm{exp}\left[i_𝒞d^4x_𝒞d^4x^{}J^{}(x)S(xx^{})J(x^{})\right],$$
(38)
where and hereafter $`x`$ denotes the space-time coordinates $`(𝐱,t)`$ for simplicity of notation. The source independent term $`Z_0[0]`$ will cancel between the numerator and the denominator in all expectation values in (21).

Furthermore, we are interested in computing Green’s functions of *finite real times* which are defined for fields in the forward ($`+`$) and backward ($``$) time branches but not in the imaginary branch. For these real-time Green’s functions the contributions to the generating functional from one source in the imaginary branch and another source in either the forward or backward branch vanish by the Riemann-Lebesgue lemma tadpole ; ivp , since the time arguments are infinitely far apart along the contour. Therefore the contour integrals of the source terms and Green’s functions in the generating functional factorize into a term in which the sources are those either in the forward and backward branches and another term in which *both* sources are in the imaginary branch tadpole ; ivp . The latter term (with both sources in the imaginary branch) cancel between numerator and denominator in expectation values and the only remnant of the imaginary branch is through the periodic boundary conditions along the full contour in the Green’s function.

Thus the generating functional for real-time Green’s functions simplifies to the following expression, defined solely along the forward and backward time branches tadpole ; ivp

$`Z[J^\pm ,J^\pm ]`$ $`=`$ $`\mathrm{exp}\left[i{\displaystyle _{\mathrm{}}^+\mathrm{}}d^4x\left\{_{\mathrm{int}}[i\delta /\delta J^+,i\delta /\delta J^+]_{\mathrm{int}}[i\delta /\delta J^{},i\delta /\delta J^{}]\right\}\right]`$ (39)
$`\times \mathrm{exp}\{i{\displaystyle _{\mathrm{}}^+\mathrm{}}d^4x{\displaystyle _{\mathrm{}}^+\mathrm{}}d^4x^{}[J^+(x)S^{++}(xx^{})J^+(x^{})+J^{}(x)`$
$`\times S^{}(xx^{})J^{}(x^{})J^+(x)S^+(xx^{})J^{}(x^{})J^{}(x)S^+(xx^{})J^+(x^{})]\},`$
with

$`S^{++}(xx^{})=S^>(𝐱𝐱{}_{}{}^{},tt^{})\mathrm{\Theta }(tt^{})+S^<(𝐱𝐱{}_{}{}^{},tt^{})\mathrm{\Theta }(t^{}t),`$
$`S^{}(xx^{})=S^>(𝐱𝐱{}_{}{}^{},tt^{})\mathrm{\Theta }(t^{}t)+S^<(𝐱𝐱{}_{}{}^{},tt^{})\mathrm{\Theta }(tt^{}),`$
$`S^+(xx^{})=S^>(𝐱𝐱{}_{}{}^{},tt^{}),`$
$`S^+(xx^{})=S^<(𝐱𝐱{}_{}{}^{},tt^{}),`$ (40)
where now $`\mathrm{}t,t^{}+\mathrm{}`$ and the superscripts $`+`$, $``$ correspond to the sources defined on the forward ($`+`$) and backward ($``$) time branches, respectively. An important issue that must be highlighted at this stage, is that derivatives with respect to sources in the forward ($`+`$) time branch correspond to insertion of operators *pre-multiplying* the density matrix $`\rho `$ and derivatives with respect to sources in the backward ($``$) branch correspond to the insertion of operators *post-multiplying* the density matrix. That this is so is a consequence of the fact that the density matrix evolves in time as $`U(t,t_0)\rho _0U^1(t,t_0)`$ with $`U(t,t_0)`$ the time evolution operator.

These four Green’s functions are not independent because of the identity

$$S^{++}(xx^{})+S^{}(xx^{})S^+(xx^{})S^+(xx^{})=0.$$
(41)
The diagonal elements in $`S^{++}(xx^{})`$ are the normal Green’s functions, representing the propagation of single electrons, whereas the off-diagonal elements are the anomalous Green’s functions, corresponding to the annihilation and creation of two electrons of opposite spins, respectively.

The functions $`S^{}(xx^{})`$, which are solutions of the homogeneous free field equation of motion, are simply related to the correlation functions of the free Nambu-Gor’kov fields $`\mathrm{\Psi }`$, $`\mathrm{\Psi }^{}`$. Indeed, taking variational derivatives of the free field generating functional $`Z_0[J]`$ with respect to $`J^\pm `$, $`J^\pm `$, one can show that

$$S_{ab}^>(xx^{})=i\mathrm{\Psi }_a(x)\mathrm{\Psi }_b^{}(x^{}),$$

$$S_{ab}^<(xx^{})=i\mathrm{\Psi }_b^{}(x^{})\mathrm{\Psi }_a(x),$$
(42)
where and hereafter $`a,b=1,2`$ denote the Nambu-Gor’kov indices. The expectation values in the expressions above are in the non-interacting thermal density matrix which corresponds to the quadratic part of the Lagrangian (Hamiltonian), i.e., the density matrix that describes free Bogoliubov quasiparticles in thermal equilibrium at inverse temperature $`\beta `$.

While the matrix elements $`S_{ab}^{}(xx^{})`$ can be obtained through the usual Bogoliubov transformation to the quasiparticle basis book:fetter , we present in the Appendix an alternative derivation of these correlation functions directly from the spinor solutions of the homogeneous equations of motion. We find that the correlation functions in (42) in the infinite volume limit are given by

$$S_{ab}^{}(xx^{})=\frac{d^3k}{(2\pi )^3}S_{ab}^{}(k,tt^{})e^{i𝐤(𝐱𝐱{}_{}{}^{})},$$
(43)
where

$$S_{ab}^>(k,tt^{})=i\left[[1n_\mathrm{F}(E_k)]𝒮_{ab}(k)e^{iE_k(tt^{})}+n_\mathrm{F}(E_k)\overline{𝒮}_{ab}(k)e^{iE_k(tt^{})}\right],$$

$$S_{ab}^<(k,tt^{})=i\left[n_\mathrm{F}(E_k)𝒮_{ab}(k)e^{iE_k(tt^{})}+[1n_\mathrm{F}(E_k)]\overline{𝒮}_{ab}(k)e^{iE_k(tt^{})}\right],$$
(44)
with

$$𝒮(k)=\left[\begin{array}{cc}u_k^2& u_kv_k^{}\\ u_kv_k& |v_k|^2\end{array}\right],\overline{𝒮}(k)=\left[\begin{array}{cc}|v_k|^2& u_kv_k^{}\\ u_kv_k& u_k^2\end{array}\right].$$
(45)
In the above expressions $`k|𝐤|`$, $`n_\mathrm{F}(\omega )=1/(e^{\beta \omega }+1)`$ is the Fermi-Dirac distribution function, $`u_k`$ and $`v_k`$ satisfying $`u_k^2+|v_k|^2=1`$ are the Bogoliubov coefficients \[see (145)\], and $`E_k`$ is the energy of free Bogoliubov quasiparticles (bogolons)

$$E_k=\sqrt{\xi _k^2+|\mathrm{\Delta }_0|^2},$$
(46)
where $`\xi _k=k^2/2m\mu `$ (see the Appendix). Using the relation $`1n_\mathrm{F}(\omega )=e^{\beta \omega }n_\mathrm{F}(\omega )`$, one can easily verify the KMS condition book:kadanoff

$$S^>(k,ti\beta t^{})=S^<(k,tt^{}).$$
(47)
Hence, the correlation functions for the Nambu-Gor’kov fields $`\mathrm{\Psi }`$, $`\mathrm{\Psi }^{}`$ that will enter in the nonequilibrium perturbative expansion are completely determined by (42)-(46).

### III.3 Feynman rules

From the generating functional of nonequilibrium Green’s functions (39), it is clear that the effective interaction Lagrangian relevant for the nonequilibrium calculations is given by

$$_{\mathrm{int}}^{\mathrm{eff}}[\mathrm{\Psi }^\pm ,\mathrm{\Psi }^\pm ,\chi ^\pm ,\chi ^\pm ]=_{\mathrm{int}}[\mathrm{\Psi }^+,\mathrm{\Psi }^+,\chi ^+,\chi ^+]_{\mathrm{int}}[\mathrm{\Psi }^{},\mathrm{\Psi }^{},\chi ^{},\chi ^{}],$$
(48)
where the fields $`\mathrm{\Psi }^\pm `$, $`\mathrm{\Psi }^\pm `$, $`\chi ^\pm `$ and $`\chi ^\pm `$ are defined on the forward ($`+`$) and backwards ($``$) time branches respectively. Consequently, this generating functional leads to the following Feynman rules that define the perturbative expansion for calculations of nonequilibrium expectation values.

* There are *two* sets of interaction vertices defined by $`_{\mathrm{int}}^{\mathrm{eff}}[\mathrm{\Psi }^\pm ,\mathrm{\Psi }^\pm ,\chi ^\pm ,\chi ^\pm ]`$: those in which the fields are in the forward ($`+`$) branch and those in which the fields are in the backward ($``$) branch. There is a relative minus sign between these two types of vertices.
* For each kind of fields there are *four* types of Green’s functions corresponding to correlations between fields defined in the forward or backward branch.

The Green’s functions of the Nambu-Gor’kov fields $`\mathrm{\Psi }`$, $`\mathrm{\Psi }^{}`$ are given by (40) in terms of the correlation functions displayed in (42) that are completely determined by (44), (45) and (46).

The Green’s functions of the pair fields $`\chi `$, $`\chi ^{}`$ can be obtained in an analogous manner. However, due to the nonpropagating nature of the pair fields, the Green’s functions are local in time and hence those of fields defined in different branches vanish identically.
* The combinatoric factors, trace over the Nambu-Gor’kov indices for fermion loops, etc., are the same as in the imaginary-time (Matsubara) formulation.

## IV Relaxation of condensate perturbations: An initial value problem

### IV.1 Equations of motion

The equations of motion for the small amplitude superconducting condensate perturbation $`\delta (x)`$ induced by the external source $`\eta (x)`$ is obtained by implementing the *tadpole method* tadpole ; ivp . This method begins by writing the pair fields $`\mathrm{\Delta }^\pm (x)`$ in the forward ($`+`$) and backward ($``$) time branches as

$$\mathrm{\Delta }^\pm (x)=\mathrm{\Delta }_0+\delta (x)+\chi ^\pm (x),$$
(49)
where $`\mathrm{\Delta }_0`$ is the homogeneous condensate in the absence of external source, $`\delta (x)`$ is the perturbation of the condensate induced by the external source which vanishes in the absence of external source, and $`\chi ^\pm (x)`$ are the noncondensate part of the pair fields in the forward and backward time branches. The external source $`\eta (x)`$ are $`c`$-number fields and hence taken the same value in both forward and backward branches. The strategy to obtain the equations of motion for small amplitude condensate perturbation $`\delta (x)`$, $`\delta ^{}(x)`$ is to consider the interaction Lagrangian $`_{\mathrm{int}}[\mathrm{\Psi }^{},\mathrm{\Psi },\chi ^{},\chi ]`$ given in (20) in perturbation theory and impose the *tadpole condition*

$$\chi ^\pm (x)_\eta =\chi ^\pm (x)_\eta =0$$
(50)
order by order in the perturbative expansion, but, consistent with linear response, only keep contributions linear in $`\delta (x)`$, $`\delta ^{}(x)`$.

Using the Feynman rules described in the previous section, to one-loop order the tadpole condition $`\chi ^+(0)_\eta =\chi ^+(0)_\eta =0`$ leads to the following expression

$$d^4x\chi ^+(0)\chi ^+(x)\left[\frac{\delta (x)}{g}+d^4x^{}\left[\mathrm{\Sigma }_{11}(xx^{})\delta (x^{})+\mathrm{\Sigma }_{12}(xx^{})\delta ^{}(x^{})\right]+𝒯(\mathrm{\Delta }_0,\mathrm{\Delta }_0^{})\eta (x)\right]=0,$$

$$d^4x\chi ^+(0)\chi ^+(x)\left[\frac{\delta ^{}(x)}{g}+d^4x^{}\left[\mathrm{\Sigma }_{22}(xx^{})\delta ^{}(x^{})+\mathrm{\Sigma }_{21}(xx^{})\delta (x^{})\right]+𝒯^{}(\mathrm{\Delta }_0,\mathrm{\Delta }_0^{})\eta ^{}(x)\right]=0,$$
(51)
where $`\mathrm{\Sigma }_{ab}(xx^{})`$ are the retarded self-energies of the pair fields and the *tadpole* $`𝒯(\mathrm{\Delta }_0,\mathrm{\Delta }_0^{})`$ denotes terms independent of the condensate perturbation $`\delta `$, $`\delta ^{}`$. The diagrams for $`\mathrm{\Sigma }_{11}(xx^{})`$ and $`\mathrm{\Sigma }_{12}(xx^{})`$ to one-loop order are depicted in Fig. 2, and those for $`\mathrm{\Sigma }_{22}(xx^{})`$ and $`\mathrm{\Sigma }_{21}(xx^{})`$ can be obtained in an analogous manner. Explicitly to one-loop order we find

$`\mathrm{\Sigma }_{11}(xx^{})`$ $`=`$ $`i\left[\mathrm{\Psi }^+(x)\sigma _{}\mathrm{\Psi }^+(x)\mathrm{\Psi }^+(x^{})\sigma _+\mathrm{\Psi }^+(x^{})\mathrm{\Psi }^+(x)\sigma _{}\mathrm{\Psi }^+(x)\mathrm{\Psi }^{}(x^{})\sigma _+\mathrm{\Psi }^{}(x^{})\right]`$
$`=`$ $`i\text{tr}\left[\sigma _{}S^{++}(xx^{})\sigma _+S^{++}(x^{}x)\sigma _{}S^+(xx^{})\sigma _+S^+(x^{}x)\right],`$
$`\mathrm{\Sigma }_{12}(xx^{})`$ $`=`$ $`i\left[\mathrm{\Psi }^+(x)\sigma _{}\mathrm{\Psi }^+(x)\mathrm{\Psi }^+(x^{})\sigma _{}\mathrm{\Psi }^+(x^{})\mathrm{\Psi }^+(x)\sigma _{}\mathrm{\Psi }^+(x)\mathrm{\Psi }^{}(x^{})\sigma _{}\mathrm{\Psi }^{}(x^{})\right]`$ (52)
$`=`$ $`i\text{tr}\left[\sigma _{}S^{++}(xx^{})\sigma _{}S^{++}(x^{}x)\sigma _{}S^+(xx^{})\sigma _{}S^+(x^{}x)\right],`$
where $`\mathrm{tr}`$ denotes the trace over the Nambu-Gor’kov indices. The expressions for $`\mathrm{\Sigma }_{22}(xx^{})`$ and $`\mathrm{\Sigma }_{21}(xx^{})`$ can be obtained, respectively, from expressions in (52) through the replacement $`\sigma _+\sigma _{}`$.

The diagrams for the tadpole $`𝒯(\mathrm{\Delta }_0,\mathrm{\Delta }_0^{})`$ are depicted in Fig. 3, from which we find

$`𝒯(\mathrm{\Delta }_0,\mathrm{\Delta }_0^{})`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Delta }_0}{g}}\mathrm{\Psi }^+(x)\sigma _{}\mathrm{\Psi }^+(x)`$ (53)
$`=`$ $`{\displaystyle \frac{\mathrm{\Delta }_0}{g}}+i\mathrm{tr}\left[\sigma _{}S^<(0)\right].`$
A simple calculation with the nonequilibrium Green’s functions of the Nambu-Gor’kov fields obtained in the previous section shows that the tadpole to one-loop order is given by

$$𝒯(\mathrm{\Delta }_0,\mathrm{\Delta }_0^{})=\mathrm{\Delta }_0\left[\frac{1}{g}\frac{d^3q}{(2\pi )^3}\frac{12n_\mathrm{F}(E_q)}{2E_q}\right],$$
(54)
where $`E_q`$ on the right-hand side of (54) is understood to be a function of $`|\mathrm{\Delta }_0|`$ \[see (46)\] and, hereafter, the momentum integral $`\frac{d^3q}{(2\pi )^3}`$ is restricted to the electronic states near the Fermi surface.

It is customary to rewrite the integral over the momentum $`𝐪`$ as being over the energy $`\xi `$ (measured from the Fermi surface) at the expense of introducing the density of states $`𝒩(\xi )`$, which is taken to be constant near the Fermi surface, and cutting the integral off at $`\pm \omega _\mathrm{D}`$ with $`\omega _\mathrm{D}`$ being the Debye energy, thus leading to

$$𝒯(\mathrm{\Delta }_0,\mathrm{\Delta }_0^{})=\mathrm{\Delta }_0\left[\frac{1}{g}𝒩(0)_0^{\omega _\mathrm{D}}\frac{d\xi }{E}[12n_\mathrm{F}(E)]\right],$$
(55)
where $`E=\sqrt{\xi ^2+|\mathrm{\Delta }_0|^2}`$ and $`𝒩(0)=mk_\mathrm{F}/2\pi ^2`$ is the density of states at the Fermi surface. Setting the external source $`\eta =0`$, the equilibrium condition for the homogeneous condensate $`𝒯(\mathrm{\Delta }_0,\mathrm{\Delta }_0^{})=0`$ becomes

$$\mathrm{\Delta }_0\left[\frac{1}{g}𝒩(0)_0^{\omega _\mathrm{D}}\frac{d\xi }{E}[12n_\mathrm{F}(E)]\right]=0.$$
(56)
In equilibrium and below the critical temperature, the condensate $`\mathrm{\Delta }_00`$. Therefore (56) leads to the finite-temperature BCS gap equation that determines $`\mathrm{\Delta }_0(T)`$:

$$g𝒩(0)_0^{\omega _\mathrm{D}}\frac{d\xi }{E}[12n_\mathrm{F}(E)]=1.$$
(57)
Introducing the space Fourier transforms for $`\delta (x)`$, $`\eta (x)`$ and $`\mathrm{\Sigma }_{ab}(xx^{})`$ as

$$\delta (𝐱,t)=\frac{d^3k}{(2\pi )^3}\delta _𝐤(t)e^{i𝐤𝐱},$$

etc., we find the equations of motion in momentum space to be given by

$$\delta _𝐤(t)+g_{\mathrm{}}^+\mathrm{}𝑑t^{}\left[\mathrm{\Sigma }_{11}(𝐤,tt^{})\delta _𝐤(t^{})+\mathrm{\Sigma }_{12}(𝐤,tt^{})\delta _𝐤^{}(t^{})\right]+g𝒯(\mathrm{\Delta }_0,\mathrm{\Delta }_0^{})\delta ^{(3)}(𝐤)g\eta _𝐤(t)=0,$$

$$\delta _𝐤^{}(t)+g_{\mathrm{}}^+\mathrm{}𝑑t^{}\left[\mathrm{\Sigma }_{22}(𝐤,tt^{})\delta _𝐤^{}(t^{})+\mathrm{\Sigma }_{21}(𝐤,tt^{})\delta _𝐤(t^{})\right]+g𝒯^{}(\mathrm{\Delta }_0,\mathrm{\Delta }_0^{})\delta ^{(3)}(𝐤)g\eta _𝐤^{}(t)=0.$$
(58)
The equations of motion obtained from the tadpole condition $`\chi ^{}(0)_\eta =\chi ^{}(0)_\eta =0`$ are the same as those given in the above expressions. While the equations of motion (58) are obtained to one-loop order, it is straightforward to conclude after a simple diagrammatic analysis that the structure of the equations of motion obtained above is general and valid to all orders in perturbation theory.

From the explicit expressions for the self-energies or by taking complex conjugation of the equation of motion for $`\delta _𝐤(t)`$, one can show that

$$\mathrm{\Sigma }_{21}(𝐤,tt^{})=\mathrm{\Sigma }_{12}^{}(𝐤,tt^{}),$$

$$\mathrm{\Sigma }_{22}(𝐤,tt^{})=\mathrm{\Sigma }_{11}^{}(𝐤,tt^{}).$$
(59)
Furthermore, rotational and parity invariance imply that the self-energies are only functions of $`k`$. Upon expressing the nonequilibrium Green’s functions in terms of the correlation functions $`S^{}`$, one finds immediately that the retarded self-energies have the following causal structure

$$\mathrm{\Sigma }_{11}(𝐤,tt^{})=\mathrm{\Sigma }_{11}^R(k,tt^{})\mathrm{\Theta }(tt^{}),$$

$$\mathrm{\Sigma }_{12}(𝐤,tt^{})=\mathrm{\Sigma }_{12}^R(k,tt^{})\mathrm{\Theta }(tt^{}).$$
(60)
Before proceeding further, we note that the invariance of the Lagrangian under the global gauge transformation (16) implies that *all* the terms in the equations of motion for $`\delta `$ ($`\delta ^{}`$) must transform as $`\delta `$ ($`\delta ^{}`$) itself under the gauge transformation. This in turn entails that the normal self-energies $`\mathrm{\Sigma }_{11}`$ and $`\mathrm{\Sigma }_{22}`$ must be invariant under the gauge transformation (16), while the anomalous ones $`\mathrm{\Sigma }_{12}`$ and $`\mathrm{\Sigma }_{21}`$ must transform as $`\mathrm{\Delta }_0^2`$ and $`\mathrm{\Delta }_0^2`$, respectively. Thus, it is proves convenient to write

$$\mathrm{\Sigma }_{11}^R(k,tt^{})=\mathrm{\Sigma }_{22}^R(k,tt^{})=\mathrm{\Sigma }_D(k,tt^{}),$$

$$\mathrm{\Sigma }_{12}^R(k,tt^{})=\mathrm{\Sigma }_{21}^R(k,tt^{})=(\mathrm{\Delta }_0/\mathrm{\Delta }_0^{})\mathrm{\Sigma }_O(k,tt^{}),$$
(61)
where both $`\mathrm{\Sigma }_D`$ and$`\mathrm{\Sigma }_O`$ are *invariant* under the gauge transformation (16). While rewriting the self-energies in this manner may seem a redundant exercise, the main point is to highlight and make explicit their transformation laws under the gauge transformation. This is an important aspect that needs to be addressed carefully in order to extract the Ward identities, an *exact* result of the underlying gauge symmetry to be explored below.

The gauge invariant self-energies $`\mathrm{\Sigma }_D`$ and $`\mathrm{\Sigma }_O`$ can be written in terms of their spectral representation as

$`\mathrm{\Sigma }_D(k,tt^{})={\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑\omega \left[\overline{A}(k,\omega )\mathrm{sin}\omega (tt^{})+i\overline{S}(k,\omega )\mathrm{cos}\omega (tt^{})\right],`$
$`\mathrm{\Sigma }_O(k,tt^{})={\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑\omega \widehat{A}(k,\omega )\mathrm{sin}\omega (tt^{}).`$ (62)
The symmetric ($`\overline{S}`$) and antisymmetric ($`\overline{A}`$, $`\widehat{A}`$) spectral functions are, respectively, even and odd functions of $`\omega `$ and are to be given by

$`\overline{S}(k,\omega )`$ $`=`$ $`{\displaystyle \frac{1}{4}}{\displaystyle }{\displaystyle \frac{d^3q}{(2\pi )^3}}\{[1n_\mathrm{F}(E_q)n_\mathrm{F}(E_p)]({\displaystyle \frac{\xi _q}{E_q}}+{\displaystyle \frac{\xi _p}{E_p}})[\delta (\omega E_qE_p)+\delta (\omega +E_q+E_p)]`$
$`[n_\mathrm{F}(E_q)n_\mathrm{F}(E_p)]({\displaystyle \frac{\xi _q}{E_q}}{\displaystyle \frac{\xi _p}{E_p}})[\delta (\omega E_q+E_p)+\delta (\omega +E_qE_p)]\},`$
$`\overline{A}(k,\omega )`$ $`=`$ $`{\displaystyle \frac{1}{4}}{\displaystyle }{\displaystyle \frac{d^3q}{(2\pi )^3}}\{[1n_\mathrm{F}(E_q)n_\mathrm{F}(E_p)](1+{\displaystyle \frac{\xi _q\xi _p}{E_qE_p}})[\delta (\omega E_qE_p)\delta (\omega +E_q+E_p)]`$
$`[n_\mathrm{F}(E_q)n_\mathrm{F}(E_p)](1{\displaystyle \frac{\xi _q\xi _p}{E_qE_p}})[\delta (\omega E_q+E_p)\delta (\omega +E_qE_p)]\},`$
$`\widehat{A}(k,\omega )`$ $`=`$ $`{\displaystyle \frac{1}{4}}{\displaystyle }{\displaystyle \frac{d^3q}{(2\pi )^3}}{\displaystyle \frac{|\mathrm{\Delta }_0|^2}{E_qE_p}}\{[1n_\mathrm{F}(E_q)n_\mathrm{F}(E_p)][\delta (\omega E_qE_p)\delta (\omega +E_q+E_p)]`$ (63)
$`+[n_\mathrm{F}(E_q)n_\mathrm{F}(E_p)][\delta (\omega E_q+E_p)\delta (\omega +E_qE_p)]\},`$
where $`p=|𝐤+𝐪|`$ and $`\xi _q=q^2/2m\mu `$.

The terms proportional to $`[1n_\mathrm{F}(E_q)n_\mathrm{F}(E_p)]`$ in the above expressions correspond to processes in which two quasiparticles are either created or destroyed, whereas the terms proportional to $`[n_\mathrm{F}(E_q)n_\mathrm{F}(E_p)]`$ arise from *Landau damping* processes corresponding to scattering of quasiparticle excitations in the medium.

In an experimental situation the dynamical evolution of the small amplitude condensate perturbation is studied by preparing a superconducting state slightly perturbed away from equilibrium by adiabatically coupling to some external source in the infinite past. Once the source is switched-off at time $`t=0`$ the perturbed condensate relaxes towards equilibrium and the relaxation dynamics is studied. As discussed above, this experimental situation can be realized within the real-time formulation described here by taking the spatial Fourier transform of the external source to be of the form

$$\eta _𝐤(t)=\eta _𝐤e^{ϵt}\mathrm{\Theta }(t),ϵ0^+.$$
(64)
The $`ϵ`$-term serves to switch on the source adiabatically from $`t=\mathrm{}`$ so as not to disturb the system too far from equilibrium in the process. If at $`t=\mathrm{}`$ the system was in an equilibrium state, then the condition of equilibrium (56) ensures that for $`t<0`$ there is a solution of the equations of motion (58) of the form

$$\delta _𝐤(t)=\delta _𝐤(0)e^{ϵt}\mathrm{for}t<0,$$
(65)
where $`\delta _𝐤(0)`$ is related to $`\eta _𝐤`$ through the equations of motion for $`t<0`$. The advantage of the adiabatic switching-on of the external source is that the time derivative of the solution (65) satisfies $`\dot{\delta }_𝐤(t<0)0`$ as $`ϵ0^+`$.

Let us introduce auxiliary quantities $`\mathrm{\Pi }_{ab}^R(k,tt^{})`$ defined as

$$\mathrm{\Sigma }_{ab}^R(k,tt^{})=_t^{}\mathrm{\Pi }_{ab}^R(k,tt^{}),$$
(66)
then upon using integration by parts, neglecting terms that vanish in the adiabatic limit $`ϵ0^+`$, and taking $`\mathrm{\Delta }_0`$ to be the equilibrium condensate and hence $`𝒯(\mathrm{\Delta }_0,\mathrm{\Delta }_0^{})=0`$, we find the equations of motion (58) become

$`[1+g\mathrm{\Pi }_{11}^R(k,0)]\delta _𝐤(t)+g\mathrm{\Pi }_{12}^R(k,0)\delta _𝐤^{}(t)g{\displaystyle _0^t}𝑑t^{}\left[\mathrm{\Pi }_{11}^R(k,tt^{})\dot{\delta }_𝐤(t^{})+\mathrm{\Pi }_{12}^R(k,tt^{})\dot{\delta }_𝐤^{}(t^{})\right]=0,`$
$`[1+g\mathrm{\Pi }_{22}^R(k,0)]\delta _𝐤^{}(t)+g\mathrm{\Pi }_{21}^R(k,0)\delta _𝐤(t)g{\displaystyle _0^t}𝑑t^{}\left[\mathrm{\Pi }_{22}^R(k,tt^{})\dot{\delta }_𝐤^{}(t^{})+\mathrm{\Pi }_{21}^R(k,tt^{})\dot{\delta }_𝐤(t^{})\right]=0.`$ (67)
The above coupled equations of motion for the condensate perturbation are now in the form of an *initial value problem* with initial conditions specified at $`t=0`$ and can be solved by Laplace transform. Introducing a two-component Nambu-Gor’kov spinor and its Laplace transform

$`\varphi _𝐤(t)=\left[\begin{array}{c}\delta _𝐤(t)\\ \delta _𝐤^{}(t)\end{array}\right],\stackrel{~}{\varphi }_𝐤(s)=\left[\begin{array}{c}\stackrel{~}{\delta }_𝐤(s)\\ \stackrel{~}{\delta }_𝐤^{}(s)\end{array}\right],`$ (68)
where $`s`$ is the Laplace variable, one can rewrite the Laplace transformed equations of motion in a compact matrix form as

$$\stackrel{~}{G}^1(k,s)\stackrel{~}{\varphi }_𝐤(s)=\frac{1}{s}\left[\stackrel{~}{G}^1(k,s)\stackrel{~}{G}^1(k,0)\right]\varphi _𝐤(0).$$
(69)
In the above equation, $`\stackrel{~}{G}^1(k,s)`$ is the inverse Green’s function (matrix) to one-loop order expressed in terms of the Laplace variable $`s`$

$$\stackrel{~}{G}^1(k,s)=\left[\begin{array}{cc}1+g\stackrel{~}{\mathrm{\Sigma }}_D(k,s)& g(\mathrm{\Delta }_0/\mathrm{\Delta }_0^{})\stackrel{~}{\mathrm{\Sigma }}_O(k,s)\\ g(\mathrm{\Delta }_0^{}/\mathrm{\Delta }_0)\stackrel{~}{\mathrm{\Sigma }}_O(k,s)& 1+g\stackrel{~}{\mathrm{\Sigma }}_D(k,s)\end{array}\right],$$
(70)
where $`\stackrel{~}{\mathrm{\Sigma }}_D`$ and $`\stackrel{~}{\mathrm{\Sigma }}_O`$ are the Laplace transforms of $`\mathrm{\Sigma }_D`$ and $`\mathrm{\Sigma }_O`$, respectively,

$`\stackrel{~}{\mathrm{\Sigma }}_D(k,s)={\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{dk_0}{k_0is}}\left[\overline{A}(k,k_0)+\overline{S}(k,k_0)\right],`$
$`\stackrel{~}{\mathrm{\Sigma }}_O(k,s)={\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{dk_0}{k_0is}}\widehat{A}(k,k_0).`$ (71)
The solution of (70) reads

$$\stackrel{~}{\varphi }_𝐤(s)=\frac{1}{s}[1\stackrel{~}{G}(k,s)\stackrel{~}{G}^1(k,0)]\varphi _𝐤(0),$$
(72)
where

$$\stackrel{~}{G}(k,s)=\frac{1}{\stackrel{~}{D}(k,s)}\left[\begin{array}{cc}1+g\stackrel{~}{\mathrm{\Sigma }}_D(k,s)& g(\mathrm{\Delta }_0/\mathrm{\Delta }_0^{})\stackrel{~}{\mathrm{\Sigma }}_O(k,s)\\ g(\mathrm{\Delta }_0^{}/\mathrm{\Delta }_0)\stackrel{~}{\mathrm{\Sigma }}_O(k,s)& 1+g\stackrel{~}{\mathrm{\Sigma }}_D(k,s)\end{array}\right],$$
(73)
with the denominator $`\stackrel{~}{D}(k,s)`$ given by

$$\stackrel{~}{D}(k,s)=[1+g\stackrel{~}{\mathrm{\Sigma }}_D(k,s)][1+g\stackrel{~}{\mathrm{\Sigma }}_D(k,s)][g\stackrel{~}{\mathrm{\Sigma }}_O(k,s)]^2.$$
(74)
The real-time evolution of the condensate perturbation $`\varphi _𝐤(t)`$ with an initial value $`\varphi _𝐤(0)`$ is now obtained from the inverse Laplace transform

$$\varphi _𝐤(t)=_{}\frac{ds}{2\pi i}e^{st}\stackrel{~}{\varphi }_𝐤(s),$$
(75)
where the Bromwich contour $``$ runs parallel to the imaginary axis in the complex $`s`$ plane to the right of all the singularities (poles and cuts) of $`\stackrel{~}{\varphi }_𝐤(s)`$ ivp . We note that there is *no* isolated pole in $`\stackrel{~}{\varphi }_𝐤(s)`$ at $`s=0`$ since the residue vanishes.

### IV.2 Self-energies in the long-wavelength, low-frequency limit

In this article we are interested to study the relaxation of long-wavelength, low-frequency fluctuations of the pair field, hence our next task is to expand the self-energies as a function of $`k,s`$ up to $`𝒪(k^2,s^2)`$. With $`\xi _q=q^2/2m\mu `$ and $`p=|𝐪+𝐤|`$ we use the following approximations

$$\xi _p\xi _q+\delta \xi ,\delta \xi =v_\mathrm{F}k\mathrm{cos}\theta ,$$

$$E_p=E_q+\xi _q\frac{\delta \xi }{E_q}+\frac{(\delta \xi )^2}{2}\frac{|\mathrm{\Delta }_0|^2}{E_q^3},$$
(76)
where $`\theta `$ is the angle between $`\widehat{𝐤}`$ and $`\widehat{𝐪}`$. We keep only up to terms of $`𝒪(k^2,s^2)`$ in the Laplace transform of the self-energies and obtain

$`\stackrel{~}{\mathrm{\Sigma }}_D(k,s)`$ $`=`$ $`𝒩(0){\displaystyle _0^{\omega _\mathrm{D}}}d\xi \{{\displaystyle \frac{12n_\mathrm{F}(E)}{2E}}[(1+{\displaystyle \frac{\xi ^2}{E^2}})(1{\displaystyle \frac{v_\mathrm{F}^2k^2}{12E^2}}{\displaystyle \frac{s^2}{4E^2}})(3{\displaystyle \frac{5\xi ^2}{E^2}}){\displaystyle \frac{v_\mathrm{F}^2k^2\xi ^2}{6E^4}}]`$
$`{\displaystyle \frac{n_\mathrm{F}(E)}{E}}{\displaystyle \frac{|\mathrm{\Delta }_0|^2}{E^2}}[{\displaystyle \frac{v_\mathrm{F}^2k^2\xi ^2}{2E^4}}+(1{\displaystyle \frac{3s^2}{2E^2}})(1{\displaystyle \frac{isE}{2v_\mathrm{F}k\xi }}\mathrm{ln}{\displaystyle \frac{is/v_\mathrm{F}k+\xi /E}{is/v_\mathrm{F}k\xi /E}})]\}+\mathrm{},`$
$`\stackrel{~}{\mathrm{\Sigma }}_O(k,s)`$ $`=`$ $`𝒩(0){\displaystyle _0^{\omega _\mathrm{D}}}d\xi {\displaystyle \frac{|\mathrm{\Delta }_0|^2}{E^2}}\{{\displaystyle \frac{12n_\mathrm{F}(E)}{2E}}[1({\displaystyle \frac{1}{4}}{\displaystyle \frac{5\xi ^2}{6E^2}}){\displaystyle \frac{v_\mathrm{F}^2k^2}{E^2}}{\displaystyle \frac{s^2}{4E^2}}]`$ (77)
$`{\displaystyle \frac{n_\mathrm{F}(E)}{E}}[{\displaystyle \frac{v_\mathrm{F}^2k^2}{2E^2}}({\displaystyle \frac{1}{3}}{\displaystyle \frac{\xi ^2}{E^2}})(1{\displaystyle \frac{3s^2}{2E^2}}+{\displaystyle \frac{s^2}{\xi ^2}})(1{\displaystyle \frac{isE}{2v_\mathrm{F}k\xi }}\mathrm{ln}{\displaystyle \frac{is/v_\mathrm{F}k+\xi /E}{is/v_\mathrm{F}k\xi /E}})]\}+\mathrm{},`$
where $`E=\sqrt{\xi ^2+|\mathrm{\Delta }_0|^2}`$ and the dots stand for terms of higher order in the ratios $`k/|\mathrm{\Delta }_0|,s/|\mathrm{\Delta }_0|`$. We note that to this order both $`\stackrel{~}{\mathrm{\Sigma }}_D(k,s)`$ and $`\stackrel{~}{\mathrm{\Sigma }}_O(k,s)`$ are even functions of the Laplace variable $`s`$. This important feature leads to the decoupling of the phase and amplitude fluctuations of the pair field in the equations of motion, as will be shown below.

### IV.3 One-loop effective action

The full equations of motion in real time (58) and their Laplace transform around the equilibrium state (69) allow us to obtain at once the *retarded effective action* in Fourier space. The retarded Green’s function kernel is obtained from $`\stackrel{~}{G}^1(k,s)`$ in (69) by the analytic continuation in the complex $`s`$ plane $`si\omega +0^+`$. In terms of the Fourier transform (retarded) of the two-component Nambu-Gor’kov spinor and the tadpole spinor

$$𝒯_𝐤(\omega )=\left[\begin{array}{c}𝒯(\mathrm{\Delta }_0,\mathrm{\Delta }_0^{})\\ 𝒯^{}(\mathrm{\Delta }_0,\mathrm{\Delta }_0^{})\end{array}\right]\delta ^{(3)}(𝐤)\delta (\omega ),$$
(78)
the retarded one-loop effective action to quadratic order in the fluctuations is therefore given by

$$S[\delta ,\delta ^{}]=S[0,0]+\frac{1}{2g}d^3k𝑑\omega \left[\varphi _𝐤^{}(\omega )G^1(k,\omega )\varphi _𝐤(\omega )+2\varphi _𝐤^{}(\omega )𝒯_𝐤(\omega )\right],$$
(79)
where $`G^1(k,\omega )=\stackrel{~}{G}^1(k,s=i\omega +0^+)`$ and $`S[0,0]`$ is a function of $`\mathrm{\Delta }_0`$, $`\mathrm{\Delta }_0^{}`$ such that

$$\frac{S[0,0]}{\mathrm{\Delta }_0^{}}=𝒯(\mathrm{\Delta }_0,\mathrm{\Delta }_0^{}),\frac{S[0,0]}{\mathrm{\Delta }_0}=𝒯^{}(\mathrm{\Delta }_0,\mathrm{\Delta }_0^{}).$$
(80)
Obviously, variational derivatives with respect to $`\delta ,\delta ^{}`$ reproduce the retarded equations of motion. We identify (79) with the one-loop effective action quadratic in the fluctuations.

## V Ward identity and static self-energies

The Ward identities, a consequence of the underlying (global) gauge invariance, are an integral part of the program to establish a connection with the Ginzburg-Landau description. Furthermore the equations of motion *must* fulfill these for consistency. In this section we show how the Ward identities emerge directly from the method described above used to obtain the equations of motion.

A straightforward diagrammatic analysis with the Feynman rules described above reveals that the generic structure of the equations of motion obtained via the tadpole method remains the same to all orders in perturbation theory. When combined with the transformation properties of $`\mathrm{\Delta }_0`$, $`\delta `$ and $`\chi `$ under the gauge transformation (16), this general form of the equations of motion allows to derive to all orders in perturbation theory the Ward identity for the tadpole $`𝒯(\mathrm{\Delta }_0,\mathrm{\Delta }_0^{})`$.

First, consider the case in which $`\delta ,\eta =0`$ the tadpole condition $`\chi ^+(0)_\eta =0`$ leads to $`𝒯(\mathrm{\Delta }_0,\mathrm{\Delta }_0^{})=0`$, which is the equilibrium condition for the homogeneous condensate. For a space-time independent shift of the condensate $`\mathrm{\Delta }_0\mathrm{\Delta }_0+\delta `$ induced by a space-time independent source $`\eta `$, the tadpole condition now leads to $`𝒯(\mathrm{\Delta }_0+\delta ,\mathrm{\Delta }_0^{}+\delta ^{})=\eta `$, which upon expanding to linear order in $`\delta `$ and $`\delta ^{}`$ becomes

$$𝒯(\mathrm{\Delta }_0,\mathrm{\Delta }_0^{})+\frac{𝒯(\mathrm{\Delta }_0,\mathrm{\Delta }_0^{})}{\mathrm{\Delta }_0}\delta +\frac{𝒯(\mathrm{\Delta }_0,\mathrm{\Delta }_0^{})}{\mathrm{\Delta }_0^{}}\delta ^{}=\eta .$$
(81)
We now compare (81) with the first equation of motion given in (58), which has the same generic structure as the full equation of motion obtained to all orders in perturbation theory, for space-time independent $`\delta `$ and $`\delta ^{}`$. Comparing the coefficients of $`\delta ,\delta ^{}`$ we recognize immediately that

$`{\displaystyle \frac{𝒯(\mathrm{\Delta }_0,\mathrm{\Delta }_0^{})}{\mathrm{\Delta }_0}}`$ $`=`$ $`{\displaystyle \frac{1}{g}}+\mathrm{\Sigma }_{11}(k=0,\omega =0),`$
$`{\displaystyle \frac{𝒯(\mathrm{\Delta }_0,\mathrm{\Delta }_0^{})}{\mathrm{\Delta }_0^{}}}`$ $`=`$ $`\mathrm{\Sigma }_{12}(k=0,\omega =0),`$ (82)
which relate the self-energies at *zero frequency* and *zero momentum* to the derivatives of the tadpole with respect to the condensate and is valid to *all orders in perturbation theory*.

The second important ingredient and which stems from the equations of motion (58) is that under a global gauge (phase) transformation (16) the tadpole $`𝒯(\mathrm{\Delta }_0,\mathrm{\Delta }_0^{})`$ transforms just as $`\delta `$, $`\mathrm{\Delta }_0`$ and $`\eta `$, i.e., $`𝒯(e^{i\theta }\mathrm{\Delta }_0,e^{i\theta }\mathrm{\Delta }_0^{})=e^{i\theta }𝒯(\mathrm{\Delta }_0,\mathrm{\Delta }_0^{})`$. Taking the gauge parameter $`\theta `$ to be infinitesimal and comparing the linear terms in $`\theta `$, we find to *all orders in perturbation theory* the Ward identity for the tadpole

$$𝒯(\mathrm{\Delta }_0,\mathrm{\Delta }_0^{})=\frac{𝒯(\mathrm{\Delta }_0,\mathrm{\Delta }_0^{})}{\mathrm{\Delta }_0}\mathrm{\Delta }_0\frac{𝒯(\mathrm{\Delta }_0,\mathrm{\Delta }_0^{})}{\mathrm{\Delta }_0^{}}\mathrm{\Delta }_0^{}.$$
(83)
Therefore upon combining (82) and the Ward identity (83), we obtain an alternative statement of the Ward identity which is an *exact* relationship between the tadpole and the self-energies at zero frequency and momentum

$$\mathrm{\Delta }_0\left[\frac{1}{g}+\mathrm{\Sigma }_{11}(0,0)\right]\mathrm{\Delta }_0^{}\mathrm{\Sigma }_{12}(0,0)=𝒯(\mathrm{\Delta }_0,\mathrm{\Delta }_0^{}),$$
(84)
where $`\mathrm{\Sigma }_{ab}(0,0)`$ stands for $`\mathrm{\Sigma }_{ab}(k=0,\omega =0)`$ for notational simplicity. Above the critical temperature, *both* the tadpole and the condensate vanish thus the above equation becomes a trivial identity. However, below the critical temperature $`\mathrm{\Delta }_00`$ and hence (84) leads to

$$\frac{g𝒯(\mathrm{\Delta }_0,\mathrm{\Delta }_0^{})}{\mathrm{\Delta }_0}=1+g\left[\mathrm{\Sigma }_{11}(0,0)\frac{\mathrm{\Delta }_0^{}}{\mathrm{\Delta }_0}\mathrm{\Sigma }_{12}(0,0)\right].$$
(85)
It is customary to choose the condensate to be real by redefining its phase via the gauge transformation (16). However, as argued above, for a condensate with an arbitrary phase, the anomalous self-energy $`\mathrm{\Sigma }_{12}`$ *must be* proportional to $`\mathrm{\Delta }_0^2`$ since in the equation of motion it multiplies $`\delta ^{}`$, which transforms under gauge transformations just as $`\mathrm{\Delta }_0^{}`$. This fact can be seen explicitly at the one-loop order in the expressions for the respective anomalous self-energy in (61) as well as in (62) with the spectral functions given by (63). Since the product $`\mathrm{\Sigma }_{12}\delta ^{}`$ must transform just as $`\delta `$ or $`\mathrm{\Delta }_0`$, therefore the phase $`\mathrm{\Delta }_0^{}/\mathrm{\Delta }_0`$ cancels the phase of $`\mathrm{\Delta }_0^2`$ in $`\mathrm{\Sigma }_{12}`$. In terms of the *gauge invariant* self-energies $`\mathrm{\Sigma }_D`$ and $`\mathrm{\Sigma }_O`$, the Ward identity (85) can be cast into an explicit gauge invariant form

$$\frac{g𝒯(\mathrm{\Delta }_0,\mathrm{\Delta }_0^{})}{\mathrm{\Delta }_0}=1+g\mathrm{\Sigma }_D(0,0)g\mathrm{\Sigma }_O(0,0).$$
(86)
We emphasize that the Ward identity (85) or (86) is an exact relationship valid *in* or *out* of equilibrium (corresponding to when the tadpole $`𝒯(\mathrm{\Delta }_0,\mathrm{\Delta }_0^{})`$ vanishes or not, respectively). In equilibrium the equilibrium condition $`𝒯(\mathrm{\Delta }_0,\mathrm{\Delta }_0^{})=0`$ implies that the self-energies in the static limit must satisfy the relation

$$1+g\mathrm{\Sigma }_D(0,0)g\mathrm{\Sigma }_O(0,0)=0.$$
(87)
The long-wavelength limit of the *static* self-energies is *defined* as

$`\mathrm{\Sigma }_D(0,0)=\underset{k0}{lim}\underset{s0}{lim}\mathrm{\Sigma }_D(k,s),`$
$`\mathrm{\Sigma }_O(0,0)=\underset{k0}{lim}\underset{s0}{lim}\mathrm{\Sigma }_O(k,s),`$ (88)
which, using (77), are found to one-loop order given by

$`\mathrm{\Sigma }_D(0,0)=𝒩(0){\displaystyle _0^{\omega _\mathrm{D}}}𝑑\xi \left[{\displaystyle \frac{12n_\mathrm{F}(E)}{2E}}\left(1+{\displaystyle \frac{\xi ^2}{E^2}}\right){\displaystyle \frac{n_\mathrm{F}(E)}{E}}{\displaystyle \frac{|\mathrm{\Delta }_0|^2}{E^2}}\right],`$
$`\mathrm{\Sigma }_O(0,0)=𝒩(0){\displaystyle _0^{\omega _\mathrm{D}}}𝑑\xi {\displaystyle \frac{|\mathrm{\Delta }_0|^2}{E^2}}\left[{\displaystyle \frac{12n_\mathrm{F}(E)}{2E}}+{\displaystyle \frac{n_\mathrm{F}(E)}{E}}\right].`$ (89)
The opposite limit $`lim_{s0}lim_{k0}\mathrm{\Sigma }_{D,O}(k,s)`$ yields the first terms in the expressions (89) but *not* the last terms proportional to $`n_\mathrm{F}(E)/E`$. This is a consequence of the non-analyticity of the Landau damping contributions to the self-energies.

Using the expression for the one-loop tadpole given by (55) it becomes clear that the identities (82) are fulfilled by the precise order of limits determined by (89), i.e., the long-wavelength limit of the *static* ($`s=0`$) self-energies.

We now show explicitly that the Ward identity (86) is fulfilled to one-loop order. The tadpole to one-loop order is given by (55). From the explicit form of the one-loop self-energies in the static limit (89), find that

$$\mathrm{\Sigma }_D(0,0)\mathrm{\Sigma }_O(0,0)=𝒩(0)_0^{\omega _\mathrm{D}}\frac{d\xi }{E}[12n_\mathrm{F}(E)],$$
(90)
which is an *exact* relationship to one-loop order, obtained for arbitrary $`\mathrm{\Delta }_0`$.

Upon collecting the above results, we find to one-loop order that

$`{\displaystyle \frac{g𝒯(\mathrm{\Delta }_0,\mathrm{\Delta }_0^{})}{\mathrm{\Delta }_0}}`$ $`=`$ $`1+g\mathrm{\Sigma }_D(0,0)g\mathrm{\Sigma }_O(0,0)`$ (91)
$`=`$ $`1g𝒩(0){\displaystyle _0^{\omega _\mathrm{D}}}{\displaystyle \frac{d\xi }{E}}[12n_\mathrm{F}(E)].`$
Therefore, the Ward identity (86) is manifestly fulfilled to one-loop order. This is an important advantage of the tadpole method of nonequilibrium field theory: the equations of motion obtained at a given order in the loop expansion are causal and guaranteed to fulfill the corresponding Ward identities to that order.

We highlight that whereas the Ward identity is independent of the limits, $`k,s0`$, the individual self-energies have different limits because of the non-analyticity of the Landau damping contribution. The relationship between the self-energies and the tadpole (82) is only valid in the static limit (89).

The results of this section, namely the Ward identity which relates the self-energies to derivatives of the tadpole will play an important role in the derivation of the effective Ginzburg-Landau theory.

## VI Time-dependent Ginzburg-Landau Theory

We are now in position to establish contact with the Ginzburg-Landau description of the long-wavelength, low-frequency excitations near the critical point. The Ginzburg-Landau description in terms of a functional of the order parameter is valid near the critical region when $`|TT_c|T_c,\mathrm{\Delta }(T)\left[T_c(T_cT)\right]^{1/2}T_c`$. The link between the Ginzburg-Landau (GL) theory and the microscopic theory of superconductivity in the static limit was established by Gor’kov gorkov59 .

Identifying the expectation value of the pair field $`\mathrm{\Delta }\mathrm{\Delta }(𝐱,t)`$ as the complex order parameter, the free energy in Ginzburg-Landau theory for an homogeneous order parameter (absence of gradients) is given by

$$(\mathrm{\Delta },\mathrm{\Delta }^{})=a(T)|\mathrm{\Delta }|^2+\frac{b}{2}|\mathrm{\Delta }|^4,$$
(92)
where near the critical temperature $`a(T)(TT_c)`$.

The linearized equation of motion of the small amplitude fluctuation of the condensate can be obtained by expanding the free energy around the equilibrium value of gap function; namely $`\mathrm{\Delta }_0`$ and $`\mathrm{\Delta }_0^{}`$. Writing $`\mathrm{\Delta }=\mathrm{\Delta }_0+\delta `$ and $`\mathrm{\Delta }^{}=\mathrm{\Delta }_0^{}+\delta ^{}`$, keeping only quadratic terms in the fluctuations:

$`(\mathrm{\Delta },\mathrm{\Delta }^{})=(\mathrm{\Delta }_0,\mathrm{\Delta }_0^{})+\delta {\displaystyle \frac{}{\mathrm{\Delta }}}+\delta ^{}{\displaystyle \frac{}{\mathrm{\Delta }^{}}}+\delta \delta ^{}{\displaystyle \frac{^2}{\mathrm{\Delta }\mathrm{\Delta }^{}}}+{\displaystyle \frac{\delta ^2}{2}}{\displaystyle \frac{^2}{\mathrm{\Delta }^2}}+{\displaystyle \frac{\delta ^2}{2}}{\displaystyle \frac{^2}{\mathrm{\Delta }^2}}.`$ (93)
The linearized equation of motion for both $`\delta `$ and $`\delta ^{}`$ can then be obtained by minimizing $``$ with respect to $`\delta ^{}`$ and $`\delta `$ respectively and are given by

$$\mathrm{\Sigma }_{11}^{\mathrm{GL}}\delta +\mathrm{\Sigma }_{12}^{\mathrm{GL}}\delta ^{}+𝒯^{\mathrm{GL}}(\mathrm{\Delta }_0,\mathrm{\Delta }_0^{})=0,$$

$$\mathrm{\Sigma }_{22}^{\mathrm{GL}}\delta ^{}+\mathrm{\Sigma }_{21}^{\mathrm{GL}}\delta +[𝒯^{\mathrm{GL}}(\mathrm{\Delta }_0,\mathrm{\Delta }_0^{})]^{}=0,$$
(94)
with

$$𝒯^{\mathrm{GL}}(\mathrm{\Delta }_0,\mathrm{\Delta }_0^{})=\mathrm{\Delta }_0\left[a(T)+b|\mathrm{\Delta }_0|^2\right],$$

$$\mathrm{\Sigma }_{11}^{\mathrm{GL}}=[\mathrm{\Sigma }_{22}^{\mathrm{GL}}]^{}=a(T)+2b|\mathrm{\Delta }_0|^2\frac{𝒯^{\mathrm{GL}}}{\mathrm{\Delta }_0},$$

$$\mathrm{\Sigma }_{12}^{\mathrm{GL}}=[\mathrm{\Sigma }_{21}^{\mathrm{GL}}]^{}=b\mathrm{\Delta }_0^2\frac{𝒯^{\mathrm{GL}}}{\mathrm{\Delta }_0^{}}.$$
(95)
The gap equation obtained from the Ginzburg-Landau free energy (92) is

$$𝒯^{\mathrm{GL}}(\mathrm{\Delta }_0,\mathrm{\Delta }_0^{})=0,$$
(96)
which determines the equilibrium value of the order parameter for $`\delta =\delta ^{}=0`$. When the gap equation (96) is fulfilled, the matrix of the second derivatives of the free energy (92) with respect to $`\delta ,\delta ^{}`$ has a zero eigenvalue with eigenvector determined by the relation

$$\delta ^{}=\frac{\mathrm{\Delta }_0^{}}{\mathrm{\Delta }_0}\delta .$$
(97)
Obviously this eigenvector with zero eigenvalue corresponds to a phase fluctuation and is the Goldstone boson, or the Anderson-Bogoliubov-Goldstone mode, associated with the broken global gauge symmetry.

The similarity between the equations of motion (94) and (58) for an homogeneous perturbation with vanishing external sources, as well as the similarities between the identities (95) and (82), (84) are now obvious and suggest the following identification

$`𝒯(\mathrm{\Delta }_0,\mathrm{\Delta }_0^{})`$ $``$ $`𝒯^{\mathrm{GL}}(\mathrm{\Delta }_0,\mathrm{\Delta }_0^{}),`$
$`{\displaystyle \frac{1}{g}}+\mathrm{\Sigma }_{11}(0,0)`$ $``$ $`\mathrm{\Sigma }_{11}^{\mathrm{GL}},`$
$`\mathrm{\Sigma }_{12}(0,0)`$ $``$ $`\mathrm{\Sigma }_{12}^{\mathrm{GL}},`$ (98)
where the self-energies on the left-hand side of the relations are understood as the long-wavelength limit of the *static* self energies as discussed in Sec. V. This similarity can be put on firmer footing by expanding the one-loop expression for the tadpole $`𝒯(\mathrm{\Delta }_0,\mathrm{\Delta }_0^{})`$ (55) in terms of a power series expansion in $`\mathrm{\Delta }_0,\mathrm{\Delta }_0^{}`$ and keeping up to third order consistently with the Ginzburg-Landau expansion.

Using the form of the tadpole given in (55) and the results obtained in Ref. ketterson, , we find ($`\mathrm{}=k_\mathrm{B}=1`$)

$$g𝒯(\mathrm{\Delta }_0,\mathrm{\Delta }_0^{})=\left[1g𝒩(0)\mathrm{ln}\left(\frac{1.14\omega _\mathrm{D}}{T}\right)\right]\mathrm{\Delta }_0+\frac{7g𝒩(0)\zeta (3)}{8\pi ^2T_c^2}\mathrm{\Delta }_0|\mathrm{\Delta }_0|^2+𝒪(\mathrm{\Delta }_0^5),$$
(99)
where $`\zeta `$ is the Riemann zeta function with $`\zeta (3)=1.202`$ and $`T_c=1.14\omega _\mathrm{D}e^{1/g𝒩(0)}`$. For $`TT_c`$ the coefficient of the linear term can be expanded as $`g𝒩(0)(TT_c)/T_c`$, therefore from the relations (82) and the definitions (61) we obtain

$`1+g\mathrm{\Sigma }_D(0,0)`$ $`=`$ $`g𝒩(0)\left({\displaystyle \frac{TT_c}{T_c}}\right)+{\displaystyle \frac{7g𝒩(0)\zeta (3)}{4\pi ^2T_c^2}}|\mathrm{\Delta }_0|^2,`$
$`g\mathrm{\Sigma }_O(0,0)`$ $`=`$ $`{\displaystyle \frac{7g𝒩(0)\zeta (3)}{8\pi ^2T_c^2}}|\mathrm{\Delta }_0|^2.`$ (100)
We now focus on the long-wavelength, low-frequency self-energies in the Ginzburg-Landau regime characterized by $`v_\mathrm{F}k,s|\mathrm{\Delta }_0|T`$. This regime corresponds to the description of the long-wavelength, low-frequency excitations near the critical point and will allow a consistent analysis of the real-time dynamics of (small) phase and amplitude fluctuations. We begin by writing the long-wavelength, low-frequency self-energies as \[see (77)\]

$`\stackrel{~}{\mathrm{\Sigma }}_D(k,s)`$ $`=`$ $`\mathrm{\Sigma }_D(0,0)+I(k,s),`$
$`\stackrel{~}{\mathrm{\Sigma }}_O(k,s)`$ $`=`$ $`\mathrm{\Sigma }_O(0,0)+J(k,s),`$ (101)
with $`\mathrm{\Sigma }_{D,O}(0,0)`$ given by (89). Since the self-energies are dimensionless functions of their arguments, the expansion in terms of $`k`$ and $`s`$ must involve ratios of these variables and the typical scales in the integrals. There are two, widely separated scales in the Ginzburg-Landau regime, namely, $`TT_c`$ and $`|\mathrm{\Delta }_0|`$ with $`|\mathrm{\Delta }_0|/T1`$. The expressions for $`I,J`$ feature energy denominators that would lead to infrared divergences if $`|\mathrm{\Delta }_0|`$ is set to zero. These divergences reflect the presence of inverse powers of $`|\mathrm{\Delta }_0|`$ in the expansions.

The leading terms in the expansion can be extracted by rescaling

$$\xi =|\mathrm{\Delta }_0|z,E=|\mathrm{\Delta }_0|ϵ,$$
(102)
with $`ϵ=\sqrt{z^2+1}`$ and introducing the dimensionless variables

$$x=\frac{|\mathrm{\Delta }_0|}{T},\kappa =\frac{v_\mathrm{F}k}{|\mathrm{\Delta }_0|},\overline{s}=\frac{s}{v_\mathrm{F}k},$$
(103)
in the integrals for $`I(k,s)`$ and $`J(k,s)`$. We obtain

$$I(k,s)=𝒩(0)x\left[\kappa ^2I_a(x)+\kappa ^2\overline{s}^2I_b(\overline{s},x)+i\overline{s}I_c(\overline{s},x)\right],$$
(104)
where

$`I_a(x)`$ $`=`$ $`{\displaystyle \frac{1}{4}}{\displaystyle _0^{\mathrm{}}}dz\{{\displaystyle \frac{2T(x,z)}{3xϵ^3}}[{\displaystyle \frac{2z^2}{ϵ^4}}{\displaystyle \frac{z^21}{4ϵ^2}}(1+{\displaystyle \frac{z^2}{ϵ^2}})]{\displaystyle \frac{C(x,z)z^2}{2ϵ^6}})\},`$
$`I_b(\overline{s},x)`$ $`=`$ $`{\displaystyle \frac{1}{4}}{\displaystyle _0^{\mathrm{}}}𝑑z\left\{{\displaystyle \frac{T(x,z)}{2xϵ^3}}\left(1+{\displaystyle \frac{z^2}{ϵ^2}}\right)+{\displaystyle \frac{3C(x,z)}{2ϵ^4}}\left[1{\displaystyle \frac{i\overline{s}ϵ}{2z}}\mathrm{ln}{\displaystyle \frac{i\overline{s}+z/ϵ}{i\overline{s}z/ϵ}}\right]\right\},`$
$`I_c(\overline{s},x)`$ $`=`$ $`{\displaystyle \frac{1}{4}}{\displaystyle _0^{\mathrm{}}}𝑑z\left[{\displaystyle \frac{C(x,z)}{2zϵ}}\mathrm{ln}{\displaystyle \frac{i\overline{s}+z/ϵ}{i\overline{s}z/ϵ}}\right],`$ (105)
with

$$T(x,z)=\mathrm{tanh}\left(\frac{xϵ}{2}\right),C(x,z)=\mathrm{cosh}^2\left(\frac{xϵ}{2}\right).$$
(106)
Similarly, for $`J(k,s)`$ we obtain

$$J(k,s)=𝒩(0)x\left[\kappa ^2J_a(x)+\kappa ^2\overline{s}^2J_b(\overline{s},x)+i\overline{s}J_c(\overline{s},x)\right],$$
(107)
where

$`J_a(x)`$ $`=`$ $`{\displaystyle \frac{1}{24}}{\displaystyle _0^{\mathrm{}}}𝑑z\left[T(x,z){\displaystyle \frac{7z^23}{xϵ^7}}C(x,z){\displaystyle \frac{2z^21}{ϵ^6}}\right],`$
$`J_b(\overline{s},x)`$ $`=`$ $`{\displaystyle \frac{1}{8}}{\displaystyle _0^{\mathrm{}}}𝑑z\left\{{\displaystyle \frac{T(x,z)}{xϵ^5}}C(x,z){\displaystyle \frac{2z^21}{z^2ϵ^4}}\left[1{\displaystyle \frac{i\overline{s}ϵ}{2z}}\mathrm{ln}{\displaystyle \frac{i\overline{s}+z/ϵ}{i\overline{s}z/ϵ}}\right]\right\},`$
$`J_c(\overline{s},x)`$ $`=`$ $`I_c(\overline{s},x).`$ (108)
An important consequence of the long-wavelength, low-frequency expansion of the self-energies in (101) together with (105), (108) is that to lowest order the self-energies are *even* functions of the Laplace variable $`s`$. This important aspect leads to the decoupling of the phase and amplitude fluctuations in the equations of motion, as can be seen as follows.

If we write the fluctuation of the order parameter around the space-time constant equilibrium solution of the tadpole equation (56) in the form

$`\mathrm{\Delta }(x,t)`$ $`=`$ $`\mathrm{\Delta }_0+\delta (x,t)`$ (109)
$``$ $`[\rho _0+\delta \rho (x,t)]e^{i\theta _0+i\delta \theta (x,t)},`$
where $`\delta \rho (x,t)`$ and $`\delta \theta (x,t)`$, respectively, are identified as the amplitude and phase fluctuations of the order parameter $`\mathrm{\Delta }(x,t)`$. It is convenient to introduce the following projection vectors

$`𝐩`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}i}}\left[\begin{array}{c}1\\ \frac{\mathrm{\Delta }_0^{}}{\mathrm{\Delta }_0}\end{array}\right],𝐩^{}={\displaystyle \frac{1}{\sqrt{2}i}}[1,{\displaystyle \frac{\mathrm{\Delta }_0}{\mathrm{\Delta }_0^{}}}],`$
$`𝐚`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left[\begin{array}{c}1\\ \frac{\mathrm{\Delta }_0^{}}{\mathrm{\Delta }_0}\end{array}\right],𝐚^{}={\displaystyle \frac{1}{\sqrt{2}}}[1,{\displaystyle \frac{\mathrm{\Delta }_0}{\mathrm{\Delta }_0^{}}}],`$ (110)
in terms of these vectors and the Nambu-Gor’kov spinors (68), the phase and amplitude fluctuations can be written as

$$\delta \theta (x,t)=\frac{1}{\sqrt{2}\mathrm{\Delta }_0}𝐩^{}\varphi (x,t),\frac{\delta \rho (x,t)}{\rho _0}=\frac{1}{\sqrt{2}\mathrm{\Delta }_0}𝐚^{}\varphi (x,t),$$
(111)
or, equivalently,

$$\stackrel{~}{\varphi }_𝐤(s)=\sqrt{2}\mathrm{\Delta }_0\left[\stackrel{~}{\delta \theta }_𝐤(s)𝐩+\frac{\stackrel{~}{\delta \rho }_𝐤(s)}{\rho _0}𝐚\right].$$
(112)
In terms of the Laplace transforms $`\stackrel{~}{\delta \theta }_𝐤(s)`$ and $`\stackrel{~}{\delta \rho }_𝐤(s)`$, the equations of motion (69) now become

$`\left[1+g\stackrel{~}{\mathrm{\Sigma }}_D(k,s)g\stackrel{~}{\mathrm{\Sigma }}_O(k,s)\right]\stackrel{~}{\delta \theta }_𝐤(s)={\displaystyle \frac{g}{s}}\left[\left(\stackrel{~}{\mathrm{\Sigma }}_D(k,s)\stackrel{~}{\mathrm{\Sigma }}_O(k,s)\right)\left(\stackrel{~}{\mathrm{\Sigma }}_D(k,0)\stackrel{~}{\mathrm{\Sigma }}_O(k,0)\right)\right]\delta \theta _𝐤(0),`$
$`\left[1+g\stackrel{~}{\mathrm{\Sigma }}_D(k,s)+g\stackrel{~}{\mathrm{\Sigma }}_O(k,s)\right]\stackrel{~}{\delta \rho }_𝐤(s)={\displaystyle \frac{g}{s}}\left[\left(\stackrel{~}{\mathrm{\Sigma }}_D(k,s)+\stackrel{~}{\mathrm{\Sigma }}_O(k,s)\right)\left(\stackrel{~}{\mathrm{\Sigma }}_D(k,0)+\stackrel{~}{\mathrm{\Sigma }}_O(k,0)\right)\right]\delta \rho _𝐤(0).`$ (113)
The decoupling of the equations of motion between phase and amplitude fluctuations is a direct consequence of the fact that self-energies are even functions of the Laplace variable to this, lowest order. In particular the equality $`\mathrm{\Sigma }_D(k,s)=\mathrm{\Sigma }_D(k,s)`$ to this order guarantees that the vectors $`𝐩`$ and $`𝐚`$ are eigenvectors of the matrix $`\stackrel{~}{G}^1(k,s)`$ in (69). The right-hand sides of (VI) are identified as the Laplace transform of the source terms $`J_\theta `$, $`J_\rho `$, respectively, that serve the purpose of setting up the initial value problem.

In particular to lowest order in the long-wavelength expansion, and after the analytic continuation $`si\omega +0^+`$, the *retarded* one-loop effective action for fluctuations around the equilibrium solution is given, up to a constant, by

$`S_{\mathrm{eff}}[\delta \theta ,\delta \rho ]`$ $`=`$ $`{\displaystyle \frac{|\mathrm{\Delta }_0|^2}{2g}}{\displaystyle d^3k𝑑\omega \left\{\delta \theta _𝐤(\omega )\left[1+g\mathrm{\Sigma }_D(k,\omega )g\mathrm{\Sigma }_O(k,\omega )\right]\delta \theta _𝐤(\omega )+2J_{\theta ,𝐤}(\omega )\delta \theta _𝐤(\omega )\right\}}`$ (114)
$`+{\displaystyle \frac{1}{2g}}{\displaystyle d^3k𝑑\omega \left\{\delta \rho _𝐤(\omega )\left[1+g\mathrm{\Sigma }_D(k,\omega )+g\mathrm{\Sigma }_O(k,\omega )\right]\delta \rho _𝐤(\omega )+2J_{\rho ,𝐤}(\omega )\delta \rho _𝐤(\omega )\right\}},`$
where $`\mathrm{\Sigma }_{D,O}(k,\omega )=\stackrel{~}{\mathrm{\Sigma }}_{D,O}(k,s=i\omega +0^+)`$ and the momentum integral must be understood to be restricted near the Fermi surface. Expanding the self-energies as in (101) together with (105), (108) and using the Ward identities (86) valid in equilibrium, one can easily show that the phase fluctuations are Goldstone modes. Using the explicit expressions for $`I(k,s)`$ and $`J(k,s)`$, we find that the effective action for the phase fluctuation is identical to that obtained in Refs. aitchison97, ; aitchison95,  to quadratic order, while the effective action for the amplitude fluctuation is a new contribution.

The approach followed here, based directly on the equations of motion in real time also allows us to obtain the effective action for the amplitude fluctuations on the same footing.

### VI.1 Phase fluctuations

The initial value problem for the phase and amplitude fluctuations described by the equations of motions (VI) can now be studied straightforwardly. For the phase fluctuation the inverse Laplace transform is obtained by integrating $`\stackrel{~}{\delta \theta }_𝐤(s)`$ along the Bromwich contour in the complex $`s`$ plane parallel to the imaginary axis and to the right of all the singularities of the Laplace transform

$$\frac{\stackrel{~}{\delta \theta }_𝐤(s)}{\delta \theta _𝐤(0)}=\frac{1}{s}\frac{F(\overline{s},x)}{F(\overline{s},x)+f(x)},$$
(115)
where

$`f(x)=I_a(x)J_a(x),`$
$`F(\overline{s},x)=\overline{s}^2\left[I_b(\overline{s},x)J_b(\overline{s},x)\right].`$ (116)
In the Ginzburg-Landau regime $`x1`$, we find from (105) and (108) that

$`f(x)`$ $``$ $`{\displaystyle \frac{x}{12}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dz}{z^3}}\left[\mathrm{tanh}{\displaystyle \frac{z}{2}}{\displaystyle \frac{z}{2}}\right]+𝒪(x^2)`$ (117)
$`=`$ $`{\displaystyle \frac{7\zeta (3)}{48\pi ^2}}x+𝒪(x^2).`$
Figure 4 displays a numerical evaluation of $`f(x)`$ with the full $`x`$ dependence. It shows clearly that $`f(x)1`$ and that $`f(x)`$ is accurately described by (117) in the Ginzburg-Landau regime.

To leading order in $`x`$, we can set $`x0`$ in the arguments when evaluating $`F(\overline{s},x)`$. After the analytic continuation $`si\omega +0^+`$ and introducing the dimensionless variable

$$\alpha =\frac{\omega }{v_\mathrm{F}k},$$
(118)
we find in the Ginzburg-Landau regime that the real and imaginary parts of $`F(\overline{s}=i\alpha +0^+,0)`$ are, respectively, given by (here we have set $`x0`$)

$`F_\mathrm{R}(\alpha )`$ $`=`$ $`{\displaystyle \frac{\alpha ^2}{8}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dz}{ϵ^2}}\left[1+{\displaystyle \frac{1}{z^2}}\left(1{\displaystyle \frac{\alpha ϵ}{2z}}\mathrm{ln}\left|{\displaystyle \frac{\alpha +\frac{z}{ϵ}}{\alpha \frac{z}{ϵ}}}\right|\right)\right]`$
$`\stackrel{\alpha 1}{}`$ $`{\displaystyle \frac{\pi ^2|\alpha |^3}{64}},`$
$`F_\mathrm{I}(\alpha )`$ $`=`$ $`{\displaystyle \frac{\pi \alpha ^3}{16}}\mathrm{\Theta }(1\alpha ^2){\displaystyle _{z_{\mathrm{min}}}^{\mathrm{}}}{\displaystyle \frac{dz}{z^3ϵ}}`$ (119)
$`\stackrel{\alpha 1}{}`$ $`{\displaystyle \frac{\pi \alpha }{32}},`$
where $`z_{\mathrm{min}}=\sqrt{\alpha ^2/(1\alpha ^2)}`$ and we have set $`x0`$ in the argument of $`F(\overline{s},x)`$. These expressions agree with those found by Aitchison et al. aitchison97  in the limit $`x0`$ comparison1 .

The real \[$`F_\mathrm{R}(\alpha )`$\] and imaginary \[$`F_\mathrm{I}(\alpha )`$\] parts of $`F(\alpha )`$ for $`0\alpha 1`$ are displayed in Fig. 5. The real part is a monotonically increasing function of $`\alpha `$ while the imaginary part only has contribution in the Landau damping cut $`1<\alpha <1`$ corresponding to $`v_\mathrm{F}k<\omega <v_\mathrm{F}k`$.

Isolated *real* poles of the Laplace transform, describing stable quasiparticle excitations, correspond to the solutions of the following equations

$$F_\mathrm{R}(\alpha )f(x)=0,F_\mathrm{I}(\alpha )=0.$$
(120)
It is clear from (117) and Fig. 5 that the first equation in (120) can be fulfilled for $`f(x)1`$ only in a region where $`F_\mathrm{I}(\alpha )0`$. Therefore in the Ginzburg-Landau regime the long-wavelength, low-frequency Laplace transform (115) has *no* isolated quasiparticle poles, the only singularity is a branch cut in the imaginary axis $`iv_\mathrm{F}k<s<iv_\mathrm{F}k`$, a consequence of Landau damping.

Having understood the analytic structure of the Laplace transform, we can now proceed to study the time evolution via the inverse Laplace transform (75) by closing the Bromwich contour wrapping around the cut $`iv_\mathrm{F}k<s<iv_\mathrm{F}k`$. We obtain

$$\delta \theta _𝐤(t)=_1^1\frac{d\alpha }{\alpha }\rho _{\mathrm{ph}}(\alpha )\mathrm{cos}(\alpha \tau )\delta \theta _𝐤(0),$$
(121)
where $`\tau =v_\mathrm{F}kt`$ and the spectral density for the phase fluctuations is given by

$$\rho _{\mathrm{ph}}(\alpha )=\frac{1}{\pi }\frac{f(x)F_\mathrm{I}(\alpha )}{\left[F_\mathrm{R}(\alpha )f(x)\right]^2+F_\mathrm{I}^2(\alpha )}.$$
(122)
The spectral density $`\rho _{\mathrm{ph}}(\alpha )`$ and the real-time evolution of $`\delta \theta _𝐤(t)`$ for several values of $`x1`$ are displayed in Figs. 6 and 7, respectively. Three important features are gleaned from these figures:

* The spectral density features a sharp peak at a value $`\alpha _{\mathrm{peak}}(x)`$ that vanishes continuously as $`x0`$. This peak would indicate quasiparticle excitations with dispersion relation $`\omega =\alpha _{\mathrm{peak}}(x)v_\mathrm{F}k`$. The group velocity $`v_\mathrm{g}(x)14\zeta (3)v_\mathrm{F}x/3\pi ^3`$ vanishes at the critical point $`x=0`$ and increases linearly with $`x`$ at least within the range $`0<x0.2`$.
* While the spectral density is not of the Breit-Wigner type (Lorentzian or resonance) and hence a true width cannot be extracted unambiguously, it is clear from the figure that qualitatively the quasiparticle excitations have a “width”. This width vanishes at the critical point and increases monotonically with $`x1`$ at least within the range consistent with a Ginzburg-Landau expansion. Therefore these quasiparticles will be Landau damped below $`T_c`$ and the relaxation time scale (the inverse of the “width”) diverges at the critical point $`x=0`$. This expectation will be confirmed below and suggests critical slowing down of long-wavelength phase fluctuations.
* The reason that despite the appearance of a sharp peak in the spectral density near the critical point the real-time dynamics does not reveal the oscillations associated with a “quasiparticle pole” is clear. As $`x0`$ *both* the damping rate and the group velocity $`v_\mathrm{g}(x)`$ vanish in such a way that the time scale for damping is either shorter or of the same order as the time scale for the oscillation.

Furthermore, (121) evaluated at $`t=0`$ leads to the sum rule

$$_1^1\frac{d\alpha }{\alpha }\rho _{\mathrm{ph}}(\alpha )=1,$$
(123)
which we have confirmed numerically for a wide range of $`x`$.

It is clear that damping becomes more pronounced for larger $`x`$ and while the peak would seem to lead to oscillations with period $`2\pi /\alpha _{\mathrm{peak}}(x)`$ there are no hints of oscillatory behavior in the real-time evolution. Phase fluctuations are strongly overdamped without featuring a propagating mode. Hence we conclude that near the critical point, in the Ginzburg-Landau regime $`x1`$, Goldstone modes or phase fluctuations are severely damped despite the fact that the spectral density features a peak that would indicate a quasiparticle “dispersion relation” $`\omega =\alpha _{\mathrm{peak}}(x)k`$. The damping becomes larger for larger $`x`$ and is solely a consequence of *collisionless* Landau damping.

For $`x1`$, we find that the nonequilibrium relaxation of the phase fluctuation is very well approximated by an exponential (see the logarithmic plot in Fig. 7). The numerical analysis clearly indicates that $`\delta \theta _𝐤(t)\delta \theta _𝐤(0)e^{\gamma _k(x)t}`$, where from (117) and (119) the damping rate $`\gamma _k(x)`$ is found to be given by

$$\gamma _k(x)\frac{14\zeta (3)v_\mathrm{F}k}{3\pi ^3}x\mathrm{for}x1.$$
(124)
This result reveals *critical slowing down* since the damping rate vanishes at the critical point. Furthermore we also see that for fixed $`x`$ the damping rate also vanishes in the long-wavelength limit, in agreement with the expectation that the relaxation time scale of Goldstone bosons should diverge in the long-wavelength limit. This is one of the novel results of this study: The long-wavelength phase fluctuations are *overdamped* by Landau damping in the Ginzburg-Landau regime but the damping rate vanishes at the critical point indicating critical slowing down.

In the Ginzburg-Landau regime and for long-wavelength, low-frequency fluctuations the nonequilibrium retarded Ginzburg-Landau effective action for phase fluctuations (Goldstone modes) is given to lowest order by (we have set $`J_\theta =0`$)

$$S_{\mathrm{ph}}^{\mathrm{GL}}[\delta \theta ]=𝒩(0)\frac{|\mathrm{\Delta }_0|}{2T_c}d^3k𝑑\omega \left\{\delta \theta _𝐤(\omega )(v_\mathrm{F}k)^2\left[F(k,\omega )f(x)\right]\delta \theta _𝐤(\omega )\right\},$$
(125)
up to an additive constant, where $`f(x)`$ and $`F(k,\omega )=F(\alpha =\omega /v_\mathrm{F}k)`$ are given by (117) and (119), respectively. The imaginary part of $`F(\omega ,k)`$ originates in Landau damping. This long-wavelength, low-frequency effective action leads to the equations of motion for long-wavelength phase fluctuations in the linearized approximation valid in the Ginzburg-Landau regime. Thus it can be genuinely called the effective time dependent Ginzburg-Landau effective action. It is nonlocal in time as a consequence of Landau damping and describes real-time relaxation which is completely overdamped. The damping rate vanishes at the critical temperature thus signaling critical slowing down.

#### Away from the Ginzburg-Landau regime: $`x1`$

While we have focused in the Ginzburg-Landau regime $`x1`$, we now study the region $`x1`$ *away* from the domain of validity of the Ginzburg-Landau expansion mainly with the purpose of comparing our results to those obtained in Refs. aitchison97, ; aitchison95, .

In this case we must keep the full $`x`$ dependence of the functions $`I_{a,b}`$ and $`J_{a,b}`$. The spectral density has the same form as in (122) but now with the replacement $`F(\alpha )F(\overline{s}=i\alpha +0^+,x)`$, whose real \[$`F_\mathrm{R}(\alpha ,x)`$\] and imaginary \[$`F_\mathrm{I}(\alpha ,x)`$\] parts can be found straightforwardly as in the previous case. In the low-temperature limit $`x1`$, we obtain

$`f(x)`$ $``$ $`{\displaystyle \frac{1}{12x}},`$
$`F_\mathrm{R}(\alpha ,x)`$ $``$ $`{\displaystyle \frac{\alpha ^2}{4x}},`$
$`F_\mathrm{I}(\alpha ,x)`$ $``$ $`{\displaystyle \frac{\pi \alpha }{8}}\left[\sqrt{1\alpha ^2}e^{x/\sqrt{1\alpha ^2}}+{\displaystyle \frac{\alpha ^2x}{2}}\left[e^x\mathrm{Ei}(x_+)+e^x\mathrm{Ei}(x_{})\right]\right],`$ (126)
where $`\mathrm{Ei}(x)`$ is the exponential integral function, $`x_\pm =x\left(\frac{1}{\sqrt{1\alpha ^2}}\pm 1\right)`$ and the exponentially small temperature corrections to $`f(x)`$ and $`F_\mathrm{R}(\alpha ,x)`$ have been neglected. The spectral density with the full $`x`$ dependence is displayed in Fig. 8 for $`x1`$ away from the regime of validity of the Ginzburg-Landau approximation. It shows clearly the emergence of a sharp quasiparticle peak, which for $`x1`$ is at $`\alpha _{\mathrm{peak}}=1/\sqrt{3}`$ in agreement with the results of Refs. aitchison97, ; aitchison95, . The real-time evolution of phase fluctuations in this regime is displayed in Fig. 9. It is clear from these figures that the sharp quasiparticle peak in the spectral density results in a real-time dynamics that is weakly *underdamped* by Landau damping. This is in contrast to the real-time dynamics in the Ginzburg-Landau regime, where Landau damping is so severe that the phase fluctuation is overdamped and hence there is no quasiparticle interpretation.

Thus the real-time evolution displayed above for this case confirms the results of Refs. aitchison97, ; aitchison95,  valid well below the critical temperature, that a narrow quasiparticle peak describes the dynamics of phase fluctuations. As in the previous case discussed above, the damping rate is found to be given by

$$\gamma _k(x)\frac{5\pi v_\mathrm{F}k}{6}e^{\sqrt{\frac{3}{2}}x}\mathrm{for}x1,$$
(127)
where the factor $`v_\mathrm{F}k`$ is again a consequence of the Goldstone nature of the phase fluctuation, leading to a vanishing damping rate in the long-wavelength limit. Because the effects of Landau damping are suppressed in the low-temperature limit ($`x1`$), the damping rate becomes smaller as the temperature decreases and hence the oscillatory behavior associated with the “quasiparticle pole” is evidenced. In this region a *local* time-dependent effective theory is a good approximate description of the nonequilibrium dynamics. On the contrary, the real-time dynamics in the Ginzburg-Landau regime ($`x1`$) is *overdamped*, dominated by Landau damping and *cannot* be accurately described by a local effective action in real time. Thus the Ginzburg-Landau dynamics is *purely dissipative*.

### VI.2 Amplitude fluctuations

For the amplitude fluctuations in the long-wavelength, low-frequency limit the equation of motion in terms of the Laplace transform requires the inverse propagator $`1+g\stackrel{~}{\mathrm{\Sigma }}_D(k,s)+g\stackrel{~}{\mathrm{\Sigma }}_O(k,s)`$. In the case of amplitude fluctuations we expect a “gap” of order $`|\mathrm{\Delta }_0|`$ in the spectrum of the quasiparticle excitations. While this pole is away from the region of validity of the long-wavelength, low-frequency approximation we can, nevertheless, obtain a qualitative if not a quantitative estimate of the dispersion relation for amplitude fluctuations.

*Isolated poles*: Since the imaginary parts of $`\mathrm{\Sigma }_D`$, $`\mathrm{\Sigma }_O`$ are nonzero only for $`v_\mathrm{F}k<\omega <v_\mathrm{F}k`$, for $`k0`$ and $`\omega |\mathrm{\Delta }_0|`$ the expected quasiparticle pole will be away from the continuum. We can find the position of this pole by looking for the solutions of

$$\mathrm{Re}\left[1+g\stackrel{~}{\mathrm{\Sigma }}_D(k,s)+g\stackrel{~}{\mathrm{\Sigma }}_O(k,s)\right]_{s=i\omega +0^+}=0.$$
(128)
This equation has solutions for values of $`\omega `$ given by $`\omega _{\mathrm{amp}}(k)`$, which determine the dispersion relation for the amplitude fluctuation.

In the Ginzburg-Landau regime the gap of the spectrum can be estimated by setting $`k=0`$ in the expressions of the self-energies and keeping the lowest order \[$`𝒪(s^2)`$\] terms in the expansion of the self-energies. After rescaling variables in the integrals, we obtain

$$1+g\stackrel{~}{\mathrm{\Sigma }}_D(0,s)+g\stackrel{~}{\mathrm{\Sigma }}_O(0,s)g𝒩(0)_0^{\mathrm{}}𝑑z\frac{T(x,z)}{ϵ^3}\left[1+\frac{z^2s^2}{4ϵ^2|\mathrm{\Delta }_0|^2}\right],$$
(129)
where use has been made of the gap equation. Keeping the lowest order in $`x1`$ in the Ginzburg-Landau regime, the integrals can be done straightforwardly and we find

$$1+g\stackrel{~}{\mathrm{\Sigma }}_D(0,s)+g\stackrel{~}{\mathrm{\Sigma }}_O(0,s)g𝒩(0)\frac{\pi x}{4}\left[1+\frac{s^2}{8|\mathrm{\Delta }_0|^2}\right],$$
(130)
which, upon the analytic continuation $`si\omega +0^+`$, suggests the gap of the spectrum in the Ginzburg-Landau regime to be given by $`2\sqrt{2}|\mathrm{\Delta }_0|`$.

Away from the Ginzburg-Landau regime and for arbitrary values of $`k`$ the equation for the dispersion relation, (128) must be solved numerically. Figure 10 displays the dispersion relation for several values of $`x`$ away from the Ginzburg-Landau regime. The values of the gap shown are consistent with the results obtained by Aitchison et al. comparison2 . We haste to emphasize, however, that the position of these single (quasi)particle poles are *away* from the regime of validity of our approximations and must only be taken as indicative and consistent with the findings of Ref. aitchison97,  but *not* as accurate dispersion relations since higher orders in the ratio $`s/|\mathrm{\Delta }_0|`$ will modify these results.

This analysis is included here with the sole purposes of (i) emphasizing that there are single (quasi)particle poles away from the Landau damping continuum, consistent with the expectation of a gap in the spectrum of small amplitude fluctuations, (ii) establishing a comparison with the results of Ref. aitchison97, , and (iii) offering at least a qualitative, if not a reliable quantitative, discussion of the terms contributing to the time evolution of amplitude fluctuations. The full dispersion relations must be obtained by keeping *all* terms in the self-energies which will involve a substantial numerical effort, a task clearly beyond the scope of this article whose focus is on the Ginzburg-Landau regime in the long-wavelength, low-frequency limit.

*Landau damping cut*: The imaginary part associated with Landau damping arises for $`k0`$ and is nonvanishing only along the Landau damping cut $`v_\mathrm{F}k<\omega <v_\mathrm{F}k`$. Focusing on the long-wavelength, low-frequency limit and in the Ginzburg-Landau regime, we find

$`1+g\stackrel{~}{\mathrm{\Sigma }}_D(k,s)+g\stackrel{~}{\mathrm{\Sigma }}_O(k,s)`$ $``$ $`2g\left[\mathrm{\Sigma }_O(0,0)+i\overline{s}x𝒩(0)I_c(\overline{s},x)\right],`$ (131)
where we have used the Ward identity (85) and neglected the contributions from $`I_{a,b}`$ and $`J_{a,b}`$, which are subdominant for $`v_\mathrm{F}k,s|\mathrm{\Delta }_0|`$. We note that in the case of phase fluctuations the terms with $`I_c`$ and $`J_c`$ cancel each other in the difference in the self-energies because $`I_c=J_c`$, however, for amplitude fluctuations these terms add up and furnish the dominant contribution.

We now proceed to study the real-time evolution of the amplitude fluctuation by integrating $`\stackrel{~}{\delta \rho }_𝐤(s)`$ along the Bromwich contour in the complex $`s`$ plane parallel to the imaginary axis and to the right of all the singularities of the Laplace transform

$$\frac{\stackrel{~}{\delta \rho }_𝐤(s)}{\delta \rho _𝐤(0)}=\frac{1}{s}\frac{H(\overline{s},x)}{H(\overline{s},x)+h},$$
(132)
where

$$xh=\mathrm{\Sigma }_O(0,0)/𝒩(0),H(\overline{s},x)=i\overline{s}I_c(\overline{s},x).$$
(133)
We obtain

$$\delta \rho _𝐤(t)=_{\mathrm{}}^{\mathrm{}}\frac{d\omega }{\omega }\rho _{\mathrm{amp}}(\omega ,k)\mathrm{cos}(\omega t)\delta \rho _𝐤(0),$$
(134)
where $`\rho _{\mathrm{amp}}(\omega ,k)`$ is the spectral density for the amplitude fluctuation

$`\rho _{\mathrm{amp}}(\omega ,k)`$ $`=`$ $`\mathrm{sgn}(\omega )Z(k)\omega _{\mathrm{amp}}^2(k)\delta [\omega ^2\omega _{\mathrm{amp}}^2(k)]+\rho _{\mathrm{amp}}^{\mathrm{cut}}(\alpha ),`$
$`\rho _{\mathrm{amp}}^{\mathrm{cut}}(\alpha )`$ $`=`$ $`{\displaystyle \frac{1}{\pi }}\mathrm{Im}\left[{\displaystyle \frac{h}{h+H(\overline{s},x)}}\right]_{\overline{s}=i\alpha +0^+}.`$ (135)
In the above expression the first term arises from quasiparticle pole $`\omega _{\mathrm{amp}}(k)`$ and the second term arises from the Landau damping cut $`v_\mathrm{F}k<\omega <v_\mathrm{F}k`$.

The quasiparticle pole will contribute an undamped oscillatory component to the time evolution given by

$$\delta \rho _𝐤^{\mathrm{pole}}(t)=Z(k)\mathrm{cos}[\omega _{\mathrm{amp}}(k)t].$$
(136)
The residue of the quasiparticle pole $`Z(k)`$ is determined by

$$Z(k)=1_1^1\frac{d\alpha }{\alpha }\rho _{\mathrm{amp}}^{\mathrm{cut}}(\alpha ).$$
(137)
Figure 12 shows the temperature dependence of $`Z(k)`$ in the long-wavelength limit $`v_\mathrm{F}k|\mathrm{\Delta }_0|`$. It reveals clearly that in the Ginzburg-Landau regime ($`x1`$) the spectral density $`\rho _{\mathrm{amp}}(k,\omega )`$ is dominated by the Landau damping cut.

Whereas the damping rate of the quasiparticle vanishes in the long-wavelength, low-frequency approximation, we expect that higher order contributions, in particular the decay into a pair of Bogoliubov quasiparticles (bogolons) ketterson , will lead to a width of the quasiparticle pole and hence a finite damping rate if the gap the quasiparticle spectrum is greater than $`2|\mathrm{\Delta }_0|`$, as seems to be the case from the previous analysis. However, a more comprehensive study of the dispersion relation is needed before reaching a quantitative conclusion. Again this is beyond the regime of validity of the long-wavelength, low-frequency approximation studied here.

In the Ginzburg-Landau regime $`x1`$, we find \[see (VI)\]

$$h=\frac{7\zeta (3)}{8\pi ^2}x,$$
(138)
and the real and imaginary parts of $`H(\overline{s}=i\alpha +0^+,0)`$ to be given by (here we have set $`x=0`$)

$`H_\mathrm{R}(\alpha )`$ $`=`$ $`{\displaystyle \frac{\alpha }{8}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dz}{zϵ}}\mathrm{ln}\left|{\displaystyle \frac{\alpha +\frac{z}{ϵ}}{\alpha \frac{z}{ϵ}}}\right|,`$
$`\stackrel{\alpha 1}{}`$ $`{\displaystyle \frac{\pi ^2|\alpha |}{16}},`$
$`H_\mathrm{I}(\alpha )`$ $`=`$ $`{\displaystyle \frac{\pi \alpha }{8}}\mathrm{\Theta }(1\alpha ^2){\displaystyle _{z_{\mathrm{min}}}^{\mathrm{}}}{\displaystyle \frac{dz}{zϵ}},`$ (139)
$`\stackrel{\alpha 1}{}`$ $`{\displaystyle \frac{\pi \alpha }{8}}\mathrm{ln}{\displaystyle \frac{2}{\alpha }},`$
respectively, where $`z_{\mathrm{min}}=\sqrt{\alpha ^2/(1\alpha ^2)}`$. The real and imaginary parts of $`H(\alpha )`$ are displayed in Fig. 11. The Landau damping contribution to the spectral density in the Ginzburg-Landau regime and the long-wavelength, low-frequency limit is accurately described by

$$\rho _{\mathrm{amp}}^{\mathrm{cut}}(\alpha )=\frac{h|H_\mathrm{I}(\alpha )|}{\left[h+H_\mathrm{R}(\alpha )\right]^2+H_\mathrm{I}^2(\alpha )}.$$
(140)
which, as displayed in Fig. 13, clearly reveals a sharp peak near $`\alpha 0`$.

The real-time evolution of the Landau damping contribution to the amplitude fluctuation $`\delta \rho _𝐤^{\mathrm{cut}}(t)`$ is displayed in Fig. 14. While for $`x0.05`$ the nonequilibrium relaxation through Landau damping can be approximated by an exponential, clearly such is not the case for $`x0.05`$. Nevertheless we still find that the Landau damping relaxation time scale becomes *longer* as $`x0`$, revealing that this contribution to the nonequilibrium relaxation seems to be critically slowed down near the critical point.

We close this section, by providing the nonlocal Landau damping contribution to the *retarded* effective action for long-wavelength, low-frequency amplitude fluctuations in the Ginzburg-Landau regime (we have set $`J_\rho =0`$)

$$S_{\mathrm{amp}}^{\mathrm{GL}}[\delta \rho ]=2𝒩(0)\frac{|\mathrm{\Delta }_0|}{T_c}d^3k𝑑\omega \delta \rho _𝐤(\omega )\left[\frac{7\zeta (3)|\mathrm{\Delta }_0|}{8\pi ^2T}+H(\omega ,k)\right]\delta \rho _𝐤(\omega ),$$
(141)
up to a additive constant, where $`H(\omega ,k)=H(\alpha =\omega /v_\mathrm{F}k)`$ given by (139). Obviously, $`S_{\mathrm{amp}}^{\mathrm{GL}}[\delta \rho ]`$ leads to the retarded equations of motion for the long-wavelength, low-frequency amplitude fluctuations.

## VII Conclusions

In this article we have focused on the real-time nonequilibrium dynamics of small phase and amplitude fluctuations of the order parameter in neutral BCS superconductors in the Ginzburg-Landau regime near the critical point. We have implemented the Schwinger-Keldysh formulation of nonequilibrium field theory combined with the novel tadpole method to obtain directly in real time the *retarded* equations of motion for small fluctuations around equilibrium. These equations allow to extract the one-loop effective action for the long-wavelength, low-frequency phase and amplitude fluctuations in the Ginzburg-Landau regime, which is characterized by $`|\mathrm{\Delta }_0(T)|/T1`$ with $`\mathrm{\Delta }_0(T)`$ the finite-temperature gap. Furthermore, the retarded equations of motion can be cast as an initial value problem to study the relaxation of nonequilibrium fluctuations directly in real time.

We studied in detail the relaxation of phase fluctuations within and away from the Ginzburg-Landau regime. Despite the fact that the spectral density features a sharp peak with a Goldstone-like dispersion relation in the Ginzburg-Landau regime the relaxation is completely *overdamped* as a consequence of Landau damping. This is consistent with a purely dissipative time-dependent Ginzburg-Landau equation. However, the effective action is nonlocal because of Landau damping. The relaxation is exponential in time with the damping rate $`\gamma _k(T)=14\zeta (3)v_\mathrm{F}k|\mathrm{\Delta }_0(T)|/3\pi ^3T`$. The factor $`v_\mathrm{F}k`$ is a consequence of the Goldstone nature of the phase fluctuations. The relaxation of phase fluctuations near the critical point features *critical slowing down*, i.e., the relaxation time scale diverges at the critical point.

Far away from the Ginzburg-Landau regime at low temperatures, the spectral density features sharp quasiparticle peaks and the nonequilibrium relaxation is *underdamped* in agreement with the results of Refs. aitchison95, ; aitchison97, . Away from the critical region, the contribution from Landau damping is negligible. The long-wavelength amplitude fluctuations are severely Landau-damped near the critical region and the relaxation time scales also feature critical slowing down.

While we have focused on the nonequilibrium dynamics of neutral BCS superconductors, as a next step we will apply these methods to the case of charged superconductors to study in detail the dynamics of the Carlson-Goldman modes as well as the effective action including gauge fields near the critical temperature. We postpone this study to a forthcoming article.

###### Acknowledgements.

We thank I. J. R. Aitchison and V. P. Gusynin for helpful and illuminating correspondence. S.M.A. would like to thank King Fahd University of Petroleum and Minerals for financial support. The work of D.B. was supported in part by the US NSF under grant PHY-9988720. The work of S.-Y.W. was supported in part by the US DOE under contract W-7405-ENG-36.

## Appendix A Plane wave solutions and the Bogoliubov transformation

In this Appendix we present an alternative derivation of the correlation functions for the fields $`\psi _\sigma `$, $`\psi _\sigma ^{}`$ directly from the plane wave solutions of the homogeneous equations of motion (i.e., in the absence of source) for the Nambu-Gor’kov field given by (31)

$$\left[i\frac{}{t}+\sigma _3\left(\frac{^2}{2m}+\mu \right)+\sigma _+\mathrm{\Delta }_0+\sigma _{}\mathrm{\Delta }_0^{}\right]\mathrm{\Psi }(𝐱,t)=0.$$
(142)
The plane wave solution can be written in the form

$$\mathrm{\Psi }(𝐱,t)=\mathrm{\Phi }(k)e^{i(\omega t𝐤𝐱)},\mathrm{\Phi }(k)=\left[\begin{array}{c}U_k\\ V_k\end{array}\right].$$
(143)
The two-component Nambu-Gor’kov spinor obeys

$$\left[\begin{array}{cc}\xi _k& \mathrm{\Delta }_0\\ \mathrm{\Delta }_0^{}& \xi _k\end{array}\right]\left[\begin{array}{c}U_k\\ V_k\end{array}\right]=\omega \left[\begin{array}{c}U_k\\ V_k\end{array}\right],$$
(144)
where $`\xi _k=k^2/2m\mu `$. This is an eigenvalue equation with the eigenvalues given by $`\omega =\pm E_k`$, where $`E_k=\sqrt{\xi _k^2+|\mathrm{\Delta }_0|^2}`$. The normalization of the positive and negative energy spinors is chosen so that $`\mathrm{\Phi }^{(\alpha )}(k)\mathrm{\Phi }^{(\beta )}(k)=\delta _{\alpha \beta }`$, where $`\alpha ,\beta =1,2`$ (not to be confused with the Nambu-Gor’kov indices) correspond to $`\omega =\pm E_k`$, respectively. Introducing the Bogoliubov coefficients $`u_k`$, $`v_k`$ given by

$$u_k=\left(\frac{\xi _k+E_k}{2E_k}\right)^{1/2},v_k=u_k\left(\frac{\mathrm{\Delta }_0^{}}{\xi _k+E_k}\right),u_kv_k=\frac{\mathrm{\Delta }_0^{}}{2E_k},$$
(145)
and satisfying $`u_k^2+|v_k|^2=1`$, we find that the positive and negative energy spinors are given by

$$\mathrm{\Phi }^{(1)}(k)=\left[\begin{array}{c}u_k\\ v_k\end{array}\right],\mathrm{\Phi }^{(2)}(k)=\left[\begin{array}{c}v_k^{}\\ u_k\end{array}\right],$$
(146)
respectively. After accounting for the interpretation of negative energy solutions as antiparticles, one can therefore write the general plane wave solution of the homogeneous equation of motion (142) as

$$\mathrm{\Psi }(𝐱,t)=\frac{1}{\sqrt{V}}\underset{𝐤}{}\left\{b_𝐤\left[\begin{array}{c}u_k\\ v_k\end{array}\right]e^{i(E_kt𝐤𝐱)}+d_𝐤^{}\left[\begin{array}{c}v_k^{}\\ u_k\end{array}\right]e^{i(E_kt𝐤𝐱)}\right\},$$
(147)
where $`V`$ is the quantization volume.

At the level of second quantization, one recognizes that (147) is the Bogoliubov transformation, in which the operator $`b_𝐤^{}`$, $`d_𝐤^{}`$ create Bogoliubov quasiparticles of momentum $`𝐤`$ (energy $`E_k`$) and obey the usual canonical anticommutation relations. The correlation functions of the Nambu-Gor’kov fields in the density matrix that describes free Bogoliubov quasiparticles in thermal equilibrium at inverse temperature $`\beta `$ are therefore found to be given by (in the continuum limit)

$$\mathrm{\Psi }_a(𝐱,t)\mathrm{\Psi }_b^{}(𝐱{}_{}{}^{},t^{})=\frac{d^3k}{(2\pi )^3}[[1n_\mathrm{F}(E_k)]𝒮_{ab}(k)e^{iE_k(tt^{})}+n_\mathrm{F}(E_k)\overline{𝒮}_{ab}(k)e^{iE_k(tt^{})}]e^{i𝐤(𝐱𝐱{}_{}{}^{})},$$

$$\mathrm{\Psi }_b^{}(𝐱{}_{}{}^{},t^{})\mathrm{\Psi }_a(𝐱,t)=\frac{d^3k}{(2\pi )^3}[n_\mathrm{F}(E_k)𝒮_{ab}(k)e^{iE_k(tt^{})}+[1n_\mathrm{F}(E_k)]\overline{𝒮}_{ab}(k)e^{iE_k(tt^{})}]e^{i𝐤(𝐱𝐱{}_{}{}^{})},$$

where $`n_\mathrm{F}(E_k)`$ is the equilibrium distribution for Bogoliubov quasiparticles of momentum $`𝐤`$

$$n_\mathrm{F}(E_k)=b_𝐤^{}b_𝐤=d_𝐤^{}d_𝐤=\frac{1}{e^{\beta E_k}+1},$$
(148)
and $`𝒮(k)`$, $`\overline{𝒮}(k)`$ are given by (45).