text:
# Hydrodynamic approach to time-dependent density functional theory; response properties of metal clusters

## I Introduction

The hydrodynamic analogy of quantum mechanics was first explored by Madelung  who transformed the single particle Schrodinger equation into a pair of hydrodynamical equations. The theory views  the electron cloud as a classical fluid moving under the action of classical Coulomb forces augmented by the forces of quantum origin. The basic dynamical variables of this theory are the particle and current densities which satisfy two fluid dynamical equations, namely, the continuity and an Euler type equation. The work of Madelung was followed by Bloch’s attempt  to develop hydrodynamical theory for many-electron systems within the realm of Thomas-Fermi (TF) theory . Although then proposed without any rigorous foundation, the theory can now be derived  from the equations of time-dependent density-functional theory (TDDFT) . It is based on the assumption that the dynamics of many-electron system can be described by considering it as a fluid of density $`\rho (𝐫,t)`$ and a velocity field $`𝐯(𝐫,t)`$ which is assumed to be curl free (that is $`𝐯(𝐫,t)=S(𝐫,t)`$, where $`S(𝐫,t)`$ is scalar velocity potential).

Using $`\rho (𝐫,t)`$ and $`S(𝐫,t)`$ as conjugate variables Bloch derived two fluid dynamical equations: the continuity equation

$$\frac{\rho }{t}+(\rho S)=0,$$
(1)
and the Euler equation

$$\frac{S}{t}=\frac{1}{2}|S|^2+\frac{\delta T_0}{\delta \rho }+v_{ext}(𝐫,t)+\frac{\rho (𝐫^{},t)}{|𝐫𝐫^{}|}𝑑𝐫^{}.$$
(2)
Here $`T_0`$ is the TF kinetic energy (KE) functional and $`v_{ext}(𝐫,t)`$ represents the external potential. These equations were subsequently used to study photoabsorption cross section and collective excitation of atoms, collective excitations  and plasmons of metal clusters  and surface plasmons  in metals.

Bloch’s theory also formed the basis of initial attempts by Ying  to extend density-functional theory (DFT) to include time-dependent (TD) external potentials. He did this by replacing the TF KE term $`T_0`$ by a general functional $`G[\rho ]`$ consisting of the KE and the exchange-correlation (XC) contribution to the total energy. In Ying’s work it is implicit that, like in the static DFT, a universal functional $`G[\rho (𝐫,t)]`$ can be written for the TD problem. The ad-hoc nature of Bloch’s theory and its extension by Ying was removed by the pioneering works of Deb and Ghosh , Bartolotti  and Runge and Gross . Runge and Gross rigorously proved the existence of a Hohenberg-Kohn  like theorem for TD potentials, and showed that the TD density $`\rho (𝐫,t)`$ can be determined by solving the hydrodynamical equations

$$\frac{\rho }{t}+𝐣=0,$$
(3)
which is the continuity equation, and the Euler’s equation

$$\frac{𝐣}{t}=𝐏[\rho (𝐫,t)].$$
(4)
Here $`𝐏`$ is the three-component density-functional of Runge and Gross and the vector $`𝐣`$ is the current density corresponding to the many-body wavefunction $`\mathrm{\Psi }(𝐫_1,\mathrm{},𝐫_N;t)`$. An explicit expression for $`𝐏[\rho (𝐫,t)]`$ in terms of the wavefunction has recently been given  using the TD differential virial theorem.

Although there have been some calculations, as mentioned above, in the past by employing hydrodynamic theory, its full potential remains unexplored. This is evidently because with the increasing computing resources one can perform  orbital based calculations like TD Hartee-Fock (HF) or TD Kohn-Sham (KS) with relative ease. Recently, however, hydrodynamical theory is being applied in situations where such orbital based calculations are still computationally difficult to implement. One such example is systems which contain thousands of atom such as nanoparticles and clusters. In these systems hydrodynamical theory becomes the method of choice. Thus the theory has been applied to study photoabsorption cross-section of metal particles , collective  and magnetoplasmon excitations  of confined electronic systems and also to study the interaction of strong laser light with atomic systems . Besides the computational ease offered by it, one is also tempted to work within the hydrodynamic formulation because it provides an intuitively appealing approach to the time dependent many-body problem.

In this paper we develop perturbation theory within the hydrodynamic formalism to calculate linear and nonlinear response properties of large systems. The motivation for this comes from our experience with the calculation of static response properties  employing density based perturbation theory  within the Hohenberg-Kohn formalism of DFT. The density based method reduces the numerical effort required for such calculations substantially while leading to reasonably accurate results  for the response properties. In the same manner, the hydrodynamical approach proves to be useful for calculating frequency dependent response properties of extended systems for which orbital based theories become quite difficult to implement because of the large number of orbitals involved.

The work presented in this paper is divided into two parts: First, we derive the generalized Bloch type equation using the concept of time-averaged energy (quasi-energy) of an electronic system subjected to a TD periodic field. For calculating optical response properties we then develop variation-perturbation (VP) theory in terms of the quasi-energy using particle and current densities as the basic variables. This is presented in section II. The perturbation theory developed here proceeds along the lines of density based stationary-state perturbation theory  by making use of the stationary nature of the time-averaged energy with respect to $`\rho (𝐫,t)`$ and $`S(𝐫,t)`$. In the second part we demonstrate the applicability of the perturbation theory developed here by calculating frequency dependent linear and nonlinear polarizabilities of inert gas atoms (section III) and comparing our numbers with the standard results obtained from the wavefunctional approach. Having demonstrated the accuracy of hydrodynamical approach, we then apply it to calculate frequency dependent response properties of metal clusters with number of atoms up to 5000 using the spherical jellium background model (SJBM).

## II Variation perturbation method in hydrodynamical theory

### A. Time-averaged energy

The central quantity around which the VP theory is developed in the static case is the ground-state energy of the system. For periodic TD hamiltonians, this role is played by the time-averaged energy or quasi-energy . Therefore, in the following we first derive an expression for the quasi-energy as a functional of particle and current densities and show that it obeys stationarity for the correct solutions of these functions. As expected from the work of Runge and Gross , stationarity with respect to the density leads to the equation of motion for the density. In addition, variation with respect to the current density gives the continuity equation.

We begin with the TD Schrodinger equation for the many-body wave function $`\mathrm{\Psi }(𝐫_1,\mathrm{},𝐫_N;t)`$ given by

$$\left(H(t)ı\frac{}{t}\right)\mathrm{\Psi }(𝐫,t)=0$$
(5)
where

$$H(t)=\widehat{T}+\widehat{V}_{ee}+\widehat{V}(t).$$
(6)
In the above equation $`\widehat{T}`$ and $`\widehat{V}_{ee}`$ are the kinetic and electron-electron interaction energy operators, respectively, and $`\widehat{V}(t)`$ denotes the TD external potential containing both the nuclear and the applied potential. For periodic hamiltonians, that is

$$H(t+T)=H(t),$$
(7)
where T is the time period, in accordance with Floquet’s theorem there exists a solution $`\mathrm{\Psi }(𝐫,t)`$ of the form

$$\mathrm{\Psi }(𝐫,t)=\mathrm{\Phi }(𝐫,t)e^{ıEt}$$
(8)
where $`\mathrm{\Phi }(𝐫_1,\mathrm{},𝐫_N;t)`$ is also periodic in time with time period T, i.e., $`\mathrm{\Phi }(t+T)=\mathrm{\Phi }(t)`$. Such a state has been termed as the steady-state of the system with $`E`$ being the corresponding quasi-energy. The equation of motion for $`\mathrm{\Phi }`$ is easily seen to be

$$\left(H(t)ı\frac{}{t}\right)\mathrm{\Phi }(𝐫,t)=E\mathrm{\Phi }(𝐫,t).$$
(9)
The corresponding expression for the quasi-energy is the time averaged expectation value

$$E[\mathrm{\Phi }]=\left\{\mathrm{\Phi }|H(t)ı\frac{}{t}|\mathrm{\Phi }\right\}.$$
(10)
The curly bracket in Eq.(10) denotes the time averaging over one period T defined as

$$\left\{fg\right\}=\frac{1}{T}_0^Tf^{}(t)g(t)𝑑t$$
(11)
The quasi-energy represents the average energy of induction  of a system subjected to a TD potential as is easily seen by the TD Hellmann-Feynman theorem .

To convert Eq.(10) into its hydrodynamical counterpart we decompose the complex steady-state many-body wavefunction in polar form, so that

