Location via proxy:   [ UP ]  
[Report a bug]   [Manage cookies]                
Next Article in Journal / Special Issue
Levy-Lieb-Based Monte Carlo Study of the Dimensionality Behaviour of the Electronic Kinetic Functional
Previous Article in Journal / Special Issue
Energetic Study of Clusters and Reaction Barrier Heights from Efficient Semilocal Density Functionals
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Geometric Derivation of the Stress Tensor of the Homogeneous Electron Gas

1
Department of Physics, Temple University, Philadelphia 19122-1801, PA, USA
2
Department of Physics, University of Missouri, Columbia, MO 65211, USA
3
Theoretical Division & Center for Integrated Nanotechnologies, Los Alamos National Laboratory, Los Alamos, New Mexico, NM 87545, USA
*
Author to whom correspondence should be addressed.
Computation 2017, 5(2), 28; https://doi.org/10.3390/computation5020028
Submission received: 4 May 2017 / Revised: 23 May 2017 / Accepted: 2 June 2017 / Published: 8 June 2017

Abstract

:
The foundation of many approximations in time-dependent density functional theory (TDDFT) lies in the theory of the homogeneous electron gas. However, unlike the ground-state DFT, in which the exchange-correlation potential of the homogeneous electron gas is known exactly via the quantum Monte Carlo calculation, the time-dependent or frequency-dependent dynamical potential of the homogeneous electron gas has not been known exactly, due to the absence of a similar variational principle for excited states. In this work, we present a simple geometric derivation of the time-dependent dynamical exchange-correlation potential for the homogeneous system. With this derivation, the dynamical potential can be expressed in terms of the stress tensor, offering an alternative to calculate the bulk and shear moduli, two key input quantities in TDDFT.

1. Introduction

Time-dependent density functional theory (TDDFT) [1,2] is the most important extension of Kohn–Sham [3] ground-state DFT to excited or time-dependent states. In this theory, only the time-dependent dynamical exchange-correlation (xc) potential must be approximated. The simplest construction of the time-dependent dynamical potential is the adiabatic approximation [4], which takes the same form of the ground-state xc potential, but with the ground-state density n 0 ( r ) replaced by the instantaneous time-dependent density n ( r , t ) ,
v x c ad ( [ n ] ; r , t ) = δ E x c [ n 0 ] / δ n 0 ( r ) | n 0 ( r ) = n ( r , t ) .
Because of its simplicity in both theoretical construction and numerical implementation, the adiabatic TDDFT has been most widely used in the study of low-lying single-particle excitations [5]. Moreover, it has been shown [6,7] that (at least for small systems) the largest source of error in the prediction of low-lying excitation energies arises from the approximation for the ground-state xc potential, but not from the frequency dependence of the dynamical xc potential. However, this approximation fails to describe multi-particle [8,9] and high-lying excitations, because the retardation effect, which is important in these situations, has been ignored. Since it disregards the frequency dependence and thus the history of the system prior to time t, the adiabatic approximation is unable to describe energy dissipation [10,11,12]. To properly treat the retardation effect, one must go beyond the adiabatic approximation and construct a frequency-dependent non-adiabatic correction for the dynamical potential.
One way to develop the frequency-dependent part [13,14] of the dynamical xc potential is through the time-dependent xc stress tensor [15,16] P μ ν of quantum hydrodynamics. In this approach, the exchange-correlation vector potential A x c ( r , ω ) (whose derivative with respect to time yields the exchange-correlation electric field) is represented as the sum of two contributions [17,18]:
i ω e c A x c , μ ( r , ω ) = μ V x c ad ( r , ω ) + 1 n 0 ν Δ P μ ν x c ( r , ω ) ,
where e is the absolute value of the electron charge, c is the speed of light and μ , ν are Cartesian indices with the convention that repeated indices are summed over. The first term V x c ad is the adiabatic xc potential, the gradient of which gives the vector xc potential, which only has implicit frequency dependence via the time-dependent electron density. Δ P μ ν is the dynamical part of the xc stress tensor, the divergence of which defines the frequency-dependent part to the xc field. This is the core part of TDDFT that allows explicit frequency dependence. The dynamical xc stress tensor is defined by:
Δ P μ ν x c ( r , ω ) = P μ ν x c ( r , ω ) P μ ν x c , ad ( r , ω ) .
Here, P μ ν x c ( r , ω ) is the full dynamical xc stress tensor, while P μ ν x c , ad ( r , ω ) is the adiabatic part constructed from the equilibrium stress tensor, i.e.,
P μ ν x c , ad ( r , t ) = P μ ν x c , eq [ n 0 ] | n 0 = n ( r , t ) ,
with n ( r , t ) being the instantaneous electron density. Numerical calculations [19,20] have shown that the frequency-dependent part of the dynamical potential often leads to a good improvement over the adiabatic approximation.
A time-dependent dynamical stress tensor consists of two contributions: the infinite-frequency part and the finite-frequency part. Tokatly [15,16] presented a geometric formulation of TDDFT based on a nonlinear coordinate transformation, the Lagrangian coordinate transformation, and calculated the infinite-frequency part for the homogeneous system. The advantage of working in a co-moving frame is that the particle density is time independent, and the expectation value of current operator identically vanishes. As a result, the many-body dynamics looks more similar to a static problem. The strong nonlocality of the stress tensor in the laboratory frame is hidden in the space-time dependence of the deformation tensor in the Lagrangian frame (or co-moving frame). We have extended [11,21] this geometric approach to the inhomogeneous system and calculated the infinite-frequency part of the stress tensor Δ P μ ν x c ( r , ω ) . In this paper, we present a detailed geometric derivation of the whole spectra of the dynamical stress tensor for the homogeneous system, aiming to provide a simpler way to construct the time-dependent dynamical xc potential [22,23,24] for the homogeneous density and a better understanding of time-dependent processes beyond the adiabatic TDDFT. Although we are giving a hydrodynamic representation of the xc field, the treatment of time-dependent or excited processes is still within TDDFT, but not orbital-free quantum hydrodynamics.

2. Results and Discussion

2.1. Homogeneous Electron Gas in a Uniform Metric Field