$$\mathrm{\Phi }(𝐫,t)=\chi (𝐫,t)e^{ıS(𝐫,t)}$$
(12)
where both $`\chi `$ and $`S`$ are real functions of $`𝐫_1,𝐫_2,\mathrm{},𝐫_N`$ and are periodic in time with time period T. Further, $`S(𝐫,t)`$ also has a purely TD component which integrates to zero over time period T (for detail see Ref.). Note that, $`S`$ is zero for the ground-state of the system. By substituting Eq.(12) in Eq.(10), the expression for the average energy becomes

$$E[\chi ,S]=\left\{\chi |\widehat{T}^{}+\widehat{V}_{ee}+\frac{S}{t}|\chi ı\chi |\frac{\chi }{t}\right\}$$
(13)
where

$$\widehat{T}^{}=\underset{i}{}\left(\frac{1}{2}_i^2+\frac{1}{2}|_iS|^2\right)$$
(14)
Since the periodicity and reality of $`\chi (𝐫,t)`$ implies that

$`\left\{\chi |{\displaystyle \frac{\chi }{t}}\right\}`$ $`=`$ $`{\displaystyle \frac{1}{2T}}{\displaystyle _0^T}𝑑t{\displaystyle \frac{}{t}(\chi ^{}\chi )𝑑𝐫}`$ (15)
$`=`$ $`0,`$
the quasi-energy is given as

$$E[\chi ,S]=\left\{\chi |\underset{i}{}\frac{1}{2}_i^2+\widehat{V}_{ee}+\widehat{V}_{ext}+\underset{i}{}\frac{1}{2}(_iS)^2+\frac{S}{t}|\chi \right\}$$
(16)
Now by invoking the Runge-Gross theorem , it can be written as a functional of the density alone. However, in the hydrodynamical formulation the density and the velocity potential $`S`$ are treated as independent variables which means that the energy above is a functional of these two quantities. An advantage of this decoupling of the density and $`S`$ is that one does not have to know the functional dependence of the $`S`$ on the density. Further, this facilitates approximating the expectation value $`\chi |\frac{1}{2}^2|\chi `$ as a functional of the density by the KE functionals well studied in static DFT. Evidently, the first three terms of the equation above can be represented by a functional of TD density as

$$\left\{F[\rho (𝐫,t)]+\left[v_0(𝐫)+v_{app}(𝐫,t)\right]\rho (𝐫,t)𝑑𝐫\right\},$$
(17)
where $`v_0(𝐫)`$ represents the static external potential and TD part of the potential is represented by $`v_{app}(𝐫,t)`$. This is because changes in $`S`$ do not affect their values. The universal functional $`\left\{F[\rho (𝐫,t)]\right\}`$ given by

$$\left\{F[\rho (𝐫,t)]\right\}=\left\{T_s[\rho (𝐫,t)]\right\}+\left\{E_H[\rho (𝐫,t)]\right\}+\left\{E_{xc}[\rho (𝐫,t)]\right\}.$$
(18)
Here $`\left\{T_s[\rho (𝐫,t)]\right\}`$, $`\left\{E_H[\rho (𝐫,t)]\right\}`$ and $`\left\{E_{xc}[\rho (𝐫,t)]\right\}`$ represent the time-averaged KE, Hartree energy and exchange and correlation (XC) energy functionals respectively for the TD system. The TD particle density is given by

$$\rho (𝐫,t)=\chi ^{}(𝐫,𝐫_2,\mathrm{},𝐫_N,t)\chi (𝐫,𝐫_2,\mathrm{},𝐫_N,t)𝑑𝐫_2\mathrm{}𝑑𝐫_N$$
(19)
So far we have written the first three terms in terms of the particle density, a quantity defined in 3D configuration space. The last two terms representing the current still have all the co-ordinates of the configuration space in them. As such any equation involving $`S(𝐫_1,\mathrm{},𝐫_N;t)`$ can not be projected on to 3D space. To do this one needs to consider some approximate form for the phase $`S`$. One such approximation for $`S`$ which is generally employed  is that it can be written as the sum of single particle phases, that is

$$S(𝐫_1,\mathrm{},𝐫_N;t)=\underset{i=1}{\overset{N}{}}S(𝐫_i;t)$$
(20)
with the same function $`S`$ representing each electron. This approximation is equivalent to assuming the velocity field of the electron fluid to be curl free as was done by Bloch  in deriving Eq.(2). With this approximation the average energy functional of Eq.(16) is given as

$`E[\rho ,S]`$ $`=`$ $`\{F[\rho (𝐫,t)]+{\displaystyle }(v_0(𝐫)+v_{app}(𝐫,t))\rho (𝐫,t)d𝐫`$ (21)
$`+`$ $`{\displaystyle \frac{1}{2}}{\displaystyle }\rho (𝐫,t)(S)^2d𝐫+{\displaystyle }{\displaystyle \frac{S}{t}}\rho (𝐫,t)d𝐫\}`$
This is the expression for the quasi-energy of a many electron system (under the approximation made above) interacting with a TD periodic potential. Since the purely TD component of $`S(𝐫,t)`$ is the same as the phase of the wavefunction, it does not contribute to the energy. In Eq.21 this is ensured by $`\rho (𝐫,t)`$ integrating to a fixed number of electrons at all times. We therefore drop it and work with only the co-ordinate dependent component of $`S(𝐫,t)`$. This is similar to separating out the overall phase of TD wavefunction  in the TD perturbation theory.

We now demonstrate the variational nature of $`E[\rho ,S]`$ with respect to $`\rho `$ and $`S`$. Making $`E[\rho ,s]`$ stationary with respect to $`\rho `$ and $`S`$ gives the Euler equation

$$\mu (t)=\frac{S}{t}+\frac{1}{2}(S)^2+v_0(𝐫)+v_{app}(𝐫,t)+\frac{\delta F}{\delta \rho }$$
(22)
where $`\mu (t)`$ is the Lagrange-multiplier ensuring that $`\rho (𝐫,t)`$ integrates to the correct number of electrons at each instant of time, and the continuity equation

$$\frac{\rho }{t}+(\rho S)=0,$$
(23)
respectively. Eq.(22) is the same as that proposed by Ying. As such if $`F[\rho ]`$ is approximated by the TF functional, it gives the Bloch’s hydrodynamical equation correctly. Further, for time independent hamiltonians it correctly reduces to the Euler equation of static DFT. All these facts demonstrate the variational nature of $`E[\rho ,S]`$ with respect to $`\rho `$ and $`S`$. Employing this we now develop the VP method in terms of the particle and the current densities. We show that the $`(2n+1)`$ theorem and its variational corollary is satisfied in terms of these variables.

### B. Perturbation theory

To develop perturbation theory we assume that $`v_{ext}(𝐫,t)`$ is relatively weak and under its action the ground-state density $`\rho ^{(0)}(𝐫)`$ changes to $`\rho ^{(0)}(𝐫)+\mathrm{\Delta }\rho (𝐫,t)`$ and the velocity potential changes to $`S^{(0)}(𝐫)+\mathrm{\Delta }S(𝐫,t)`$. The particle density change $`\mathrm{\Delta }\rho `$ satisfies the normalization condition

$$\mathrm{\Delta }\rho (𝐫,t)𝑑𝐫=0.$$
(24)
However no such condition is required for the change in the velocity potential $`\mathrm{\Delta }S`$. The changes $`\mathrm{\Delta }\rho `$ and $`\mathrm{\Delta }S`$ are expanded in perturbation series as

$`\mathrm{\Delta }\rho `$ $`=`$ $`{\displaystyle \underset{j}{}}\rho ^{(j)}`$
$`\mathrm{\Delta }S`$ $`=`$ $`{\displaystyle \underset{j}{}}S^{(j)},`$ (25)
where $`\rho ^{(j)}`$ and $`S^{(j)}`$ correspond to the $`jth`$ order terms in the perturbation parameter. The energy corresponding to $`\rho ^{(0)}+\mathrm{\Delta }\rho `$ and $`S^{(0)}+\mathrm{\Delta }S`$ is given by

$`E[\rho ^{(0)}+\mathrm{\Delta }\rho ,S^{(0)}+\mathrm{\Delta }S]`$ $`=`$ $`\{F[\rho ^{(0)}+\mathrm{\Delta }\rho ]+{\displaystyle }(v_0+v_{app})(\rho ^{(0)}+\mathrm{\Delta }\rho )d𝐫`$ (26)
$`+`$ $`{\displaystyle \frac{1}{2}}{\displaystyle (S^{(0)}+\mathrm{\Delta }S)(S^{(0)}+\mathrm{\Delta }S)(\rho ^{(0)}+\mathrm{\Delta }\rho )𝑑𝐫}`$
$`+`$ $`{\displaystyle }{\displaystyle \frac{(S^{(0)}+\mathrm{\Delta }S)}{t}}(\rho ^{(0)}+\mathrm{\Delta }\rho )d𝐫\}`$
Using Eq.(26) we now obtain the energy changes to different orders in perturbation parameter employing an approach identical to the one adopted in Ref. for time independent density based perturbation theory. The resulting expressions for average energies to different orders are:

$$E^{(1)}=\left\{v_{app}^{(1)}(𝐫,t)\rho ^{(0)}(𝐫)𝑑𝐫\right\},$$
(27)
$`E^{(2)}`$ $`=`$ $`\{{\displaystyle \frac{1}{2}}{\displaystyle }{\displaystyle \frac{\delta ^2F[\rho ]}{\delta \rho (𝐫,t)\delta \rho (𝐫^{},t)}}\rho ^{(1)}(𝐫,t)\rho ^{(1)}(𝐫^{},t)d𝐫d𝐫^{}+{\displaystyle }v_{app}^{(1)}(𝐫,t)\rho ^{(1)}(𝐫)d𝐫`$ (28)
$`+`$ $`{\displaystyle }{\displaystyle \frac{S^{(1)}(𝐫,t)}{t}}\rho ^{(1)}(𝐫,t)d𝐫+{\displaystyle \frac{1}{2}}{\displaystyle }(S^{(1)}S^{(1)})\rho ^{(0)}(𝐫)d𝐫\},`$
$`E^{(3)}`$ $`=`$ $`\{{\displaystyle \frac{1}{6}}{\displaystyle }{\displaystyle \frac{\delta ^3F[\rho ]}{\delta \rho (𝐫,t)\delta \rho (𝐫^{},t)\delta \rho (𝐫^{\prime \prime },t)}}\rho ^{(1)}(𝐫,t)\rho ^{(1)}(𝐫^{},t)\rho ^{(1)}(𝐫^{\prime \prime },t)d𝐫d𝐫^{}d𝐫^{\prime \prime }`$ (29)
$`+`$ $`{\displaystyle }(S^{(1)}S^{(1)})\rho ^{(1)}(𝐫,t)d𝐫\},`$
$`E^{(4)}`$ $`=`$ $`\{{\displaystyle \frac{1}{2}}{\displaystyle }{\displaystyle \frac{\delta ^2F[\rho ]}{\delta \rho (𝐫,t)\delta \rho (𝐫^{},t)}}\rho ^{(2)}(𝐫,t)\rho ^{(2)}(𝐫^{},t)d𝐫d𝐫^{}`$ (30)
$`+`$ $`{\displaystyle \frac{1}{6}}{\displaystyle \frac{\delta ^3F[\rho ]}{\delta \rho (𝐫,t)\delta \rho (𝐫^{},t)\delta \rho (𝐫^{\prime \prime },t)}\rho ^{(2)}(𝐫,t)\rho ^{(1)}(𝐫^{},t)\rho ^{(1)}(𝐫^{\prime \prime },t)𝑑𝐫𝑑𝐫^{}𝑑𝐫^{\prime \prime }}`$
$`+`$ $`{\displaystyle \frac{1}{24}}{\displaystyle \frac{\delta ^4F[\rho ]}{\delta \rho (𝐫,t)\delta \rho (𝐫^{},t)\delta \rho (𝐫^{\prime \prime },t)\delta \rho (𝐫^{\prime \prime \prime },t)}\rho ^{(1)}(𝐫,t)\rho ^{(1)}(𝐫^{},t)\rho ^{(1)}(𝐫^{\prime \prime },t)\rho ^{(1)}(𝐫^{\prime \prime \prime },t)}`$
$`\times `$ $`d𝐫d𝐫^{}d𝐫^{\prime \prime }d𝐫^{\prime \prime \prime }`$
$`+`$ $`{\displaystyle \frac{1}{2}}{\displaystyle (S^{(2)}S^{(2)})\rho ^{(0)}(𝐫)𝑑𝐫}+{\displaystyle \frac{1}{2}}{\displaystyle (S^{(1)}S^{(1)})\rho ^{(2)}(𝐫,t)𝑑𝐫}`$
$`+`$ $`{\displaystyle }(S^{(1)}S^{(2)})\rho ^{(1)}(𝐫,t)d𝐫+{\displaystyle }{\displaystyle \frac{S^{(2)}}{t}}\rho ^{(2)}(𝐫,t)d𝐫\}`$
and

$`E^{(5)}`$ $`=`$ $`\{{\displaystyle \frac{1}{2}}{\displaystyle }{\displaystyle \frac{\delta ^3F[\rho ]}{\delta \rho (𝐫,t)\delta \rho (𝐫^{},t)\delta \rho (𝐫^{\prime \prime },t)}}\rho ^{(2)}(𝐫,t)\rho ^{(2)}(𝐫^{},t)\rho ^{(1)}(𝐫^{\prime \prime },t)d𝐫d𝐫^{}d𝐫^{\prime \prime }`$ (31)
$`+`$ $`{\displaystyle \frac{1}{24}}{\displaystyle \frac{\delta ^4F[\rho ]}{\delta \rho (𝐫,t)\delta \rho (𝐫^{},t)\delta \rho (𝐫^{\prime \prime },t)\delta \rho (𝐫^{\prime \prime \prime },t)}\rho ^{(1)}(𝐫,t)\rho ^{(1)}(𝐫^{},t)\rho ^{(1)}(𝐫^{\prime \prime },t)\rho ^{(2)}(𝐫^{\prime \prime \prime },t)}`$
$`\times d𝐫d𝐫^{}d𝐫^{\prime \prime }d𝐫^{\prime \prime \prime }`$
$`+`$ $`{\displaystyle \frac{1}{120}}{\displaystyle \frac{\delta ^5F[\rho ]}{\delta \rho (𝐫,t)\delta \rho (𝐫^{},t)\delta \rho (𝐫^{\prime \prime },t))\delta \rho (𝐫^{\prime \prime \prime },t)\delta \rho (𝐫^{\prime \prime \prime \prime },t)}\rho ^{(1)}(𝐫,t)\rho ^{(1)}(𝐫^{},t)\rho ^{(1)}(𝐫^{\prime \prime },t)}`$
$`\times \rho ^{(1)}(𝐫^{\prime \prime \prime },t)\rho ^{(1)}(𝐫^{\prime \prime \prime \prime },t)d𝐫d𝐫^{}d𝐫^{\prime \prime }d𝐫^{\prime \prime \prime }d𝐫^{\prime \prime \prime \prime }`$
$`+`$ $`{\displaystyle }(S^{(2)}S^{(2)})\rho ^{(1)}(𝐫,t)d𝐫+{\displaystyle }(S^{(1)}S^{(2)})\rho ^{(2)}(𝐫,t)d𝐫\}`$
In deriving these equations (Eq.(27)-(31)) we have made use of the fact that for the ground-state $`S^{(0)}=0`$ which implies that the current density $`\rho ^{(0)}S^{(0)}`$ and the time-derivative $`\frac{S^{(0)}}{t}`$ vanish for the ground-state. In addition we also use the first-order

$$\frac{\rho ^{(1)}}{t}+(\rho ^{(0)}S^{(1)})=0$$
(32)
$$\mu ^{(1)}(t)=\frac{S^{(1)}}{t}+v_{app}^{(1)}(𝐫,t)+\frac{1}{2}\frac{\delta ^2F[\rho ]}{\delta \rho (𝐫,t)\delta \rho (𝐫^{},t)}\rho ^{(1)}(𝐫^{},t)𝑑𝐫^{}$$
(33)
and the second-order

$$\frac{\rho ^{(2)}}{t}+(\rho ^{(0)}S^{(2)})+(\rho ^{(1)}S^{(1)})=0$$
(34)
$`\mu ^{(2)}(t)`$ $`=`$ $`{\displaystyle \frac{S^{(2)}}{t}}+{\displaystyle \frac{1}{2}}(S^{(1)}S^{(1)})`$ (35)
$`+`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{\delta ^2F[\rho ]}{\delta \rho (𝐫,t)\delta \rho (𝐫^{},t)}\rho ^{(2)}(𝐫^{},t)𝑑𝐫^{}}`$
$`+`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{\delta ^3F[\rho ]}{\delta \rho (𝐫,t)\delta \rho (𝐫^{},t)\delta \rho )𝐫^{\prime \prime },t)}\rho ^{(1)}(𝐫^{},t)\rho ^{(1)}(𝐫^{\prime \prime },t)𝑑𝐫^{}𝑑𝐫^{\prime \prime }}`$
continuity and Euler equations obtained by expanding Eqs.(21) and (22). Expressions for the average energies (Eqs.(27)-(31)) clearly demonstrates the $`(2n+1)`$ rule of perturbation theory. Energies up to order 3 are determined completely by $`\rho ^{(1)}`$ and $`S^{(1)}`$. Similarly $`\rho ^{(2)}`$ and $`S^{(2)}`$ give energy up to the fifth-order. This is the $`(2n+1)`$ theorem of hydrodynamic perturbation theory in terms of the particle and the current densities. Moreover, even-order corollary of this theorem also holds true. Thus making E<sup>(2)</sup> stationary with respect to $`S^{(1)}`$ and $`\rho ^{(1)}`$, leads to Eqs.(32) and (33). Similarly stationarity of E<sup>(4)</sup> with respect to $`\rho ^{(2)}`$ and $`S^{(2)}`$ gives the correct perturbation equations (Eqs.(34) and (35)) for $`\rho ^{(2)}`$ and $`S^{(2)}`$. The stationary nature of the even-order energies gives a variational method to obtain approximate solutions for the corresponding induced densities and currents.

With this we complete the development of VP method in terms of particle and current densities in hydrodynamic formulation of TDDFT. In the next section we demonstrate the applicability of this theory by calculating the frequency dependent linear and nonlinear polarizabilities of inert gas atoms and comparing the results obtained with their wavefunctional counterparts. We then apply the formalism to calculate frequency dependent polarizabilities and plasmon frequencies of alkali metal clusters of large sizes.

## III. Application of hydrodynamical formalism

To demonstrate the applicability of the formalism developed above we begin this section with the calculation of frequency dependent linear and nonlinear polarizabilities of inert gas atoms. In the present formulation we can calculate the nonlinear polarizabilities corresponding to the degenerate four wave mixing (DFWM) and DC-Kerr  effect which are directly related to the fourth-order energy changes. On the other hand, unlike the orbital based Kohn-Sham approach, $`(2n+1)`$ theorem cannot be exploited to calculate  the coefficients for the third-harmonic generation and electric field induced second harmonic generation processes from only the second-order induced densities. In the following we will present the results for the nonlinear coefficients corresponding to the DFWM process only.

As pointed out earlier, we perform our calculations using the variational property of the even-order energies. To this end we choose an appropriate variational form for the induced particle and current densities. For an applied potential of the form

$$v_{app}^{(1)}(𝐫,t)=v_{app}^{(1)}(𝐫)\mathrm{cos}\omega t$$
(36)
with the spatial part given by

$$v_{app}^{(1)}(𝐫)=r\mathrm{cos}\theta $$
(37)
where $``$ is amplitude of the applied field, the time dependence of $`\rho `$ and $`S`$ at various orders can easily be inferred from Eqs(32)-(34) as

$`\rho ^{(1)}(𝐫,t)`$ $`=`$ $`\rho ^{(1)}(𝐫)\mathrm{cos}\omega t`$
$`\rho ^{(2)}(𝐫,t)`$ $`=`$ $`\rho _2^{(2)}(𝐫)\mathrm{cos}2\omega t+\rho _0^{(2)}(𝐫)`$ (38)
and

$`S^{(1)}(𝐫,t)`$ $`=`$ $`S^{(1)}(𝐫)\mathrm{sin}\omega t`$
$`S^{(2)}(𝐫,t)`$ $`=`$ $`S^{(2)}(𝐫)\mathrm{sin}2\omega t.`$ (39)
Note that unlike the second-order particle density, the corresponding current has no constant term. This is consistent with the fact that the current arises due to the flow of electrons which causes the density to be time dependent. The spatial part of the induced particle and current densities are determined variationally by minimizing the appropriate even-order energies. For this purpose we choose the forms of the induced particle densities similar to the ones used previously  for the calculation of static response properties. These are

$$\rho ^{(1)}(𝐫)=\mathrm{\Delta }_1(r)\rho ^{(0)}(𝐫)\mathrm{cos}\theta $$
(40)
and

$`\rho _2^{(2)}(𝐫)`$ $`=`$ $`\left(\mathrm{\Delta }_2^2(r)+\mathrm{\Delta }_3^2(r)\mathrm{cos}^2\theta \right)\rho ^{(0)}(𝐫)+\lambda _2\rho ^{(0)}(𝐫)`$
$`\rho _0^{(2)}(𝐫)`$ $`=`$ $`\left(\mathrm{\Delta }_2^0(r)+\mathrm{\Delta }_3^0(r)\mathrm{cos}^2\theta \right)\rho ^{(0)}(𝐫)+\lambda _0\rho ^{(0)}(𝐫)`$ (41)
with

$$\mathrm{\Delta }_i(r)=\underset{j}{}a_j^ir^j$$
(42)
where $`\rho ^{(0)}(𝐫)`$ is the ground-state density, $`a_j^i`$ are the variational parameters and $`\lambda `$s are fixed by the normalization condition for $`\rho ^{(2)}`$. To ensure satisfaction of the normalization condition at all times, $`\rho _2^{(2)}`$ and $`\rho _0^{(2)}`$ are each normalized separately. These forms for the induced particle densities are motivated by the exact solutions  for the hydrogen atom in a static field and have been shown  to lead to accurate static polarizabilities. On the basis of the continuity equation at each order and Eqs.(40)-(42), we choose

$$S^{(1)}(𝐫)=\mathrm{\Delta }_1^s(r)\mathrm{cos}\theta $$
(43)
and

$$S^{(2)}(𝐫)=\left(\mathrm{\Delta }_2^s(r)+\mathrm{\Delta }_3^s(r)\mathrm{cos}^2\theta \right),$$
(44)
where

$$\mathrm{\Delta }_i^s(r)=\underset{j}{}b_j^ir^j$$
(45)
with $`b_j^i`$ being the variational parameters to be determined by minimizing the average energy of respective orders.

Application of hydrodynamical equations also requires approximating the KE and XC energy functionals. Based on our experience with the calculation of static linear and nonlinear polarizabilities we approximate them by their static forms. Thus for the KE, we use the von Weizsacker  functional.

$$T_s[\rho ]=\frac{1}{8}\frac{\rho \rho }{\rho }𝑑𝐫$$
(46)
For a discussion on the rationale behind choosing this functional, we refer the reader to the literature . For the exchange energy functional we use the adiabatic local-density approximation (ALDA) given by the Dirac exchange functional

$`E_x[\rho ]`$ $`=`$ $`C_x{\displaystyle \rho ^{\frac{4}{3}}(𝐫)𝑑𝐫}`$
$`C_x`$ $`=`$ $`{\displaystyle \frac{3}{4}}\left({\displaystyle \frac{3}{\pi }}\right)^{\frac{1}{3}}.`$ (47)
The correlation energy functional within the ALDA is represented by the Gunnarsson-Lundquist (GL) parametrization . Thus

$$E_c[\rho ]=ϵ_c(\rho )\rho (𝐫)𝑑𝐫$$
(48)
with

$$ϵ_c(\rho )=0.0333\left[(1+x^3)\mathrm{ln}(1+\frac{1}{x})+\frac{1}{2}xx^2\frac{1}{3}\right]$$
(49)
where

$$x=\frac{r_s}{11.4}$$
(50)
and $`r_s=(\frac{3}{4\pi }\frac{1}{\rho })^{\frac{1}{3}}`$ measures the radius in atomic units of a sphere which encloses one electron.

The theory presented here treats the non-interacting KE exactly for single orbital systems (hydrogen and helium atoms). We have checked this by calculating the linear and nonlinear polarizabilities of H and He atoms. They match well with the corresponding wavefunctional results. The real test of the theory is therefore when it is applied to systems with more than one orbitals. We now present the results of these calculations by first discussing the frequency dependent polarizability numbers for the inert gas atoms of neon and argon. The ground-state electronic densities of these atoms are obtained by employing the van Leeuwen and Baerends (LB)  potential. We use this potential as the orbitals generated by it have the correct asymptotic nature so that they lead to accurate values for the static  and frequency dependent  response properties. Here, instead of using the orbitals, we are using the ground-state densities generated by this potential.