Making use of the equation of motion for the current, the dynamical stress tensor may be calculated [24] from the current-current response function [2]. Here, we present a simpler approach for the dynamical stress tensor of the system based on the geometric formulation of the many-body dynamics.
The dynamical properties of the stress tensor of the electron gas are connected to the deformation of the ground-state density, which is measured by the covariant metric tensor or Green’s deformation tensor g μ ν (for brevity, the time argument is not displayed explicitly). The metric tensor is defined by [15,16]:
g μ ν r a ξ μ · r a ξ ν , g μ ν ξ μ r a · ξ ν r a ,
where the new coordinate ξ is defined as the initial position of a fluid element. Because the instantaneous position r of the fluid element is the vector sum of the initial position ξ and the elastic displacement u of the fluid element, we have:
ξ = r u ( r ) .
Note that g μ ν g ν λ = g μ λ δ μ λ and g = d e t ( g μ ν ) = [ d e t ( g μ ν ) ] 1 . The linearization of the metric tensor leads to:
g μ ν = δ μ ν + δ g μ ν , δ g μ ν = ν u μ + μ u ν ,
g μ ν = δ μ ν + δ g μ ν , δ g μ ν = δ g μ ν .
The combination ν u μ + μ u ν is twice the strain tensor of elasticity theory, as defined by [25]. However, we find it more convenient to define:
u μ ν δ g μ ν , u μ ν δ g μ ν .
From Equations (5)–(9), we see that the assumption of g μ ν being homogeneous is equivalent to the assumption of the strain field u μ ν being homogeneous. In other words, the displacement field is a linear function of r .
A useful property [16] of the metric determinant g is:
δ g g = δ ln g = Tr ( δ ln g ) = g μ ν δ g μ ν = g μ ν δ g μ ν .
The square root of the metric determinant g is the ratio between the volume element in Euclidean coordinates and the corresponding volume element in the curvilinear coordinates, as seen from Equation (5). These relationships can be obtained with some simple algebra.
The Hamiltonian of the electron gas in a homogeneous metric field g μ ν can be written as:
H ˜ ^ = k g μ ν k μ k ν 2 m a ^ ˜ k a ^ ˜ k + 1 2 V g k , q v ( q ) a ^ ˜ k n ^ ˜ q a ^ ˜ k q .
where V is the invariant d-dimensional volume of the system and n ^ ˜ q k a ^ ˜ k q a ^ ˜ k . The first term is the kinetic contribution in the co-moving frame, while the second is the interaction potential term, with v ( q ) being the Fourier transform of the electron-electron Coulomb interaction. The norm of vectors is given by:
q g μ ν q μ q ν ,
so the kinetic energy term is simply k 2 / 2 m. The factor g 1 / 2 in front of the potential term arises from the matrix element of the interaction between “plane waves” in the curvilinear coordinates (see Equation (19) below).
In the Hamiltonian of Equation (11), the operators a ^ ˜ k and a ^ ˜ k refer to plane wave states in coordinates ξ : ϕ k ( ξ ) = e i k μ ξ μ / V . This Hamiltonian has been constructed in the following manner. First, we start with the homogeneous electron gas in coordinates r on a 3 D torus of size L × L × L = V . The one-particle wave functions are plane waves ϕ k ( r ) = e i k r / V with k satisfying the usual quantization condition. The operator a ^ ˜ k creates a particle in this state. Second, we make a coordinate transformation (6) from r to ξ . This transformation is implicitly linear if the metric g μ ν is homogeneous. The transformed domain is still a torus. The normalized wave functions have the form:
ψ ˜ k ( ξ ) = g 1 / 4 V e i k · r ( ξ ) ,
where we have used the fact that:
d ξ = g 1 / 2 d r ,
as indicated by Equation (5). For a linear transformation, this can also be written as:
ψ ˜ k ( ξ ) = g 1 / 4 V e i k ˜ μ ξ μ ,
where:
k ˜ μ = ξ μ r ν k ν , d k ˜ = g 1 / 2 d k .
The Hamiltonian (11) is simply the standard electron-gas Hamiltonian in which the operators a ^ k have been replaced by the operators g 1 / 4 a ^ ˜ k ˜ , i.e.,
a ^ k g 1 / 4 a ^ ˜ k ˜ ,
Thus, wherever we have a plane wave of wave vector k , we can have the transformed wave function (15). This substitution gives:
H ˜ ^ = g 1 / 2 k ˜ k 2 2 m a ^ ˜ k ˜ a ^ ˜ k ˜ + g 2 V k ˜ , k ˜ , q ˜ v ( q ) a ^ ˜ k ˜ + q ˜ a ^ ˜ k ˜ q ˜ a ^ ˜ k ˜ a ^ ˜ k ˜ .
To show that this is indeed equivalent to Equation (11), notice that k 2 = g μ ν k ˜ μ k ˜ ν = g μ ν k ˜ μ k ˜ ν = k ˜ or, more compactly, | k | = k ˜ . Therefore, we have,
H ˜ ^ = g 1 / 2 k ˜ g μ ν k ˜ μ k ˜ ν 2 m a ^ ˜ k ˜ a ^ ˜ k ˜ + g 2 V k ˜ , k ˜ , q ˜ v ( q ˜ ) a ^ ˜ k ˜ + q ˜ a ^ ˜ k ˜ q ˜ a ^ ˜ k ˜ a ^ ˜ k ˜ .
Finally, recall that, in the thermodynamic limit,
k ˜ . . . = V ( 2 π ) d d k ˜ . . . = V ( 2 π ) d d k g 1 / 2 . . . = 1 g 1 / 2 k . . .
The integration variables k ˜ , k ˜ , etc., are immaterial and can be renamed k , k , etc. The factors g 1 / 2 combine to produce the pre-factors observed in Equation (11).
To write the Hamiltonian in “first quantization”, we must keep in mind that the operators a ^ ˜ k create electrons in plane wave states e i k μ ξ μ / V in the new coordinates with the original normalization. The operator that has matrix elements k μ k μ / 2 m in this basis is:
T ^ = 1 2 m i ξ i μ ξ i μ .
This is the first quantized expression for the kinetic energy. Similarly, for the potential energy, we can verify that the operator with matrix element v ( q ) / ( V g 1 / 2 ) is:
W ^ = 1 2 i j e 2 ξ i ξ j .
This is the first quantized expression for the potential energy.

3. Time-Dependent Transformation

The homogeneous coordinate transformation from r to ξ defined by Equation (6) is time independent. Now, we ask what happens if we let our homogeneous coordinate transformation be time-dependent. At first sight, we should simply replace g μ ν by g μ ν ( t ) and g by g ( t ) in Equation (11). However, this is not the complete solution. The problem is that the operators a ^ ˜ k are explicitly time-dependent, and this gives a contribution to the time evolution that is not taken into account by the Hamiltonian of Equation (11). Equation (11) must therefore be amended if the Hamiltonian is to generate the correct time evolution. Since a ^ ˜ k creates a particle in the state e i k · ξ , we can write a ^ ˜ k = e i k · ξ . Performing time derivative on both sides of the state, we can deduce that:
i a ^ ˜ k t = k μ ξ ˙ μ a ^ ˜ k
where ξ ˙ ξ / t is the explicit derivative of the coordinate transformation ξ = ξ ( r , t ) with respect to time.
This requires an additional single-particle term in the Hamiltonian. Taking into account the anti-commutation relation [ a ^ ˜ ^ k , a ^ ˜ k ] + = 1 , we see that the additional term must have the form:
ξ ˙ μ k k μ a ^ ˜ k a ^ ˜ k .
Therefore, the complete time-dependent Hamiltonian must include a vector-potential term [25]:
H ˜ ^ ( t ) = k ( k μ + m ξ ˙ μ ) ( k μ + m ξ ˙ μ ) 2 m m 2 ξ ˙ μ ξ ˙ μ a ^ ˜ k a ^ ˜ k + 1 2 V g k , q v ( q ) a ^ ˜ k n ˜ ^ q a ^ ˜ k q .
In the homogeneous electron gas, the vector potential term is inconsequential (i.e., it gives only an overall phase change), although it is significantly important [11] in the analysis of the inhomogeneous electron gas.
A specific example of coordinate transformation is the transformation to a Lagrangian frame or co-moving frame, as first introduced by Tokatly [15,16] and later used by Vignale and co-workers [11,21,26,27] to calculate the dynamical stress tensor of an inhomogeneous system in the infinite-frequency limit. Let j ( r ) be the current density (expectation value of the current-density operator) and v ( r ) = j ( r ) / n ( r ) the velocity field. Then, the transformation from r (laboratory frame) to ξ (co-moving frame) is defined by the solution of the differential equation:
r ( ξ , t ) t = v ( r ( ξ , t ) ) , r ( 0 ) = ξ .
This defines r as a function of ξ and t and, upon inversion (assuming the latter to be possible) ξ as a function of r and t. The explicit derivative ξ ˙ ( r , t ) tells us that we should modify the initial condition by d ξ = ξ ˙ ( r , t ) d t in order that the trajectory defined by Equation (23) arrives at r at time t + d t . Performing the variations in tand ξ at constant r gives:
r μ x i ν d ξ ν + r μ t d t = 0 ,
leading to:
ξ ˙ μ = ξ μ r ν v ν ( ξ , t ) v ˜ μ .
Therefore, the Hamiltonian of Equation (22) takes the form:
H ˜ ^ ( t ) = k ( k μ m v ˜ μ ) ( k μ m v ˜ μ ) 2 m m 2 v ˜ μ v ˜ μ a ^ ˜ k a ^ ˜ k + 1 2 V g k , q v ( q ) a ^ ˜ k n ˜ ^ q a ^ ˜ k q .

4. Macroscopic Stress Tensor