In Figs. 1 and 2 we show the linear polarizabilities $`\frac{\alpha (\omega )}{\alpha (0)}`$ for neon and argon atoms, respectively, as a function of $`\omega `$, obtained by hydrodynamic approach. We compare these results (represented by open squares) with those obtained  within the Kohn-Sham formalism shown in the figures by filled squares. As is evident, the frequency dependence matches quite well with the KS approach. This along with the zero frequency results demonstrate that the hydrodynamic theory is capable of giving reasonable estimate of dynamic polarizabilities in the optical range. Notice though that in the present approach the increase of $`\alpha (\omega )`$ with respect to $`\omega `$ is slightly less.

To further quantify our results, we have fitted the frequency dependent polarizabilities with the formulae

$$\alpha (\omega )=\alpha (0)\left(1+C_2\omega ^2\right)$$
(51)
for small frequencies (up to $`\omega =0.05`$ a.u.). In Table I we give the results for $`C_2`$ obtained from both the hydrodynamic and the orbital based calculations . As anticipated from the discussion above, the values of C<sub>2</sub> obtained from hydrodynamic formulation are close to but slightly smaller than their wavefunctional counterparts.

Next, to study the performance of hydrodynamic approach in calculation of nonlinear response properties we calculate the coefficient corresponding to the DFWM phenomenon. These results are presented in Table II for two different frequencies, $`\lambda =10550`$Å ($`\omega 0.0433`$ a.u.) and $`\lambda =6943`$Å ($`\omega 0.0657`$ a.u.). Here also we compare the present results with the corresponding numbers obtained  by the TD Kohn-Sham method. From this Table we again observe that the results for DFWM coefficients at both the frequencies are also quite close to, but lower than, the corresponding wavefunctional numbers. Notice that the numbers for hyperpolarizabilities are also underestimated by the hydrodynamic approach and the maximum deviation from the TD Kohn-Sham number is about 10%. However, with the increase in frequency the difference between the hydrodynamical and the wavefunctional results will get larger. Nonetheless, the results obtained are reasonably accurate for frequencies up to $`\omega 0.05`$ a.u.. This is quite encouraging since the experimental measurement of nonlinear coefficients fall within this frequency regime and also the numerical effort required for the hydrodynamical calculation is substantially less in comparison to the wavefunctional approach.

Motivated by these results, in the next section we apply the hyrodynamic approach to calculate frequency dependent response properties and plasmon frequencies of metal clusters. As is well known, orbital based calculation for such systems are computationally demanding because of the large number of orbitals involved as the cluster size grows.

### Response properties of metal clusters

The hydrodynamic approach developed above is particularly useful for systems where an orbital based theory cannot be applied. Clusters are one such class of systems. These are made up of tens to thousands of atoms with properties distinct from the bulk properties of the constituent material. Further, various properties of clusters evolve in a well defined manner as their size grows. This was demonstrated by Knight et al. in their study of alkali metal clusters. Since then metal clusters have been studied  quite extensively. One of the simplest model that describes average properties of these systems correctly is the spherical jellium background model (SJBM) . In small size clusters, Kohn-Sham LDA equations can be easily solved within this model. On the other hand, for large clusters one switches over to the density based theories ; the mostly applied one has been the extended Thomas-Fermi (ETF) theory . In this paper also we use the density obtained by the ETF to calculate frequency dependent response properties (both linear and nonlinear) of alkali metal clusters with the number of atoms up to 5000. In the past, dynamic linear polarizabilities of these clusters have been studied using the time-dependent Kohn-Sham theory  but because of the difficulties mentioned above, the size up to which one could go has been limited.

In the ETF method the ground-state density is obtained by minimizing the energy functional

$$E[\rho ]=T_s^{ETF}[\rho ]+E_H[\rho ]+E_{xc}^{LDA}[\rho ]+V_I(𝐫)\rho (𝐫)𝑑𝐫+E_I$$
(52)
where $`V_I`$ and $`E_I`$ are the potential and the total electrostatic energy, respectively, of the ionic background. The functional $`T_s^{ETF}`$ is the non-interacting KE functional included up to the fourth-order in density gradients. It is given as

$$T_s^{ETF}[\rho ]=T_s^{(0)}[\rho ]+T_s^{(2)}[\rho ]+T_s^{(4)}[\rho ]$$
(53)
where

$`T_s^{(0)}`$ $`=`$ $`(3\pi ^2)^{2/3}{\displaystyle \rho ^{\frac{5}{3}}𝑑𝐫}`$
$`T_s^{(2)}`$ $`=`$ $`{\displaystyle \frac{1}{72}}{\displaystyle \frac{|\rho |^2}{\rho }𝑑𝐫}`$
$`T_s^{(4)}`$ $`=`$ $`{\displaystyle \frac{1}{540(3\pi ^2)^{\frac{3}{2}}}}{\displaystyle \rho ^{\frac{1}{3}}\left[\left(\frac{^2\rho }{\rho }\right)^2\frac{9}{8}\left(\frac{^2\rho }{\rho }\right)\left(\frac{\rho }{\rho }\right)^2+\frac{1}{3}\left(\frac{\rho }{\rho }\right)^4\right]𝑑𝐫},`$ (54)
E$`{}_{H}{}^{}[\rho ]`$ is the Hartree energy and for E$`{}_{}{}^{LDA}{}_{xc}{}^{}`$ is the LDA XC energy. For this we use the GL parametrization . The energy above is minimized by taking the variational form  for the density to be

$$\rho (r)=\rho _0\left[1+exp\left(\frac{rR}{\alpha }\right)\right]^\gamma $$
(55)
where $`R`$, $`\alpha `$ and $`\gamma `$ are the variational parameters and $`\rho _0`$ is fixed by the normalization condition for each set of these parameters. The density so obtained gives results which are quite close  to the results of more accurate KS calculations of several properties. We use the ground-state densities obtained by this method as the input for the calculation of response properties. We perform our calculations for sodium clusters with $`r_s=4.0`$, where $`r_s`$ is Wigner-Seitz radius of metal. Before presenting our results we point out that that the KE functional used to obtain the ground-state density and that for calculating the response properties are different. This is because whereas ETF functional is good for the total energies, it does not give the changes in the energies accurately. As mentioned earlier, for this purpose we use the von Weizsacker  functional.

First we present the results for static linear polarizabilities. Although polarizability of metal clusters has been investigated extensively in the past , these studies have been restricted to clusters with number of atoms up to 200 because of the use of the orbitals in the calculations. We perform our study for clusters up to 5000 atoms. The variational forms for the induced particle and the current densities are chosen to be similar to those used in the atomic case. In Fig.3 we show plot of static polarizability $`\alpha (0)`$ in the units of $`R_0^3`$ as a function of $`R_0`$ (where $`R_0=r_sN^{\frac{1}{3}}`$, denotes the radius of cluster). It clearly shows that results of our calculation match quite well with the results obtained by Kohn-Sham approach for small (up to 196 atoms) clusters. As the size of cluster grows the polarizability approaches the classical limit of $`R_0^3`$ (that is $`\frac{\alpha (0)}{R_0^3}1`$). This is exhibited rather clearly in Fig.3.

Having obtained the static polarizabilities of alkali metal clusters accurately, we next study the dynamic response properties of metal clusters focussing our attention particularly on the dipole resonance. The classical theory of dynamic polarizability predicts a single dipole resonance at the frequency given by  (in a.u.)

$$\omega _{Mie}=\sqrt{\frac{1}{r_s^3}},$$
(56)
which is equal to $`1/\sqrt{3}`$ times the bulk plasma frequency. The TDDFT results for the dipole resonance follow  the Mie result only in a qualitative way. The resonance peak corresponding to the photo absorption spectra of these clusters exhaust about 70%-90% of the dipole sum rule and is red-shifted by about 10%-20% from the classical Mie formula .

In our work we estimate these resonance peaks from the frequency dependent polarizabilities by approximately locating the frequency at which $`\alpha (\omega )`$ becomes very large. These results are presented in Table III along with the results obtained by Brack  using the RPA sum rules. It is clear from Table III that the dipole resonance frequencies of metal clusters obtained by hydrodynamical approach to TDDFT are quite accurate over the range of clusters studied. Further, the accuracy is better for larger clusters. We also find that with the increase in particle size the dipole resonance frequency approaches the classical Mie resonance frequency, which in the present case of $`r_s=4.0`$ is 0.125 a.u..