We now define the macroscopic stress tensor operator:
P ^ μ ν 2 V g δ H ˜ ^ δ g μ ν .
Making use of the identity (10) and:
δ q δ g μ ν = q μ q ν 2 q ,
we can easily obtain:
P ^ μ ν = 1 V g k k μ k ν m a ^ k a ^ k + 1 2 V g k , q [ q μ q ν q v ( q ) + g μ ν v ( q ) ] a ^ k n ^ q a ^ k q ,
where v ( q ) = d v ( q ) / d q . For Coulomb interaction in d spatial dimensions ( d = 2 , 3 ), v q = A d / q d 1 and q v q = ( d 1 ) v q . Thus, we have:
P ^ μ ν = 1 V g k k μ k ν m a ^ k a ^ k + 1 2 V g k , q v ( q ) × g μ ν ( d 1 ) q μ q ν q 2 a ^ k n ^ q a ^ k q ,
The stress tensor P ˜ μ ν ( t ) in the curvilinear reference frame, not to be confused with the stress tensor operator, is the expectation value of P ^ μ ν in the state | ψ ˜ ( t ) that evolves under the Hamiltonian H ˜ ^ :
P ˜ μ ν ( t ) ψ ˜ ( t ) | P ^ μ ν | ψ ˜ ( t ) .
Notice that, by construction, this object is independent of position. In the flat space of the laboratory frame, the stress tensor is obtained by transforming back to the laboratory frame, i.e.,
P μ ν ( t ) = ξ α r μ ξ β r ν P ˜ α β ( t ) .
In the linearized theory, we can write:
P ˜ μ ν ( t ) = P ˜ μ ν e q + P ˜ μ ν ( 1 ) ( t ) ,
where P ˜ μ ν e q is calculated in the equilibrium state of H ˜ ^ (with Euclidean metrics) and is independent of time:
P ˜ μ ν e q = 1 V k k μ k ν m n k + 1 2 V k , q q μ q ν q v ( q ) + δ μ ν v ( q ) × ρ 2 ( k , q ) ,
where n k a ^ k a ^ k e q is the momentum occupation number in the equilibrium state and:
ρ 2 ( k , q ) a ^ k n ^ q a ^ k q e q
is a two-particle density matrix related to the equilibrium exchange-correlation hole. This can be rewritten in terms of the static structure factor S ( q ) as follows:
P ˜ μ ν e q = 1 V k k μ k ν m n k + n 2 q q μ q ν q v ( q ) + δ μ ν v ( q ) × [ S ( q ) 1 ] ,
where n = N / V is the electron density of the homogeneous system, invariant in the curvilinear reference frame. In the special case of the Coulomb interaction ( v ( q ) 1 / q d 1 ), it reduces to:
P ˜ μ ν e q = 1 V k k μ k ν m n k + n 2 q v ( q ) δ μ ν ( d 1 ) q μ q ν q 2 × [ S ( q ) 1 ] = 2 T + W d V δ μ ν ,
where T = k ( k 2 / 2 m ) n k is the kinetic energy, and W = ( N / 2 ) q v ( q ) [ S ( q ) 1 ] is the potential energy of the homogeneous electron gas.
Then, we transform the stress tensor (Equation (36)) from the co-moving frame to the Laboratory or static frame. This leads to the expression for the stress tensor in the laboratory frame,
P μ ν ( t ) = P ˜ μ ν e q γ u μ P ˜ γ ν e q γ u ν P ˜ μ γ e q + P ˜ μ ν ( 1 ) ( t ) = 2 T + U d V δ μ ν + u μ ν + P ˜ μ ν ( 1 ) ( t ) .
Here, the first term is the adiabatic contribution to the stress tensor, while the second arises from the non-adiabatic or frequency-dependent contribution. This completes our formal derivation of the stress tensor. In the following, we will focus on the nonadiabatic part.

5. Tensor of Elasticity

A central concept in linear elasticity theory is the rank-four tensor of elasticity, Q ˜ μ ν α β , which connects the stress tensor to the strain tensor. In the frequency domain, there is a simple relationship between the stress tensor and the strain field,
P ˜ μ ν ( 1 ) ( ω ) = Q ˜ μ ν α β ( ω ) u α β ( ω ) ,
the real-time version of which being:
P ˜ μ ν ( 1 ) ( t ) = t Q ˜ μ ν α β ( t t ) u α β ( t ) d t .
To understand the structure of the tensor of elasticity, let us go back to Equation (29), which gives the explicit form of the stress tensor operator for arbitrary homogeneous metric fields. We see that a change in metric δ g α β = u α β is immediately reflected in a change in the form of P ^ μ ν . This is a frequency-independent contribution to the tensor of elasticity. It is the analogue of the “diamagnetic term”, which appears in the current response to a vector potential field. An additional contribution to Q comes from the change of the wave function in response to the metric field. This is a truly retarded (i.e., frequency-dependent) contribution, which can be expressed in terms of the stress-stress response function. It is expected to vanish at high frequency. The tensor of elasticity is thus the sum of the two parts,
Q ˜ μ ν α β ( ω ) = Q ˜ μ ν α β + Δ Q ˜ μ ν α β ( ω ) ,
where:
Q ˜ μ ν α β δ P ˜ ^ μ ν δ g α β g = 1 ,
and:
Δ Q ˜ μ ν α β ( ω ) V 2 P ^ μ ν ; P ^ α β ω .
( P ^ μ ν = P ˜ ^ μ ν | g = 1 .) Both the average and the response function are calculated in the equilibrium state of H ˜ ^ , with a Euclidean metric ( g = 1 ) . Notice that the tensor of elasticity in the frequency domain has the same physical dimensions as the stress tensor in real time.

5.1. High-Frequency Limit

From Equation (29) and Identities (10) and (28) and noting that:
δ g μ ν δ g α β = δ g μ ν δ g α β = g μ α g ν β + g μ β g ν α 2 ,
we calculate Q ˜ μ ν α β directly from its definition (41) and obtain:
Q ˜ μ ν α β = δ α β 2 V k k μ k ν m n k + n 2 q { v ( q ) [ δ μ ν δ α β δ μ α δ ν β + δ μ β δ ν α 2 ] + v ( q ) [ δ μ ν q α q β 2 q + δ α β q μ q ν q ] + q μ q ν q α q β 2 q 2 v ( q ) v ( q ) q } × [ S ( q ) 1 ] .
Specializing to the case of the Coulomb interaction and taking into account the additional terms of Equation (37), we see that the complete tensor of elasticity at high frequency in the laboratory frame is given by:
Q μ ν α β = 2 T + W 2 d V ( δ μ α δ ν β + δ μ β δ ν α ) + Q ˜ μ ν α β ,
where Q ˜ μ ν α β arises from the coordination transformation from the co-moving frame to the laboratory frame. Due to the rotational symmetry, the tensor of elasticity can actually be expressed in terms of just two constants, the bulk modulus K and the shear modulus μ , in the following manner:
Q μ ν α β = K 2 δ μ ν δ α β + μ 2 δ μ α δ ν β + δ μ β δ ν α 2 d δ μ ν δ α β ,
where:
K = 2 d 2 Q μ μ α α
μ = 2 ( d 1 ) ( d + 2 ) Q μ ν μ ν 1 d Q μ μ α α .
Specializing to the high-frequency regime and substituting Equation (44) into Equation (45), we obtain:
Q μ μ α α = ( d + 2 ) T V + d + 1 2 W V
Q μ ν μ ν = ( d + 2 ) T V + d 3 2 W V
Therefore, the high-frequency elastic constants are:
K = 2 ( d + 2 ) d 2 T V + d + 1 d 2 W V
μ = 2 d T V ( d 1 ) d ( d + 2 ) W V .
Here, the kinetic (T) and potential (W) energies of the electron gas can be extracted via the coupling constant integration from the correlation energy, which was calculated [28,29] from Monte Carlo methods and accurately parametrized [30,31].

5.2. Finite-Frequency Spectra