Next we discuss the results obtained for the nonlinear polarizabilities of metal clusters by the present approach. To the best of our knowledge hyperpolarizabilities for these systems have not been calculated before the present work. In these calculations we are restricted to clusters with maximum of only 300 atoms due to computational difficulties. Since electrons in metallic cluster are highly delocalized, we expect that the nonlinear response of these system should be quite large and increase rapidly with the size of the clusters. To ascertain how does the static hyperpolarizability $`\gamma (0)`$ scale with the size of clusters, in Fig.4 we plot $`\gamma (0)`$ versus $`\alpha (0)`$ on a log-log scale. The line in this figure represents the best fit to hyperpolarizability versus polarizability numbers obtained by us. It is seen from figure that for the clusters studied by us, the hyperpolarizability is linearly proportional to the linear polarizability. Since $`\alpha (0)`$ scales linearly with N, $`\gamma (0)`$ also varies in the same manner. This is a surprising result since in the atomic case we have seen that $`\gamma (0)`$ increases much more rapidly than $`\alpha (0)`$ does. This could be because the electrons in metal clusters are much more mobile and therefore screen the applied potential very efficiently.

We have also studied the frequency dependent hyperpolarizabilities $`\gamma (\omega ;\omega ,\omega ,\omega )`$ and found it increasing with $`\omega `$. Variation of $`\gamma (\omega ;\omega ,\omega ,\omega )`$ with $`\omega `$ is shown in Fig. 5. It becomes quite large (by an order of magnitude in comparison to the static result) at approximately half the dipolar resonance frequency obtained from $`\alpha (\omega )`$ (Table III). This demonstrates the inherent consistency of the theory.

Our study above has been done for clusters with number of atoms up to 300. However, the trends obtained should continue as the size grows. Slow increase of $`\gamma `$ with N is consistent with the fact that the classically $`\gamma `$ for spherical metal particle is zero. One reason why computation becomes difficult for large clusters, we suspect, is that the variational forms chosen for the second-order particle and current densities may not be appropriate for very large clusters. As such investigations for larger clusters relegated to the future studies.

## Concluding remarks

In this paper we have developed the time-dependent perturbation theory for periodic (in time) hamiltonian in terms of the particle and current densities of electrons. For this we have employed the hydrodynamic equations of TDDFT. Application of the theory developed requires that the energy functional be approximated. We have demonstrated that with the von Weizsacker  functional for the KE and the ALDA for the XC energy, the theory leads to reasonably accurate results for dynamic response properties, both linear and nonlinear of atoms. Having established that we have applied the theory to study response properties of metallic clusters with number of atoms up to 5000 within the SJBM. Of particular interest is how the hyperpolarizability varies with the size of these clusters. Although it is zero classically, our study shows that it increases linearly with the number of atoms in the cluster.

Acknowledgement: We thank Prof. M. Brack for sending us his program to calculate ground-state density using the ETF approach.

### Table Captions

Table I $`C_2`$ for inert gas atoms obtained by using hydrodynamical and wavefunctional approaches.

Table II DFWM coefficient $`\gamma (\omega ;\omega ,\omega ,\omega )`$ (in atomic units) for inert gas atoms using hydrodynamic approach.

Table III Estimate of dipole resonance frequencies (in atomic units) of some metal clusters by using hydrodynamic approach.

## Table I

| Atoms | $`C_2`$ | $`C_2`$ |
| --- | --- | --- |
|  | (hydrodynamical) | (wavefunctional)<sup>(a)</sup> |
| He | 1.12 | 1.12 |
| Ne | 0.82 | 1.04 |
| Ar | 2.16 | 2.65 |
| Kr | 2.79 | - |

(a) Ref.

## Table II

| Atoms | $`\lambda =6943\AA `$ | | $`\lambda =10550\AA `$ | |
| --- | --- | --- | --- | --- |
|  | hydrodynamic | wavefunctional<sup>(a)</sup> | hydrodynamic | wavefunctional<sup>(a)</sup> |
| He | 44.47 | 44.57 | 43.50 | 43.58 |
| Ne | 83.38 | 94.06 | 82.58 | 91.65 |
| Ar | 1095.2 | 1226.1 | 1039.3 | 1154.8 |
| Kr | 2392.7 | - | 2229.5 | - |
| Xe | 5821.1 | - | 5295.3 | - |

(a) Ref.

## Table III

| N | Present | RPA<sup>(a)</sup> |
| --- | --- | --- |
| 8 | 0.100 | 0.113 |
| 100 | 0.105 | 0.1198 |
| 500 | 0.1165 | 0.1219 |
| 1000 | 0.121 | 0.1226 |
| 5000 | 0.1223 | 0.1236 |

(a) Ref.

### Figure Captions

Fig.1 Plot of $`\alpha (\omega )/\alpha (0)`$ as a function of frequency $`\omega `$ for neon. The open and closed squares represent hydrodynamical and wavefunctional  results respectively.

Fig.2 Plot of $`\alpha (\omega )/\alpha (0)`$ as a function of frequency $`\omega `$ for argon. The open and closed squares represent hydrodynamical and wavefunctional  results respectively.

Fig.3 Static polarizability $`\alpha (0)`$ in the units of $`R_0^3`$ of alkali metal clusters as a function of $`R_0`$. The filled squares represents the results of Kohn-Sham calculations .

Fig.4 Plot of $`\gamma (0)`$ versus $`\alpha (0)`$ of alkali metal clusters.

Fig.5 Plot of $`\gamma (\omega ;\omega ,\omega ,\omega )`$ in the units of $`R_0^3`$ as a function of frequency $`\omega `$ for a metal cluster with 100 atoms.


format: text

text:
# Deriving effective transport equations for non-Abelian plasmas

## 1 Introduction

Understanding the physics of a quark-gluon plasma at high temperature or density – which is expected to be detected in the up-coming heavy ion experiments at RHIC and LHC – requires reliable theoretical tools to describe transport phenomena in non-Abelian plasmas in- or out-of-equilibrium. Such methods are also of relevance for the physics of the early universe, if the baryon asymmetry could finally be understood within an electroweak framework.

In this talkbbbInvited talk given by DFL at 5th workshop on QCD, Villefranche-sur-Mer, 3-7 Jan 2000. an approach is discussed which allows for a systematic derivation of effective kinetic equations for non-Abelian plasmas.$`^\mathrm{?}`$$`^\text{-}`$$`^\mathrm{?}`$ The interest in developing a formalism based on a classical transport theory for non-Abelian gauge fields $`^\mathrm{?}`$ is that main properties of a hot quantum plasma can already be understood within simple classical terms. Indeed, the soft non-Abelian gauge fields – those having a huge occupation number – can be treated as classical fields, while the hard gauge field modes can be treated as classical particles. This method, based on ‘integrating-out’ fluctuations about some mean fields, has been known since long for Abelian plasmas $`^\mathrm{?}`$, but a consistent extension to the non-Abelian case was longtime missing.

We apply our formalism to a hot plasma close to equilibrium. While the HTL effective theory is recovered in the Vlasov approximation $`^\mathrm{?}`$, the inclusion of leading order corrections due to fluctuations results in Bödeker’s effective theory.$`^\mathrm{?}`$ For a more extensive discussion of this and related effective theories we refer to the discussion in ref. $`^\mathrm{?}`$. Some further approaches to study quark-gluon dynamics are considered in ref. $`^\mathrm{?}`$.

## 2 Microscopic transport equations

The starting point for a classical transport theory is to consider an ensemble of classical point particles carrying a non-Abelian charge $`Q_a`$, where the colour index runs from $`a=1`$ to $`N^21`$ for a SU($`N`$) gauge group. These particles interact self-consistently amongst each others, that is, through the classical gauge fields created by the particles. Their classical equations of motion are known as the Wong equations $`^\mathrm{?}`$,

$$m\frac{d\widehat{x}^\mu }{d\tau }=\widehat{p}^\mu ,m\frac{d\widehat{p}^\mu }{d\tau }=g\widehat{Q}^aF_a^{\mu \nu }\widehat{p}_\nu ,m\frac{d\widehat{Q}^a}{d\tau }=gf^{abc}\widehat{p}^\mu A_\mu ^b\widehat{Q}^c.$$
(1)
Here, $`A_\mu `$ denotes the gauge field, $`F_{\mu \nu }^a=_\mu A_\mu ^a_\nu A_\mu ^a+gf^{abc}A_\mu ^bA_\nu ^c`$ the corresponding field strength and $`\widehat{z}(\widehat{x},\widehat{p},\widehat{Q})(\tau )`$ the world line of a particle. An ensemble of such particles can be described by their one-particle distribution function $`f_i𝑑\tau \delta (z\widehat{z}_i)`$. Making use of (1), its kinetic equation is $`^\mathrm{?}`$