Now, we turn to the finite-frequency dependent part of the elasticity tensor. In the following, we will evaluate the response function (42). We recall that:
A ^ ; B ^ ω i 0 d t [ A ^ ; B ^ ] e i ω t .
For this purpose, we split the stress tensor as the sum of the kinetic and potential parts and write the response function as:
P ^ μ ν ; P ^ α β ω = T ^ μ ν + W ^ μ ν ; T ^ α β + W ^ α β ω
where:
T ^ α β = 1 V k k α k β m a ^ k a ^ k
and:
W ^ α β = k , q q α q β q v ( q ) + δ α β ( q ) a ^ k n ^ q a ^ k q
are the kinetic and potential parts of the macroscopic stress tensor operator (at g = 1 ). Generally, the finite-frequency part of the tensor of elasticity, Δ Q ˜ μ ν α β ( ω ) , like the viscoelastic constants Δ K ( ω ) and Δ μ ( ω ) , is a complex function of frequency. Because of the Kramers–Krönig dispersion relation [2], we only need to consider the imaginary part of the response function P μ ν ; P α β ω .
Let us first evaluate the imaginary part, m P μ ν ; P α β ω . Making use of the identity [2]:
A ^ ; B ^ ω = i A ^ ˙ ; B ^ ω + [ A ^ , B ^ ] ω
= i A ^ ; B ^ ˙ ω + [ A ^ , B ^ ] ω ,
where A ^ ˙ is the time derivative of A ^ , and noting that the expectation values of the commutator terms are all real, we first write the imaginary part of the response function as:
m P ^ μ ν ; P ^ α β ω = 1 ω 2 m T ^ ˙ μ ν ; T ^ ˙ α β ω + 1 ω e [ T ^ ˙ μ ν ; W ^ α β ω W ^ μ ν ; T ^ ˙ α β ω ] + m W ^ α β ; W ^ α β ω .
Since T ^ ˙ α β ; W ^ α β ω = W ^ α β ; T ^ ˙ α β ω , Equation (59) may be further simplified to:
m P ^ μ ν ; P ^ α β ω = 1 ω 2 m T ^ ˙ μ ν ; T ^ ˙ α β ω + 2 ω e T ^ ˙ μ ν ; W ^ α β ω + m W ^ α β ; W ^ α β ω .
From the Heisenberg equation of motion for the field operator:
i a ^ ˙ k = k 2 2 m a ^ k + 1 V q , k v ( q ) a ^ k a ^ k + q a ^ k q
we obtain:
T ^ ˙ α β = i V q , k v ( q ) k α k β m ( a ^ k + q n ^ q a ^ k a ^ k n ^ q a ^ k q ) ,
which may be re-expressed in terms of the current operator as:
T ^ ˙ α β = i V q v ( q ) ( j ^ q α q β + q α j ^ q β ) n ^ q .
Here, j ^ q α is the α -component of the current operator defined by:
j ^ q α k k α m a ^ k + q 2 a ^ k q 2 .
Because [ n ^ q , a k q ] = N ^ , a constant of motion, the interaction part of the stress tensor of Equation (56) may be simplified as:
W ^ α β = 1 2 V q v ( q ) δ α β ( d 1 ) q α q β q 2 n ^ q n ^ q .
Substituting Equations (63) and (65) into the response function of Equation (60), we obtain an expression for the tensor of elasticity in terms of four-point functions. These four-point response functions contain the interaction up to second order in the Coulomb interaction. We conclude that dissipation and, hence, retardation are absent up to first order in the Coulomb interaction. In particular, they are absent at the exact exchange level. We stress that this conclusion is only valid for the homogeneous electron gas. In a non-uniform electronic system, the time derivative of the kinetic energy is non-zero even in the absence of electron-electron interactions. As a result, one can have retardation and dissipation even at the exchange-only level.
Calculation of the four-point response functions is difficult. As a first step, we resort to the mode-coupling approximation. This approximation (at T = 0 ) has been previously discussed in the literature [22,23] and can be stated as follows:
m A ^ B ^ ; C ^ D ^ ω 1 π 0 ω d ω [ m A ^ ; C ^ ω m B ^ ; D ^ ω ω + m A ^ ; D ^ ω m B ^ ; C ^ ω ω ] .
Thus, we reduce the calculation of a very complicated four-point response function to the calculation of simpler two-point response functions.
Making use of the mode-decoupling approximation, together with various relations of response functions [2]:
j ^ q α ; j ^ q β ω = χ α β J ( q , ω ) δ q q ,
n ^ q ; n ^ q ω = ( q 2 / ω 2 ) χ L ( q , ω ) δ q q ,
n ^ q ; j ^ q α ω = ( q α / ω ) χ L ( q , ω ) δ q q ,
j ^ q α ; n ^ q ω = ( q α / ω ) χ L ( q , ω ) δ q q ,
where χ α β J = [ χ L ( q , ω ) q α q β / q 2 + χ T ( q , ω ) ( δ α β q α q β / q 2 ) ] is the current-current response, we express T ^ ˙ α β ; T ^ ˙ α β in terms of the longitudinal (L) and transverse (T) current-current response functions:
m T ^ ˙ μ ν ; T ^ ˙ α β ω 1 π V 2 q v ( q ) 2 q 4 0 ω d ω { 4 ω ω ( ω ω ) 2 × m χ μ ν L ( q , ω ) m χ α β L ( q , ω ω ) + 2 ( d 1 ) ( ω ω ) 2 m χ μ ν T ( q , ω ) m χ α β L ( q , ω ω ) } ,
with the help of the convolution rule [32],
f ( ω ω ) × g ( ω ) = 0 d ω f ( ω ω ) g ( ω ) = 0 d ω g ( ω ω ) f ( ω ) = g ( ω ω ) × f ( ω ) .
Similarly, we obtain:
e T ^ ˙ μ ν ; W ^ α β ω 1 π V 2 q v ( q ) 2 q 4 0 ω d ω 2 ( 2 d ) ω 2 ( ω ω ) × m χ μ ν L ( q , ω ) m χ α β L ( q , ω ω ) .
and:
m W ^ μ ν ; W ^ α β ω 1 2 π V 2 q v ( q ) 2 q 4 0 ω d ω d 2 3 d + 3 ω 2 ( ω ω ) × m χ μ ν L ( q , ω ) m χ α β L ( q , ω ω ) .
Next, we evaluate the imaginary part of the component, m P ^ α α ; P ^ α α ω . According to Equations (54)–(56), we write:
m P ^ α α ; P ^ α α ω = m 2 T ^ + W ^ ; 2 T ^ + W ^ ω .
where T ^ = ( 1 / 2 m ) k k 2 a ^ k a ^ k is the kinetic energy operator and W ^ = ( 1 / 2 ) k , q v q a ^ k n ^ q a ^ k q the potential operator. For the uniform electron gas, the equilibrium Hamiltonian can be written as H ^ = T ^ + W ^ + V b , where V b is the external potential due to the background. Because V b is a constant of motion, Equation (75) can be simplified to:
m P ^ α α ; P ^ α α ω = m W ^ ; W ^ ω ,
where:
m W ^ ; W ^ ω = 1 4 q , q v q v q m n ^ q n ^ q ; n ^ q n ^ q ω .
Employing the mode-decoupling approximation of Equation (66) and the relation of Equation (68), we obtain:
m W ^ μ ν ; W ^ α β ω 1 2 q v q 2 q 4 0 ω d ω ω 2 ( ω ω ) × m χ L ( q , ω ) m χ L ( q , ω ω ) .
Finally, we can express the imaginary parts of the tensor of elasticity and the viscosity spectra in terms of the longitudinal and transverse current-current response functions. The result is:
m K ˜ ( ω ) = 1 2 d 2 q v ( q ) 2 0 ω d ω π q 2 ω 2 m χ L ( q , ω ) × q 2 ( ω ω ) 2 m χ L ( q , ω ω ) ,
for d = 2 (i.e., 2 D ) or 3 (i.e., 3 D ),
m μ ˜ ( ω ) = q v ( q ) 2 0 ω d ω π [ 8 15 q 2 ω 2 m χ L ( q , ω ) + 2 5 q 2 ω 2 m χ T ( q , ω ) ] q 2 ( ω ω ) 2 × m χ L ( q , ω ω ) ,
for 3D systems, and:
m μ ˜ ( ω ) = q v ( q ) 2 0 ω d ω π [ 9 16 q 2 ω 2 m χ L ( q , ω ) + 1 2 q 2 ω 2 m χ T ( q , ω ) ] q 2 ( ω ω ) 2 × m χ L ( q , ω ω ) ,
for 2D systems.
The xc kernel f x c L ( T ) can be immediately obtained from the relations:
m K ˜ ( ω ) = n 0 2 m f x c L ( ω ) 2 ( d 1 ) d m f x c T ( ω ) .
m μ ˜ ( ω ) = n 0 2 m f x c T ( ω ) ,
with the help of Equation (19) and Equations (79)–(81). To reproduce the correct high-frequency ( ω ) limit, Nifosì, Conti and Tosi (NCT) [22,23] introduced a frequency-dependent prefactor:
g x ( ω ) = 1 + ( 0.5 ω / 2 ϵ F ) 1 + ( ω / 2 ϵ F ) ,
which has the property of g x 1 / 2 as ω and g x 1 as ω 0 ( ϵ F = k F 2 2 / 2 m is the Fermi energy). The final expression for the xc kernel is given by:
m f xcL ( T ) ( ω ) = g x ( ω ) 0 ω d ω π d d q ( 2 π ) d n 0 2 v 2 ( q ) [ a L ( T ) q 2 ω 2 × m χ L ( q , ω ) + b L ( T ) q 2 ω 2 m χ T ( q , ω ) ] × q 2 ( ω ω ) 2 m χ L ( q , ω ω ) ,
where a L = 23 / 30 , a T = 8 / 15 , b L = 8 / 15 , b T = 2 / 5 for 3 D and a L = 11 / 16 , a T = 9 / 16 , b L = b T = 1 / 2 for 2 D . Apart from the prefactor g x ( ω ) , Equation (85) was derived by several authors [23,33,34] based on different approaches.
The infinite frequency limit of the elastic constants, combined with the viscosity spectra, determines the frequency-dependent elastic constants by virtue of the dispersion relations and, thus, the frequency-dependent dynamical xc potential. In the following, we will study these two kernels in the low-frequency limit.