$$p^\mu \left(\frac{}{x^\mu }gf^{abc}A_\mu ^bQ^c\frac{}{Q^a}gQ_aF_{\mu \nu }^a\frac{}{p_\nu }\right)f(x,p,Q)=0.$$
(2)
This Boltzmann equation is collisionless, since $`df/d\tau =0`$ (Liouville’s theorem). However, it contains effectively collisions due to the long range interactions as shall become clear in the sequel. Eq. (2) is completed by the Yang-Mills equation $`D_\mu F^{\mu \nu }=J^\nu `$, with the current $`J[f]`$ obtained self-consistently from the particles. In this classical picture, $`f`$ is a deterministic quantity once all initial conditions have been fixed. This is, however, not possible, and $`f`$ has to be considered as a stochastic (fluctuating) quantity instead.

## 3 Coarse graining

Of particular interest for most physical applications are the macroscopic properties of such a plasma of particles. In this case, a more appropriate quantity to study would be a coarse-grained distribution function $`\overline{f}`$, whose dynamical equation is expected to be of the Boltzmann-Langevin type,

$$p^\mu \left(\overline{D}_\mu gQ_a\overline{F}_{\mu \nu }^a_p^\nu \right)\overline{f}=C[\overline{f}]+\zeta .$$
(3)
(We introduced the shorthands $`D_\mu [A]f[_\mu gf^{abc}Q_cA_{\mu ,b}_a^Q]f`$, $`_\mu /x^\mu `$, $`_\mu ^p/p^\mu `$ and $`_a^Q/Q^a`$). A coarse-grained distribution function obtains from the microscopic one after a suitable smearing-out over some characteristic volume, which, when applied to (2), results in a collision integral $`C[f]`$ and a related source for stochastic noise $`\zeta `$. Physically speaking, these terms are expected because the coarse-graining integrates-out the modes within a coarse-graining volume. The interactions of such modes can result into additional effective interactions for the remaining long range modes. Also, the stochastic fluctuations of $`f`$ are smeared-out, resulting in a noise term $`\zeta `$.

## 4 Effective transport equations

An effective kinetic equation, like (3), can be derived systematically within classical transport theory, starting with (2). To that end, we perform an ensemble average in phase space which allows to split $`f`$ and the gauge fields into their mean field part and a fluctuating part

$`f(x,p,Q)`$ $`=`$ $`\overline{f}(x,p,Q)+\delta f(x,p,Q),`$ (4)
$`A^\mu (x)`$ $`=`$ $`\overline{A}^\mu (x)+a^\mu (x),`$ (5)
where $`\overline{f}=f`$ and $`\overline{A}=A`$. The mean value of the statistical fluctuations vanish by definition, $`\delta f=0`$ and $`a=0`$. This separation into mean fields and statistical, random fluctuations corresponds effectively to a split into long/short wavelength modes associated to the mean fields/fluctuations. Performing the ensemble average of (2) and the Yang-Mills equations gives $`^{\mathrm{?},\mathrm{?}}`$

$`p^\mu \left(\overline{D}_\mu gQ_a\overline{F}_{\mu \nu }^a_p^\nu \right)f`$ $`=`$ $`\eta +\xi .`$ (6)
$`\overline{D}_\mu \overline{F}^{\mu \nu }+J_{\text{fluc}}^\nu `$ $`=`$ $`\overline{J}^\nu .`$ (7)
In (6) and (7), we collected all terms quadratic or cubic in the fluctuations into

$`\eta `$ $``$ $`gQ_ap^\mu \left[(\overline{D}_\mu a_\nu \overline{D}_\nu a_\mu )^a+gf^{abc}a_\mu ^ba_\nu ^c\right]_p^\nu \delta f(x,p,Q),`$ (8)
$`\xi `$ $``$ $`gp^\mu f^{abc}Q^c\left[_a^Qa_\mu ^b\delta f(x,p,Q)+ga_\mu ^aa_\nu ^b_p^\nu \overline{f}(x,p,Q)\right],`$ (9)
$`J_{\text{fluc}}^{a,\nu }`$ $``$ $`g\left[f^{dbc}\overline{D}_{ad}^\mu a_{b,\mu }a_c^\nu +f^{abc}a_{b,\mu }\left((\overline{D}^\mu a^\nu \overline{D}^\nu a^\mu )_c+gf^{cde}a_d^\mu a_e^\nu \right)\right].`$ (10)
In addition, we need the dynamical equations for the fluctuations. They are derived from subtracting (2) from (6) as

$`p^\mu \left(\overline{D}_\mu gQ_a\overline{F}_{\mu \nu }^a_p^\nu \right)\delta f`$ $`=`$ $`gQ_a(\overline{D}_\mu a_\nu \overline{D}_\nu a_\mu )^ap^\mu _\nu ^p\overline{f}`$ (11)
$`+gp^\mu a_{b,\mu }f^{abc}Q_c_a^Q\overline{f}+\eta +\xi \eta +\xi `$
$`\left(\overline{D}^2a^\mu \overline{D}^\mu (\overline{D}_\nu a^\nu )\right)^a`$ $`+`$ $`2gf^{abc}\overline{F}_b^{\mu \nu }a_{c,\nu }+J_{\text{fluc}}^{a,\mu }J_{\text{fluc}}^{a,\mu }=\delta J^{a,\mu }.`$ (12)
Eqs. (11) and (12) are the master equations for all higher order correlation functions of fluctuations. (This hierarchy of dynamical equations is similar to the BBGKY hierarchy.) All the equal-time correlators of $`\delta f`$ can be derived from the basic definition of the Gibbs ensemble average in phase space.$`^\mathrm{?}`$ The most important correlator is the quadratic one $`\delta f\delta f_{t=t^{}}`$ which can be expressed in terms of $`\overline{f}`$ and a two particle correlation function.$`^{\mathrm{?},\mathrm{?}}`$

Some few comments. ($`i`$) The set of dynamical eqs. (6) - (12) is equivalent to the original set of microscopic equations. It is exact in that no approximations have been performed so far, and can be applied to both in- or out-of-equilibrium situations. ($`ii`$) The original equations as well as their split into mean gauge fields and fluctuations have been shown to be consistent with the non-Abelian mean field gauge symmetry.$`^{\mathrm{?},\mathrm{?}}`$ ($`iii`$) In order to obtain effective equations for the mean fields only, one has to integrate-out the fluctuations. This amount to solving explicitly the dynamical equations of $`\delta f`$ and the gauge field fluctuations for given initial conditions and to express the correlator terms explicitly as functions of the mean quantities, making use of the basic equal-time correlators. In this light, our procedure can be seen as a derivation of collision integrals. ($`iv`$) The dynamical equations for the fluctuations are non-linear as they involve fluctuations up to cubic order. Solving them will require some approximations. For small gauge coupling an expansion in powers of $`g`$ can be performed. If the fluctuations remain small within a coarse-graining volume, the second moment approximation can be used, which amount ultimately to a linearisation of the dynamics of the fluctuations.$`^\mathrm{?}`$ These types of approximations are systematic and consistent with the non-Abelian gauge symmetry of the mean fields. ($`v`$) In the Abelian limit $`\xi `$ vanishes identically, while $`\eta `$ reduces to the well-known Balescu-Lennard collision integral.$`^\mathrm{?}`$ ($`vi`$) In the limit where fluctuations are fully neglected, the set of equations are known as the (non-Abelian) Vlasov equations.$`^{\mathrm{?},\mathrm{?}}`$

## 5 The hot non-Abelian plasma close to equilibrium

We now turn to the specific example of a hot thermal plasma close to equilibrium. ‘Hot’ implies that particle masses can be neglected to leading order $`m(T)T`$, and that the gauge coupling, as a function of temperature, is very small $`g(T)1`$. The non-Abelian plasma is characterised by two dimensionless parameters, the gauge coupling $`g`$ and the plasma parameter $`ϵ`$, which obtains as the ratio of the volume per particle over a Debye volume. The Debye (or screening) length describes the polarisation of the plasma. A small plasma parameter is mandatory for a kinetic description to be viable. Physically, this means that a large number of particles interact coherently within a Debye volume (quasi-particle behaviour). For a classical plasma, $`ϵ`$ is an independent parameter related to the mean particle number and the gauge coupling.$`^\mathrm{?}`$ For a quantum plasma, the mean particle number can not be fixed arbitrarily, and $`ϵ`$ scales as $`ϵg^31`$. This explains why a kinetic description, for quantum plasmas, is intimately linked to the weak coupling limit.

To simplify our analysis we shall perform some approximations, all of which can be understood as a systematic expansion in $`g`$. After ensemble averaging, the distribution function $`f`$ is effectively coarse-grained over a Debye volume. Fluctuations of $`\overline{f}`$ within a Debye volume are parametrically suppressed by powers of $`g`$. We expand $`f`$ about the equilibrium distribution function $`\overline{f}^{\mathrm{eq}}`$ to leading order in $`g`$ as

$$f=\overline{f}^{\mathrm{eq}}+g\overline{f}^{(1)}+\delta f.$$
(13)
Solving (6) using (13) for $`\delta f=0`$ reproduces the HTL effective theory.$`^\mathrm{?}`$ Two further approximations are employed beyond the HTL level, when $`\delta f0`$. We shall discard cubic correlator terms in (6), which are additionally suppressed by a power in $`g`$, in favour of quadratic ones. This corresponds to neglecting three-particle collisions in favour of binary ones. We also employ the second moment approximation – setting $`\eta =\eta `$ and $`\xi =\xi `$ in (11) – which means that the effects of collisions are neglected for the dynamics of the fluctuations. All these approximations can be systematically improved to higher order.

Let us consider colour excitations, described by the colour current density

$$𝒥_a^\mu (x,𝐯)=gv^\mu 𝑑Q𝑑p_0d|𝐩|𝐩^2\mathrm{\Theta }(p_0)\delta (p^2)Q_af(x,p,Q),$$
(14)
and $`v^\mu p^\mu /p_0=(1,𝐯),𝐯^2=1`$. The current $`J`$ is obtained after an angle average over the directions of $`𝐯`$ as $`J(x)=\frac{d\mathrm{\Omega }_𝐯}{4\pi }𝒥(x,𝐯)`$. The colour measure $`dQ`$ is normalised $`𝑑Q=1`$ and contains the Casimir constraints of the gauge group, such as $`𝑑QQ_aQ_b=C_2\delta _{ab}`$ with $`C_2`$ the quadratic group Casimir (for more details, see ref. $`^{\mathrm{?},\mathrm{?}}`$). The approximate dynamical equations for the fluctuations become

$`\left(v^\mu \overline{D}_\mu \delta 𝒥^\rho \right)_a=m_D^2v^\rho v^\mu \left(\overline{D}_\mu a_0\overline{D}_0a_\mu \right)_agf_{abc}v^\mu a_\mu ^b\overline{𝒥}^{c,\rho },`$ (15)
$`\left(\overline{D}^2a^\mu \overline{D}^\mu (\overline{D}a)\right)_a+2gf_{abc}\overline{F}_b^{\mu \nu }a_{c,\nu }=\delta J_a^\mu ,`$ (16)
where the Debye mass $`m_D`$ obtains from the equilibrium distribution function and the group representation of the particles. Solving for the fluctuations in the present approximations yields a closed expression for $`\delta 𝒥`$, and an iterative expansion in powers of the background fields for $`a`$.$`^{\mathrm{?},\mathrm{?}}`$ The seeked-for dynamical equation for the mean quasi-particle colour excitations is

$`v^\mu \overline{D}_\mu \overline{𝒥}^0+m_D^2v^\mu \overline{F}_{\mu 0}`$ $`=`$ $`C_{\mathrm{lin}}[\overline{𝒥}^0]+\zeta (x,𝐯),`$ (17)
where the linearised collision integral $`C_{\mathrm{lin}}`$ can now be evaluated explicitly, using the expressions for $`a`$ and $`\delta 𝒥`$. The Yang-Mills equation remains unchanged at this order. For the collision integral, one finally obtains to linear order in the mean current, and to leading logarithmic order (LLO) in $`(\mathrm{ln}1/g)^11`$

$`C_{\mathrm{lin}}[\overline{𝒥}_a^0](x,𝐯)`$ $`=`$ $`gf_{abc}a_\mu ^b(x)\delta 𝒥^{c,\mu }(x,𝐯)|_{\overline{A}=0,\mathrm{linear}\mathrm{in}\overline{J},\mathrm{LLO}}`$ (18)
$`=`$ $`\gamma {\displaystyle \frac{d\mathrm{\Omega }^{}}{4\pi }I(𝐯,𝐯^{})\overline{𝒥}_a^0(x,𝐯^{})}`$
with the kernel $`I(𝐯,𝐯^{})=\delta ^2(𝐯𝐯^{})\frac{4}{\pi }(𝐯.𝐯^{})^2/\sqrt{1(𝐯.𝐯^{})^2}`$. Notice that the collision integral is local in coordinate space, but non-local in the angle variables. The rate $`\gamma =\frac{g^2T}{4\pi }\mathrm{ln}1/g`$ is (twice) the hard gluon damping rate. In obtaining (18), we introduced an infra-red cut-off $`gm_D`$ for the elsewise unscreened magnetic sector. This result has been first obtained by Bödeker in ref. $`^\mathrm{?}`$. In addition, the source for stochastic noise $`\zeta `$ in (15) can be identified as

$`\zeta _a(x,𝐯)`$ $`=`$ $`gf_{abc}a_\mu ^b(x)\delta 𝒥^{c,\mu }(x,𝐯)|_{\overline{A}=0,\overline{J}=0}`$ (19)
$`\zeta ^a(x,𝐯)\zeta ^b(y,𝐯^{})`$ $`=`$ $`2\gamma m_D^2TI(𝐯,𝐯^{})\delta ^{ab}\delta (xy),`$ (20)
and the self-correlator (20) has been evaluated from (19) to LLO. Notice that its strength is determined by the kernel of the dissipative process. Next, we show the close relationship to the fluctuation-dissipation theorem (FDT).

## 6 Fluctuation-dissipation theorem

It is well-known that dissipation in a (quasi-stationary) plasma is intimately linked to the fluctuations, a link which is given by the FDT. For the coarse-grained kinetic equation (3) one may ask how the spectral functions of $`\zeta `$ and $`f`$ have to be related to $`C[f]`$ in order to satisfy the FDT. This question has been considered in ref. $`^\mathrm{?}`$. Within classical transport theory, the FDT is implemented in a straightforward way.$`^\mathrm{?}`$ The pivotal element is the kinetic entropy $`S[f]`$. We consider small deviations of $`f`$ from the equilibrium, $`f=f^{\mathrm{eq}}+\mathrm{\Delta }f`$. The entropy, stationary at equilibrium, is expanded to quadratic order in $`\mathrm{\Delta }f`$, $`S=S_{\mathrm{eq}}+\mathrm{\Delta }S`$. We define the thermodynamical force as $`F(z)=\delta (\mathrm{\Delta }S)/\delta (\mathrm{\Delta }f(z))`$, where $`z(x,𝐩,Q)`$. Given $`C[f](z)`$, it can be shown $`^{\mathrm{?},\mathrm{?}}`$ that

$$\zeta (z)\zeta (z^{})=\left(\frac{\delta C(z)}{\delta F(z^{})}+\frac{\delta C(z^{})}{\delta F(z)}\right)$$
(21)
is the noise self-correlator for $`\zeta `$ compatible with the FDT. For the particular example studied above, (20) follows from inserting (18) into (21). This proves that Bödeker’s effective theory is compatible with the FDT. Furthermore, the correlator $`\mathrm{\Delta }f\mathrm{\Delta }f|_{t=t^{}}`$ can be derived along the same lines and agrees, to leading order, with $`\delta f\delta f|_{t=t^{}}`$ from the Gibbs ensemble average. This guarantees that our formalism is consistent with FDT.

## 7 Discussion

We have presented a systematic approach to derive effective transport equations for non-Abelian plasmas. The procedure amounts to the ‘integrating-out’ of fluctuations about some mean fields. Most importantly, the formalism is applicable for both in- and out-of-equilibrium situations, simply because the main statistical information as encoded in the equal-time correlators of fluctuations does not depend on the mean distribution function being thermal, or not. At the same time, the formalism – and certain approximations to it – respects the underlying non-Abelian mean field symmetry, which is at the basis for any reliable computation.

For the close-to-equilibrium plasma, we have shown how Bödekers effective theory emerges to leading logarithmic accuracy. The collision integral as well as the necessary noise source have been derived explicitely from the microscopic theory. In addition, it is established on general grounds that the formalism (and hence Bödeker’s effective theory) is consistent with the fluctuation-dissipation theorem.

As a final comment we point out that the effective theory is the same for a classical or a semi-classical/quantum plasma, differing only in the equilibrium distribution function (Maxwell-Boltzmann vs. Bose-Einstein or Fermi-Dirac), and hence in the corresponding value for the Debye mass. The sole ‘quantum’ effect which entered our computation resides in the non-classical statistics of the particles, which is all that is needed to correctly describe a hot quantum non-Abelian plasma close to equilibrium at the present order of accuracy.

## References