5.3. Low-Frequency Limit of m f x c L ( T )

As discussed in Section 5.2, in the high-frequency limit, the longitudinal part of the xc kernel in the mode-coupling approximation is too small by a factor of two. To correct this error, NCT phenomenologically introduced a prefactor g x ( ω ) in the expression (85). This prefactor tends to one as ω 0 , leaving the low-frequency limit of the kernel unchanged. However, Qian and Vignale [35] (QV) have shown that in the low-frequency limit, the xc kernel m f x c L ( ω ) actually has two contributions: the “direct” term (D) and the “exchange” term (EX). Precisely, we have:
m f x c L ( ω ) = m f x c L D ( ω ) + m f x c L Ex ( ω ) .
The direct term is just the mode-coupling approximation (85) given by:
m f x c L D ( ω ) = ω e 4 m 2 k F 45 π 3 n 0 3 3 + 3 + λ 2 λ tan 1 λ ,
while the exchange term, which is missing in the mode-coupling approximation, is given by:
m f x c L Ex ( ω ) = ω e 4 m 2 k F 45 π 3 n 0 3 { 1 1 λ sin 1 ( λ ( 1 + λ ) 1 λ tan 1 λ + 1 λ 2 + λ π 2 tan 1 1 λ 2 + λ } .
Here, λ = 2 k F / k s , and k s = 4 k F / π a 0 is the Thomas–Fermi screening wave vector with a 0 = 2 / m e 2 being the Bohr radius. Substituting Equations (87) and (88) into Equation (86) yields the xc kernel in the low-frequency limit:
m f x c L ( ω ) = ω e 4 m 2 k F 45 π 3 n 0 3 { 5 λ + 5 λ tan 1 λ 2 λ sin 1 λ 1 + λ 2 + 2 λ 2 + λ 2 [ π 2 tan 1 1 λ 2 + λ 2 ] } .
Furthermore, QV have proven that as ω 0 ,
m f x c T ( ω ) / ω = d 2 ( d 1 ) m f x c L ( ω ) / ω .
To recover both the low- and high-frequency limits, based on the relations (86)–(90), we propose a new dimensionless prefactor:
g x ( ω ) = 1 + b 2 ( ω / 2 ϵ F ) 2 R + b ( ω / 2 ϵ F ) 2 ,
where:
R = lim ω m f x c L ( T ) D ( ω ) m f x c L ( T ) ( ω ) = 3 λ + 3 λ tan 1 λ / { 5 λ + 5 λ tan 1 λ 2 λ sin 1 λ 1 + λ 2 + 2 λ 2 + λ 2 × π 2 tan 1 1 λ 2 + λ 2 } .
The parameter b is determined by a fit to the zero-frequency limit of the bulk modulus, which can be calculated accurately from the xc energy per electron (see Equation (105)).

6. Exchange-Correlation Viscoelasticity Constants

The xc part of the stress tensor is defined as the difference between interacting and non-interacting stress tensors, minus the Hartree part. According to [15,16,17], the pure dynamical xc stress tensor Δ P μ ν (see Equation (3)) can be expressed in terms of the viscoelastic constants K x c ( ω ) and μ x c ( ω ) ,
Δ P μ ν x c = μ x c ( n 0 , ω ) ( u μ ν δ μ ν u α α / 3 ) + δ μ ν 2 [ K x c ( n 0 , ω ) K x c 0 ( n 0 ) ] u α α ,
where u α α = 2 · u . Similarly, the xc contributions to the viscoelastic constants, K x c ( ω ) and μ x c ( ω ) , are defined as the difference between the viscoelastic constants of interacting and non-interacting systems of the same electron density n 0 ( r ) , minus the Hartree potential part,
K x c ( ω ) = K ( ω ) K 0 , K 0 = 2 ( d + 2 ) d 2 T 0 ,
μ x c ( ω ) = μ ( ω ) μ 0 , μ 0 = 2 d T 0 ,
where K 0 and μ 0 are the viscoelastic constants of the non-interacting system, and:
T 0 = 3 10 ( 3 π 2 ) 2 / 3 n 0 5 / 3
is the non-interacting kinetic energy density for the 3D system.
Since the finite-frequency part arises from the second-order (Coulomb) interaction, it vanishes for a non-interacting system. Thus, the spectra of the xc viscoelastic constants can be expressed as:
K x c ( ω ) = Δ K x c ( ω ) + K x c ,
μ x c ( ω ) = Δ μ x c ( ω ) + μ x c ,
where:
Δ K x c ( ω ) = P d ω π m K ˜ ( ω ) ω ω ,
Δ μ x c ( ω ) = P d ω π m μ ˜ ( ω ) ω ω .
Here, “P” represents the Cauchy principal value. K x c and μ x c are the infinite-frequency parts, which can be calculated from Equations (51) and (52) as:
K x c = 2 ( d + 2 ) d 2 T c V + d + 1 d 2 W x c V ,
μ x c = 2 d T c V d 1 d ( d + 2 ) W x c V ,
Here, T c = T T 0 is the kinetic part of correlation with T being the kinetic energy density of the interacting system. W x c = W U H is the potential part of correlation, and U H is the Hartree potential energy. T c and W x c are related to the coupling-constant average of the xc energy density n ϵ x c via the coupling-constant integration formula [36], i.e.,
T c ( [ n 0 ] ; r ) = n 0 r s r s ϵ c ( r s ) ,
W x c ( r ) = n 0 ϵ x c + n 0 r s r s ϵ c ( r s ) ,
where r s ( r ) = a 0 1 ( 3 / 4 π n 0 ) 1 / 3 is the Seitz radius and a 0 is the Bohr radius (notice that the coupling constant λ controls the strength of the interaction λ e 2 / | r r | ).
Finally, the low-frequency limit of the bulk modulus K x c 0 ( n 0 ) can be calculated from:
K x c 0 ( n 0 ) = n 0 2 2 ( n 0 ϵ x c ) / n 0 2 ,
where the correlation energy per electron ϵ c ( r s ) can be obtained from the quantum Monte Carlo calculation, as pointed out above. In the next section, we will study the xc kernels within the random phase approximation (RPA).

7. Exchange-Correlation Kernel f x c L ( T ) within RPA

In this section, we calculate the xc kernel of Equation (85) within the RPA (random phase approximation). The frequency-dependent bulk and shear moduli obtained in this manner are exact in the high density limit ( r s 0 ) and provide a reasonable estimate at finite density. The shear modulus in the ω 0 limit is a very important local ingredient [11,21,37,38] of the dynamical xc potential. We first evaluate m K ˜ ( ω ) and m μ ˜ ( ω ) from Equations (82) and (83). Finally, the frequency-dependent bulk and shear moduli are calculated by substituting m K ˜ ( ω ) and m μ ˜ ( ω ) into Equations (97) and (98). We only focus on 3 D systems. 2 D systems can be treated in a similar manner.

7.1. Exchange-Correlation Kernel

The xc kernel can be calculated from RPA by replacing the response functions χ L ( T ) ( q , ω ) in Equation (85) with the RPA response functions χ L ( T ) RPA ( q , ω ) . In RPA, the response functions of the interacting system can be expressed in terms of the non-interacting response functions χ L ( T ) 0 ( q , ω ) as:
1 χ L ( T ) RPA ( q , ω ) + n / m = 1 χ L ( T ) 0 ( q , ω ) + n / m q 2 ω 2 v q L ( T ) ,
where for the Coulomb system, v q L = 4 π e 2 / q 2 for 3 D , v q L = 2 π e 2 / q for 2 D , and v q T = 0 . This immediately leads to:
m χ T RPA ( q , ω ) = m χ T 0 ( q , ω ) .
The non-interacting transverse m χ T 0 ( q , ω ) is given by [23]:
m χ T 0 ( q , ω ) = N ( 0 ) π ϵ F 8 m q ˜ [ ( 1 ν + 2 ) 2 θ ( 1 | ν + | ) ( 1 ν 2 ) 2 θ ( 1 | ν | ) ,
where N ( 0 ) = k s / ( 4 π e 2 ) = m k f / π 2 2 , ν ± = ω ¯ / q ¯ ± q ¯ / 2 , with q ¯ q / k F and ω ¯ ω / ( 2 ϵ F ) .
The current-current response χ L ( q , ω ) is related to the density-density response χ ( q , ω ) by:
χ L ( q , ω ) + n m = ω 2 q 2 χ ( q , ω ) ,
which is valid for both interacting and non-interacting systems. From Equations (106) and (109), we obtain:
m χ L RPA ( q , ω ) = ( ω 2 / q 2 ) m χ RPA ( q , ω ) = m χ 0 ( q , ω ) / | ϵ RPA ( q , ω ) | 2 .
Here, χ 0 ( q , ω ) is the Lindhard function, which can be written as:
χ 0 ( q , ω ) = e χ 0 ( q , ω ) + i m χ 0 ( q , ω ) .
The real part of χ 0 ( q , ω ) is given by:
e χ 0 ( q , ω ) = N ( 0 ) ( 1 2 1 ν 2 4 q ¯ lin | ν + 1 ν 1 | + 1 ν + 2 4 q ¯ lin | ν + + 1 ν + 1 | ) ,
while the imaginary part is given by:
m χ 0 ( q , ω ) = N ( 0 ) π 4 q ¯ [ ( 1 ν 2 ) θ ( 1 ν 2 ) ( 1 ν + 2 ) θ ( 1 ν + 2 ) ] .
Here, θ ( x ) is a step function having the property that θ ( x ) = 1 for x > 1 and 0 for x 0 . The dielectric function within the RPA is given by:
ϵ RPA ( q , ω ) = 1 v q L e χ 0 ( q , ω ) + i v q L m χ 0 ( q , ω ) .
Now, we turn to the evaluation of the xc kernel m f x c L ( T ) ( ω ) . Substituting Equation (107) and the first equality of Equation (110) into Equation (85) leads to:
m f x c L ( T ) ( ω ) = g x ( ω ) 4 0 ω 0 ω d q q 2 v ( q ) 2 ( 2 π ) 3 n 0 2 × [ a L m χ RPA ( q , ω ) + b L ( T ) q 2 ω 2 × m χ T 0 ( q , ω ) ] m χ RPA ( q , ω ) .
We can see from Figure 5.6 of [2] that χ L RPA ( q , ω ) terms make noticeable contributions to m f x c L ( T ) ( ω ) of Equation (115) only in a narrow range of ω for a given q. When we increase ω , m χ L 0 ( q , ω ) vanishes identically. However, as ω is close to some value of frequency, say Ω p , the real part of the dielectric function vanishes, i.e.,
ϵ 1 RPA ( q , Ω p ) = 1 v q L e χ 0 ( q , Ω p ) = 0 .
From the second equality of Equation (110), we can see that when e ϵ RPA ( q , Ω p ) vanishes, m χ L RPA ( q , ω ) becomes a δ -function δ ( ω Ω p ) . Because the δ -function δ ( ω Ω p ) is nonzero sharply at Ω p , it cannot be properly treated numerically. Therefore, we need to treat it separately. According to [2], we can write the RPA density-density response as the sum of two contributions:
m χ L RPA ( q , ω ) = m χ L 0 ( q , ω ) | ϵ RPA ( q , ω ) | 2 A ( q ) δ ( ω Ω p ) ,
where A ( q ) is defined by:
A ( q ) π v ( q ) | Ω p ϵ 1 RPA ( q , Ω p ) | ,
with Ω p = / Ω p and ϵ 1 RPA ( q , Ω p ) being the real part of ϵ RPA ( q , Ω p ) . Because this δ -function occurs at small q ¯ in the high-density region where RPA is valid (e.g, q ¯ 0 . 5 for r s = 1 ), expanding [2] the Lindhard function around q = 0 , up to order of q 2 , and then substituting it into the condition for the existence of δ -function (Equation (116)) yield:
Ω p 2 ω p 2 + 3 5 q 2 k F 2 ( / m ) 2 ,
where ω p = 4 π n 0 e 2 / m is the q = 0 plasmon frequency. Clearly, Ω p ω p . The exact Ω p -q relation can be obtained by solving Equation (116) for Ω p at each q. Inserting Equation (117) into Equation (115) leads to the final expression:
m f x c L ( T ) ( ω ) = g x ( ω ) 4 0 ω d ω 0 d q q 2 v ( q ) 2 ( 2 π ) 3 n 0 2 × { a L ( T ) m χ 0 ( q , ω ) | ϵ RPA ( ( q , ω ) | 2 a L ( T ) A ( q ) × δ ( ω Ω p ) + b L ( T ) q 2 ω 2 m χ T 0 ( q , ω ) } × [ m χ 0 ( q , ω ω ) | ϵ RPA ( ( q , ω ω ) | 2 A ( q ) × δ ( ω ω Ω p ) ]
We can separate the continuum electron-hole part from the plasmon contributions by neglecting the plasmon density-density response in Equation (120). The result is:
m f x c L ( T ) e h ( ω ) = g x ( ω ) 4 0 ω d ω 0 d q q 2 v ( q ) 2 ( 2 π ) 3 n 0 2 × { a L ( T ) m χ 0 ( q , ω ) | ϵ RPA ( ( q , ω ) | 2 + b L ( T ) q 2 ω 2 × m χ T 0 ( q , ω ) } m χ 0 ( q , ω ω ) | ϵ RPA ( ( q , ω ω ) | 2
In the following, we will calculate all other terms involving the plasmon contribution. Since the plasmon contribution vanishes when the plasmon enters the electron-hole continuum, the cutoff radius q c is determined by the condition [2] A ( q c ) = 0 . By inserting Equation (114) into Equation (118) and taking the q q c limit, an explicit k F q c relation is obtained as:
k s 2 k F 2 q ¯ c 2 1 2 + ( 2 + q c ) ln 2 + q ¯ c q ¯ c = 1 ,
where q ¯ c = q c / k F . For a given wavevector, we can find a corresponding cutoff radius q c by solving this equation. Define the largest q that satisfies ϵ 1 RPA ( Ω p , q ) = 0 (Equation (116)) as q p and the largest q that satisfies ϵ 2 RPA ( Ω p , q ) = 0 as q m . Then, we have q m < q p < q c .
The simplest term is the plasmon-plasmon contribution:
m f x c L ( T ) p p ( ω ) = g x ( ω ) 4 0 ω d ω 0 d q q 2 v ( q ) 2 ( 2 π ) 3 n 0 2 a L ( T ) × A 2 ( q ) δ ( ω Ω p ) δ ( ω ω Ω p ) .
Performing integration over ω yields:
m f x c L ( T ) p p ( ω ) = g x ( ω ) 4 0 ω d ω 0 d q q 2 v ( q ) 2 ( 2 π ) 3 n 0 2 a L ( T ) × A 2 ( q ) θ ( ω Ω p ) δ ( ω 2 Ω p ) .
which, through the change of variable q Ω p , can be written as:
m f x c L ( T ) p p ( ω ) = g x ( ω ) 4 ω Ω p ( q m ) d Ω Ω p q 2 v ( q ) 2 ( 2 π ) 3 n 0 2 a L ( T ) × A 2 ( q ) θ ( ω Ω p ) δ ( ω / 2 Ω p ) / 2 .
where Ω p ( q ) = d Ω p ( q ) / d q and q = q ( Ω p ) . Note that since q p < q c and Ω p ( q ) = 0 for q p < q q c , the upper limit of the q-integral should be Ω p ( q p ) , instead of Ω p ( q c ) . Furthermore, for q > q m , ϵ 2 RPA ( Ω p , q ) increases rapidly from zero. In this region, the δ -function is broadened, because ϵ 2 RPA ( Ω p , q ) > 0 , while ϵ 1 RPA ( Ω p , q ) remains zero in this region. Unlike the sharp δ -function, the broadened δ -function can be captured by numerical integration. Therefore, the contribution in the range of q m < q < q p has been already included in the electron-hole part (Equation (121)). Thus, the upper limit of the q-integral actually is q m .
m f x c L ( T ) p p ( ω ) = 2 g x 1 Ω p q 2 v ( q ) 2 ( 2 π ) 3 n 0 2 a L ( T ) A 2 ( q ) × θ ( ω 2 ω p ) θ ( 2 Ω p ω ) | Ω p = ω / 2 .
Note that the θ -function in Equation (125) has been dropped, because it is replaced by a more strict constraint imposed by the two θ -functions in Equation (126).
The electron-plasmon contribution is:
m f x c L ( T ) e p ( ω ) = g x ( ω ) 4 0 ω d ω 0 q c d q q 2 v ( q ) 2 ( 2 π ) 3 n 0 2 A ( q ) × { a L ( T ) m χ 0 ( q , ω ) | ϵ RPA ( ( q , ω ) | 2 + b L ( T ) q 2 ω 2 × m χ T 0 ( q , ω ) } δ ( ω ω Ω p )
leading to:
m f x c L ( T ) e p ( ω ) = g x ( ω ) 4 0 q c d q q 2 v ( q ) 2 ( 2 π ) 3 n 0 2 A ( q ) × { a L ( T ) m χ 0 ( q , ω Ω p ) | ϵ RPA ( ( q , ω Ω p ) | 2 + b L ( T ) q 2 ω 2 × m χ T 0 ( q , ω Ω p ) } δ ( ω Ω p )
The plasmon-electron contribution is:
m f x c L ( T ) p e ( ω ) = g x ( ω ) 4 0 ω d ω 0 q c d q q 2 v ( q ) 2 ( 2 π ) 3 n 0 2 A ( q ) × a L ( T ) δ ( ω Ω p ) m χ 0 ( q , ω ω ) | ϵ RPA ( ( q , ω ω ) | 2 ,
which, upon the integration over ω , leads to:
m f x c L ( T ) p e ( ω ) = g x ( ω ) 4 0 q c d q q 2 v ( q ) 2 ( 2 π ) 3 n 0 2 A ( q ) a L ( T ) × m χ 0 ( q , ω Ω p ) | ϵ RPA ( ( q , ω Ω p ) | 2 θ ( ω Ω p ) .
The imaginary part of the xc kernel is the sum of four contributions:
m f x c L ( T ) ( ω ) = m f x c L ( T ) e h ( ω ) + m f x c L ( T ) p p ( ω ) + m f x c L ( T ) e p ( ω ) + m f x c L ( T ) p e ( ω ) .
This completes our study of the xc kernel in the low-density limit within the RPA. Our formulation within the RPA should be accurate in the high-density regime. In the low-density regime, the RPA is not valid any more [39]. The correlation in the low-density regime presents a big challenge in TDDFT and even the ground-state DFT [40,41,42].

8. Conclusions

In conclusion, we have presented a derivation of the stress tensor of the uniform electron gas based on the geometric formulation of TDDFT. The derivation is based on the geometric reformulation of TDDFT proposed by Tokatly [11,15,16]. In this new formulation, we treated a dynamical problem as a static one by performing Lagrangian nonlinear coordination transformation of the system from the time-dependent laboratory frame to the co-moving frame. In this way, we can calculate the dynamical potential with the knowledge of the ground-state DFT. This offers a simpler way to calculate the dynamical xc potential.
Our derivation leads to an alternative form of the time-dependent dynamical xc potential for the homogeneous system and thus sheds light on a better treatment of the difficult problem for inhomogeneous systems. To better understand the properties, we have also studied the xc kernel in the high- and low-frequency limit and proposed a new interpolation formula to correct the exchange part of the kernel in the high-frequency limit (91). Our new kernel should be more accurate than that previously proposed by Nifosì, Conti and Tosi [22,23], because our new formula has been designed not only to recover the correct high-frequency limit, but also to recover the RPA low-frequency limit derived by Qian and Vignale [35], as well. Although in this work we limit ourselves to the linear response regime, this geometric approach can be extended to the nonlinear regime [43]. Because this work involves the ground-state DFT via infinite-frequency or high-frequency limits, it should be of broad interest to developers and users of TDDFT and DFT.

Acknowledgments

J.T. acknowledges support from the NSF under Grant No. CHE 1640584 and the Temple start-up via John Perdew. G.V. was supported by the DOE under Grant No. DE-FG02-05ER46203. J.X.Z. was supported by the Center for Integrated Nanotechnologies, a U.S. DOE Basic Energy Sciences user facility.

Author Contributions

J.T. drafted the manuscript. G.V. and J.Z. participated in writing the manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Runge, E.; Gross, E.K.U. Density-Functional Theory for Time-Dependent Systems. Phys. Rev. Lett. 1984, 52, 997–1000. [Google Scholar] [CrossRef]
  2. Giuliani, G.F.; Vignale, G. Quantum Theory of the Electron Liquid; Cambridge University Press: Cambridge, UK, 2005. [Google Scholar]
  3. Kohn, W.; Sham, L.J. Self-Consistent Equations Including Exchange and Correlation Effects. Phys. Rev. 1965, 140, A1133–A1138. [Google Scholar] [CrossRef]
  4. Zangwill, A.; Soven, P. Resonant Photoemission in Barium and Cerium. Phys. Rev. Lett. 1980, 45, 204–207. [Google Scholar] [CrossRef]
  5. Bauernschmit, R.; Ahlrichs, R. Treatment of electronic excitations within the adiabatic approximation of time dependent density functional theory. Chem. Phys. Lett. 1996, 256, 454–464. [Google Scholar] [CrossRef]
  6. Petersilka, M.; Gossmann, U.J.; Gross, E.K.U. Excitation Energies from Time-Dependent Density-Functional Theory. Phys. Rev. Lett. 1996, 76, 1212–1215. [Google Scholar] [CrossRef] [PubMed]
  7. Van Gisbergen, S.J.A.; Kootstra, F.; Schipper, P.R.T.; Gritsenko, O.V.; Snijders, G.; Baerends, E.J. Density-functional-theory response-property calculations with accurate exchange-correlation potentials. Phys. Rev. A 1998, 57, 2556–2571. [Google Scholar] [CrossRef]
  8. Tozer, D.J.; Handy, N.C. On the determination of excitation energies using density functional theory. Phys. Chem. Chem. Phys. 2000, 2, 2117–2121. [Google Scholar] [CrossRef]
  9. Tozer, D.J.; Amos, R.D.; Handy, N.C.; Roos, B.O.; Serrano-Andres, L. Does density functional theory contribute to the understanding of excited states of unsaturated organic compounds? Mol. Phys. 1999, 97, 859–868. [Google Scholar] [CrossRef]
  10. D’Agosta, R.; Vignale, G. Relaxation in time-dependent current-density-functional theory. Phys. Rev. Lett. 2006, 96. [Google Scholar] [CrossRef] [PubMed]
  11. Tao, J.; Vignale, G.; Tokatly, I.V. Time-dependent density functional theory: Derivation of gradient-corrected dynamical exchange-correlational potentials. Phys. Rev. B 2007, 76. [Google Scholar] [CrossRef]
  12. Ullirich, C.A. Time-dependent density-functional theory beyond the adiabatic approximation: Insights from a two-electron model system. J. Chem. Phys. 2006, 125. [Google Scholar] [CrossRef] [PubMed]
  13. Vignale, G.; Kohn, W. Current-dependent exchange-correlation potential for dynamical linear response theory. Phys. Rev. Lett. 1996, 77, 2037–2040. [Google Scholar] [CrossRef] [PubMed]
  14. Vignale, G.; Kohn, W. Electronic Density Functional Theory; Dobson, J.F., Das, M.P., Vignale, G., Eds.; Plenum Press: New York, NY, USA, 1996. [Google Scholar]
  15. Tokatly, I.V. Quantum many-body dynamics in a Lagrangian frame: I. Equations of motion and conservation laws. Phys. Rev. B 2005, 71. [Google Scholar] [CrossRef]
  16. Tokatly, I.V. Quantum many-body dynamics in a Lagrangian frame: II. Geometric formulation of time-dependent density functional theory. Phys. Rev. B 2005, 71. [Google Scholar] [CrossRef]
  17. Vignale, G.; Ullrich, C.A.; Conti, S. Time-dependent density functional theory beyond the adiabatic local density approximation. Phys. Rev. Lett. 1997, 79, 4878–4881. [Google Scholar] [CrossRef]
  18. Van Faassen, M.; de Boeij, P.L.; van Leeuwen, R.; Berger, J.A.; Snijders, J.G. Ultranonlocality in time-dependent current-density-functional theory: Application to conjugated polymers. Phys. Rev. Lett. 2002, 88. [Google Scholar] [CrossRef] [PubMed]
  19. Van Faassen, M.; de Boeij, P.L.; van Leeuwen, R.; Berger, J.A.; Snijders, J.G. Application of time-dependent current-density-functional theory to nonlocal exchange-correlation effects in polymers. J. Chem. Phys. 2003, 118. [Google Scholar] [CrossRef]
  20. Aulbur, W.G.; Jonsson, L.W.; Wilkins, J.W. Solid State Physics; Ehrenfeich, H., Spaepeu, F., Eds.; Academic Press: New York, NY, USA, 2000. [Google Scholar]
  21. Tao, J.; Vignale, G. Time-dependent density-functional theory beyond the local-density approximation. Phys. Rev. Lett. 2006, 97. [Google Scholar] [CrossRef] [PubMed]
  22. Conti, S.; Nifosi, R.; Tosi, M.P. The exchange-correlation potential for current-density functional theory of frequency-dependent linear response. J. Phys. Condens. Matter 1997, 9. [Google Scholar] [CrossRef]
  23. Nifosi, R.; Conti, S.; Tosi, M.P. Dynamic exchange-correlation potentials for the electron gas in dimensionality D = 3 and D = 2. Phys. Rev. B 1998, 58. [Google Scholar] [CrossRef]
  24. Conti, S.; Vignale, G. Elasticity of an electron liquid. Phys. Rev. B 1999, 60. [Google Scholar] [CrossRef]
  25. Landau, L.D.; Lifshitz, E.M. Lifshitz, Fluid Mechanics, 2nd ed.; Pergamon Press: Oxford, UK, 1987. [Google Scholar]
  26. Tao, J.; Gao, X.; Vignale, G.; Tokatly, I.V. Linear continuum mechanics for quantum many-body systems. Phys. Rev. Lett. 2009, 103. [Google Scholar] [CrossRef] [PubMed]
  27. Gao, X.; Tao, J.; Vignale, G.; Tokatly, I.V. Continuum mechanics for quantum many-body systems: Linear response regime. Phys. Rev. B 2010, 81. [Google Scholar] [CrossRef]
  28. Ceperley, D.M.; Alder, B.J. Ground state of the electron gas by a stochastic method. Phys. Rev. Lett. 1980, 45. [Google Scholar] [CrossRef]
  29. Ortiz, G.; Ballone, P. Correlation energy, structure factor, radial distribution function, and momentum distribution of the spin-polarized uniform electron gas. Phys. Rev. B 1994, 50. [Google Scholar] [CrossRef]
  30. Vosko, S.H.; Wilk, L.; Nusair, M. Accurate spin-dependent electron liquid correlation energies for local spin density calculations: A critical analysis. Can. J. Phys. 1980, 58, 1200–1211. [Google Scholar] [CrossRef]
  31. Perdew, J.P.; Wang, Y. Accurate and simple analytic representation of the electron-gas correlation energy. Phys. Rev. B 1992, 45. [Google Scholar] [CrossRef]
  32. Arfken, G.B.; Weber, H.J. Mathematical Methods for Physicists; Academic Press: Cambridge, MA, USA, 1985. [Google Scholar]
  33. Hasegawa, M.; Watabe, M.J. Theory of plasmon damping in metals. I. General formulation and application to an electron gas. Phys. Soc. Jpn. 1969, 27, 1393–1414. [Google Scholar] [CrossRef]
  34. Neilson, D.; Swierkowski, L.; Sj̎lander, A.; Szymanski, J. Dynamical theory for strongly correlated two-dimensional electron systems. Phys. Rev. B 1991, 44. [Google Scholar] [CrossRef]
  35. Qian, Z.; Vignale, G. Dynamical exchange-correlation potentials for an electron liquid. Phys. Rev. B 2002, 65. [Google Scholar] [CrossRef]
  36. Perdew, J.P.; Wang, Y. Pair-distribution function and its coupling-constant average for the spin-polarized electron gas. Phys. Rev. B 1992, 46. [Google Scholar] [CrossRef]
  37. Ullrich, C.A.; Burke, K.J. Excitation energies from time-dependent density-functional theory beyond the adiabatic approximation. Chem. Phys. 2004, 121. [Google Scholar] [CrossRef] [PubMed]
  38. Berger, J.A.; de Boeij, P.L.; van Leeuwen, R. Analysis of the viscoelastic coefficients in the Vignale-Kohn functional: The cases of one- and three-dimensional polyacetylene. Phys. Rev. B 2005, 71. [Google Scholar] [CrossRef]
  39. Tao, J.; Vignale, G. Analytic expression for the diamagnetic susceptibility of a uniform electron gas. Phys. Rev. B 2006, 74. [Google Scholar] [CrossRef]
  40. Tao, J.; Perdew, J.P.; Staroverov, V.N.; Scuseria, G.E. Climbing the density functional ladder: Nonempirical meta–generalized gradient approximation designed for molecules and solids. Phys. Rev. Lett. 2003, 91. [Google Scholar] [CrossRef] [PubMed]
  41. Seidl, M.; Perdew, J.P.; Kurth, S. Density functionals for the strong-interaction limit. Phys. Rev. A 2000, 62. [Google Scholar] [CrossRef]
  42. Tao, J.; Mo, Y. Accurate Semilocal Density Functional for Condensed-Matter Physics and Quantum Chemistry. Phys. Rev. Lett. 2016, 117. [Google Scholar] [CrossRef] [PubMed]
  43. Tokatly, I.V. Time-dependent deformation functional theory. Phys. Rev. B 2007, 75. [Google Scholar] [CrossRef]

Share and Cite

MDPI and ACS Style

Tao, J.; Vignale, G.; Zhu, J.-X. Geometric Derivation of the Stress Tensor of the Homogeneous Electron Gas. Computation 2017, 5, 28. https://doi.org/10.3390/computation5020028

AMA Style

Tao J, Vignale G, Zhu J-X. Geometric Derivation of the Stress Tensor of the Homogeneous Electron Gas. Computation. 2017; 5(2):28. https://doi.org/10.3390/computation5020028

Chicago/Turabian Style

Tao, Jianmin, Giovanni Vignale, and Jian-Xin Zhu. 2017. "Geometric Derivation of the Stress Tensor of the Homogeneous Electron Gas" Computation 5, no. 2: 28. https://doi.org/10.3390/computation5020028

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop