y}_{i}\frac{{\mathcal{G}}_{0}(y_{1},x_{2})({\mathcal{G}}_{0}(y_{2},x_{1})}{4|{\bf
x}_{1}-{\bf x}_{2}||{\bf y}_{1}-{\bf
y}_{2}|}\large[{\mathcal{G}}_{0}(x_{1},y_{1}){\mathcal{G}}_{0}(x_{2},y_{1})-{\mathcal{G}}_{0}(x_{1},y_{2}){\mathcal{G}}_{0}(x_{2},y_{1})\large]\;.$

And again the quantity inside the square brackets vanishes due to the same
reason as before. As for the cancellation of the last two graphs in Fig. 4, we
first use the bottom portion of (109) to obtain

$0=\frac{\delta}{\delta n}\left[\begin{picture}(136.0,20.0)(0.0,17.0) \raise
20.0pt\hbox to0.0pt{\kern 0.0pt\makebox(0.0,0.0)[lc]{\large$\frac{1}{2}$}\hss}
\ignorespaces \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \raise 20.0pt\hbox to0.0pt{\kern
84.0pt\makebox(0.0,0.0)[lc]{\large$-\frac{1}{2}$}\hss} \ignorespaces
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \end{picture}\right]=\begin{picture}(116.0,20.0)(0.0,17.0)
\raise 20.0pt\hbox to0.0pt{\kern 0.0pt\makebox(0.0,0.0)[lc]{\large$-$}\hss}
\ignorespaces \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\raise 20.0pt\hbox to0.0pt{\kern 60.0pt\makebox(0.0,0.0)[lc]{\large$+$}\hss}
\ignorespaces \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \end{picture}\;,$

and then

$0=\frac{1}{2}\left[\begin{picture}(116.0,20.0)(0.0,17.0) \raise 20.0pt\hbox
to0.0pt{\kern 0.0pt\makebox(0.0,0.0)[lc]{\large$-$}\hss} \ignorespaces
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \raise
20.0pt\hbox to0.0pt{\kern 60.0pt\makebox(0.0,0.0)[lc]{\large$+$}\hss}
\ignorespaces \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \end{picture}\right]\left(\begin{picture}(17.0,8.0)(0.0,0.0)
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\end{picture}\right)^{-1}\left[\begin{picture}(116.0,20.0)(0.0,17.0) \raise
20.0pt\hbox to0.0pt{\kern 0.0pt\makebox(0.0,0.0)[lc]{\large$-$}\hss}
\ignorespaces \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\raise 20.0pt\hbox to0.0pt{\kern 60.0pt\makebox(0.0,0.0)[lc]{\large$+$}\hss}
\ignorespaces \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \end{picture}\right]\;,$

or

$-\frac{1}{2}\;\;\begin{picture}(60.0,20.0)(78.0,17.0) \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){}
\end{picture}\;\;=\;\frac{1}{2}\;\;\begin{picture}(96.0,20.0)(0.0,-3.0)
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){}
\end{picture}\;\;-\;\begin{picture}(60.0,20.0)(0.0,-3.0) \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \end{picture}\;,$ (111)

when there is only one electron in the system.

Equation (111) means that at the single electron limit, the last two graphs in
Fig. 4 are equivalent to the last three graphs in Fig. 5. In other words, at
the single electron limit and at zero temperature, the set of order $U^{2}$
diagrams in $\Gamma[n]$ is equivalent to the set of order $U^{2}$ diagrams in
$\beta W[J_{0}]$. We shall pause at this point and elucidate the general
situation by looking at these cancellations via Hugenholtz diagrams. Negele
and Orland (1988)

With a two-body interaction term such as the Coulomb interaction, typical
Feynman diagrams treat the direct (Coulomb) and exchange matrix element
separately. It is not surprising that one may simplify the calculation by
combining the direct and exchange matrix elements into a single
antisymmetrized matrix element. The basic idea is to combine the following two
scenarios into one

$\displaystyle\begin{picture}(50.0,40.0)(-5.0,15.0) \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\raise-5.0pt\hbox to0.0pt{\kern-5.0pt\makebox(0.0,0.0)[lc]{$\alpha$}\hss}
\ignorespaces \raise-5.0pt\hbox to0.0pt{\kern
45.0pt\makebox(0.0,0.0)[rc]{$\rho$}\hss} \ignorespaces \raise 45.0pt\hbox
to0.0pt{\kern-5.0pt\makebox(0.0,0.0)[lc]{$\gamma$}\hss} \ignorespaces \raise
45.0pt\hbox to0.0pt{\kern 45.0pt\makebox(0.0,0.0)[rc]{$\delta$}\hss}
\ignorespaces\end{picture}-\begin{picture}(50.0,40.0)(-5.0,15.0)
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \raise-5.0pt\hbox
to0.0pt{\kern-5.0pt\makebox(0.0,0.0)[lc]{$\alpha$}\hss} \ignorespaces
\raise-5.0pt\hbox to0.0pt{\kern 45.0pt\makebox(0.0,0.0)[rc]{$\rho$}\hss}
\ignorespaces \raise 45.0pt\hbox
to0.0pt{\kern-5.0pt\makebox(0.0,0.0)[lc]{$\gamma$}\hss} \ignorespaces \raise
45.0pt\hbox to0.0pt{\kern 45.0pt\makebox(0.0,0.0)[rc]{$\delta$}\hss}
\ignorespaces\end{picture}$ $\displaystyle\equiv$
$\alpha$$\rho$$\gamma$$\delta$
$\displaystyle\left(\gamma\delta|\upsilon|\alpha\rho\right)-\left(\gamma\delta|\upsilon|\rho\alpha\right)$
$\displaystyle\equiv$ $\displaystyle\\{\gamma\delta|\upsilon|\alpha\rho\\}\;,$
(112)

where

$\left(\gamma\delta|\upsilon|\alpha\rho\right)\equiv\int
d\tau_{x}d\tau_{y}d{\bf x}d{\bf y}\phi_{\gamma}^{*}({\bf
x})\phi^{*}_{\delta}({\bf y})\upsilon(x,y)\phi_{\alpha}({\bf
x})\phi_{\rho}({\bf y})\;,$

and

$\left\\{\gamma\delta|\upsilon|\alpha\rho\right\\}\equiv\int
d\tau_{x}d\tau_{y}d{\bf x}d{\bf y}\phi_{\gamma}^{*}({\bf
x})\phi^{*}_{\delta}({\bf y})\upsilon(x,y)\left[\phi_{\alpha}({\bf
x})\phi_{\rho}({\bf y})-\phi_{\rho}({\bf x})\phi_{\alpha}({\bf y})\right]\;,$

with $\phi_{\alpha}({\bf x})$ being the single-particle wave function
described earlier in Eq. (93). The resulting diagrams with bullet dots as the
new vertices are called Hugenholtz diagrams. The appearance of those vertex
matrix elements comes from the following. When one evaluates a Feynman
diagram, a propagator ${\mathcal{G}}_{0}(x,x^{\prime})$ going from
$x^{\prime}$ to $x$ connects two vertices located at $x^{\prime}$ and $x$
respectively. From Eq. (95), we know that

${\mathcal{G}}_{0}(x,x^{\prime})=\sum_{\alpha}\phi_{\alpha}({\bf
x})\phi_{\alpha}^{*}({\bf
x}^{\prime})e^{-(\varepsilon_{\alpha}-\mu)(\tau_{x}-\tau_{x^{\prime}})}\left\\{\begin{array}[]{l
r}(-n_{\alpha})&{\rm if\leavevmode\nobreak\ }\tau_{x}\leq\tau_{x^{\prime}}\\\
(1-n_{\alpha})&{\rm if\leavevmode\nobreak\
}\tau_{x}>\tau_{x^{\prime}}\end{array}\right.\;.$

Upon integration of the space-time coordinates, we see that
$\phi_{\alpha}({\bf x})$ will be integrated with vertex $x$, whither the
propagator leads, while $\phi^{*}_{\alpha}({\bf x}^{\prime})$ will be
integrated with vertex $x^{\prime}$, whence the propagator originates.
Therefore, when evaluating a Feynman diagram, one may associate each
propagator going from time $\tau^{\prime}$ to time $\tau$ with a state index
$\alpha$ with

${\mathcal{G}}_{0}(\alpha,\tau-\tau^{\prime})=e^{-(\varepsilon_{\alpha}-\mu)(\tau-\tau^{\prime})}\left[(1-n_{\alpha})\theta(\tau-\tau^{\prime}-\eta)-n_{\alpha}\theta(\tau^{\prime}-\tau+\eta)\right]\;.$

Each vertex will contribute a numerical factor equivalent to its vertex matrix
element. Note that each vertex carries a time index. At zero temperature, for
a vertex with time index $\tau$, if the two incoming lines originate from
vertices with times larger than or equal to $\tau$, the vertex matrix element
associated with $\tau$ becomes zero in the single electron limit. This is
because both incoming lines each carry only the $\alpha=1$ index. Upon
antisymmetrization, the Hugenholtz vertex matrix element becomes zero.

For an arbitrary Hugenholtz diagram with $n$ vertices, one may always name
their time index such that
$\tau_{1}\leq\tau_{2}\leq\tau_{3}\leq\cdots\leq\tau_{n}$. In this case, the
vertex associated with $\tau_{1}$ gives zero matrix element since its two
incoming lines must come from two other time indices that are larger than or
equal to $\tau_{1}$. Consequently, each Hugenholtz diagram yields value zero
at zero temperature when there is only one electron present. As a matter of
fact, the order $U$ diagrams of $W[J_{0}]$ are represented by the Hugenholtz
diagram \begin{picture}(24.0,6.0)(0.0,-3.0) \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \end{picture} . On the other hand, the first two graphs of
Fig. 4 (order $U^{2}$ terms of $\Gamma[n]$) correspond to the Hugenholtz
diagram \begin{picture}(20.0,5.0)(0.0,-3.0) \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \end{picture}
, while the last two graphs of Fig. 4 (order $U^{2}$ terms of $\Gamma[n]$) are
equivalent to the last three graphs of Fig. 5 (order $U^{2}$ terms of $\beta
W[J_{0}]$) and correspond to the Hugenholtz diagram
\begin{picture}(38.0,10.0)(0.0,-3.0) \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \end{picture}
(see reference Negele and Orland, 1988). Therefore, we see that all the order
$U^{2}$ terms cancel each other out in the effective action $\Gamma[n]$.

If one were to expand the effective action (80) in powers of $U$, our approach
reduces to performing the inversion method using $e^{2}$ as the expansion
parameter.Valiev and Fernando (1997); Fukuda et al. (1995) In this case,
$W_{l\geq 1}[J_{0}]$ in Eqs. (47-73) contains exactly all the order $U^{l}$
diagrams in $W[J_{0}]$ and can be expressed as Hugenholtz diagrams of order
$U^{l}$ as well. Since all the Hugenholtz diagrams give value zero, all the
derivatives of $W_{l}$ (all the order $u^{l}$ diagrams) with respect to the
density vanish as well. This implies that all the $J_{l\geq 1}$ vanish and
consequently in the effective action the sum of Hartree energy and the
exchange-correlation energy equals zero at zero temperature when there is only
one electron present.

The perfect cancellation of the Hartree energy and exchange-correlation energy
means that one is the negative of the other. This means that
$\sum_{i=1}^{\infty}\Gamma_{i}=-\frac{1}{2}n{\scriptstyle\circ}\,U{\scriptstyle\circ}n$.
Employing Eq. (77) for the most general case, we obtain

$\tilde{J}_{0}-J=\frac{\delta\left(\sum_{i=1}\Gamma_{i}[n]\right)}{\delta
n}=-\frac{1}{2}\frac{\delta\left(n{\scriptstyle\circ}\,U{\scriptstyle\circ}n\right)}{\delta
n}=-U{\scriptstyle\circ}n\;.$

Since $J_{0}=\tilde{J}_{0}+U{\scriptstyle\circ}n$, we find that $J_{0}=J$ at
the single electron limit at zero temperature. The fact that $J_{0}=J$ means
that the final effective potential $\varepsilon_{1}-\mu-J{\cdot}n$ is nothing
but the lowest eigenenergy of the following single-particle Hamiltonian

$\hat{h}({\bf x})=-\frac{\nabla^{2}}{2m}+\upsilon_{\rm ion}({\bf x})-\mu+J$

less the expectation value of $J$, which is exactly what one expected.

###  IV.4 Screening of Coulomb Interaction

In the physical limit, $n=i\varphi$. In our formulation, lumps of charge
fluctuation around the configuration $i\varphi$ interact with each other via

$U{\scriptstyle\circ}\tilde{\cal
D}_{0}^{-1}{\scriptstyle\circ}\,U=\left(U-U{\scriptstyle\circ}D_{0}{\scriptstyle\circ}\,U\right)\;.$

As will be described in the discussion section,
$U{\scriptstyle\circ}(i\phi)=ib$ plays the role of the photon field here.
Therefore,

$-\phi{\scriptstyle\circ}\left(U-U{\scriptstyle\circ}D_{0}{\scriptstyle\circ}\,U\right){\scriptstyle\circ}\phi=(i\phi{\scriptstyle\circ}\,U){\scriptstyle\circ}\tilde{\mathcal{D}}_{0}^{-1}{\scriptstyle\circ}(U{\scriptstyle\circ}i\phi)=(ib){\scriptstyle\circ}\left(U^{-1}-D_{0}\right){\scriptstyle\circ}(ib)\;.$

Since $U(x,y)=\delta(\tau_{x}-\tau_{y})U({\bf x},{\bf y})$ and
${\mathcal{G}}_{0}(x,y)$ only depends on the relative time difference
$\tau_{x}-\tau_{y}$, we expect $D_{0}(x,y)$ to depend only on
$\tau\equiv\tau_{x}-\tau_{y}$. Furthermore, since ${\mathcal{G}}_{0}(x,y)$ is
antiperiodic in $\tau_{x}$, $\tau_{y}$, and $\tau_{x}-\tau_{y}$, one expects
that ${\mathcal{G}}_{0}(x,y){\mathcal{G}}_{0}(y,x)$ to be periodic in
$\tau_{x}-\tau_{y}$.

Introducing the spatial momenta ${\mathbf{p}}$ and ${\mathbf{q}}$ conjugate to
the spatial variables ${\bf x}$ and ${\bf y}$, we may write

$\left[U^{-1}-D_{0}\right](\nu_{n},{\mathbf{p}},{\mathbf{q}})=\int
e^{i{\mathbf{p}}{\bf x}+i{\mathbf{q}}{\bf y}}d{\bf x}d{\bf
y}\int_{0}^{\beta}\,d\tau\;e^{i\nu_{n}\tau}\left[U^{-1}({\bf x},{\bf
y})\delta(\tau)-{\mathcal{G}}_{0}(x,y){\mathcal{G}}_{0}(y,x)\right]\;.$

When $U({\bf x},{\bf y})={e^{2}\over|{\bf x}-{\bf y}|}$, we find that
$U^{-1}(x,y)=-{1\over 4\pi e^{2}}\delta(\tau_{x}-\tau_{y})\nabla^{2}_{\bf
x}\delta({\bf x}-{\bf y})$ and

$U^{-1}(\nu_{n},{\mathbf{p}},{\mathbf{q}})={(2\pi)^{3}\over 4\pi
e^{2}}{\mathbf{q}}^{2}\delta({\mathbf{p}}+{\mathbf{q}})={L^{3}\over 4\pi
e^{2}}{\mathbf{q}}^{2}\delta_{{\mathbf{p}}+{\mathbf{q}}}\;,$

where $L^{3}$ is the spatial volume of the system. Let us now write
$D_{0}(\nu_{n},{\mathbf{p}},{\mathbf{q}})$ as

$D_{0}(\nu_{n},{\mathbf{p}},{\mathbf{q}})=\int e^{i{\mathbf{p}}{\bf
x}+i{\mathbf{q}}{\bf y}}d{\bf x}d{\bf
y}\int_{0}^{\beta}\,d\tau\;e^{i\nu_{n}\tau}\,{\mathcal{G}}_{0}(x,y){\mathcal{G}}_{0}(y,x)\;.$

In momentum space,
$U^{-1}\propto{{\mathbf{q}}}^{2}\delta({\mathbf{p}}+{\mathbf{q}})$ and it is
well known that this type of Coulomb interaction leads to infrared divergence
(that occurs near ${{\mathbf{q}}}\to 0$). In the presence of $-D_{0}$, one may
wish to calculate the zero momentum contribution of $-D_{0}$. Let us define
$D_{0}(\nu_{n})\equiv D_{0}(\nu_{n},{\mathbf{p}}=0,{\mathbf{q}}=0)$. Using Eq.
(95), we find

$D_{0}(\nu_{n})=-\sum_{\alpha,\alpha^{\prime}}\delta_{\alpha,\alpha^{\prime}}\int_{0}^{\beta}d\tau
e^{i\nu_{n}\tau}\left[n_{\alpha}(1-n_{\alpha^{\prime}})\right]=-\beta\delta_{\nu_{n},0}\sum_{\alpha}n_{\alpha}(1-n_{\alpha})\;.$
(113)

From Eq. (101), we see that when $x=x^{\prime}$,

$\beta\sum_{\alpha}n_{\alpha}(1-n_{\alpha})=-\int d{\bf
x}\frac{\partial{\mathcal{G}}_{0}(x,x)}{\partial\mu}=\frac{\partial
N_{e}}{\partial\mu}\;.$

We therefore have

$D_{0}(\nu_{n})=-\delta_{\nu_{n},0}{L^{3}}\frac{\partial\bar{n}}{\partial\mu}\;,$

where $\bar{n}=N_{e}/L^{3}$ is the average electron density.

Therefore at the stationary limit, where $\nu_{n}=0$, the inverse propagator
with nearly zero momentum behaves as

$L^{3}\left[\frac{q^{2}}{4\pi
e^{2}}+\frac{\partial\bar{n}}{\partial\mu}\right]\Rightarrow\frac{L^{3}}{4\pi
e^{2}}\left[q^{2}+4\pi e^{2}\frac{\partial\bar{n}}{\partial\mu}\right]\;,$

analogous to the Thomas-Fermi results for electric charge screening.

Except during the intermediate steps of computing the excitation spectrum, one
deals with the time independent system. Note that $4\pi
e^{2}\partial\bar{n}/\partial\mu$ plays the role of $m^{2}$ in the screened
potential $e^{-mr}/r$. This shows that the static interaction between charge
fluctuations in the long wave length limit is a screened one instead of the
bare Coulomb interaction. It should be noted that the use of this screened
propagator dates back half a century. DuBois DuBois (1959) replaced the bare
Coulomb interaction by the screened interaction in his study of electron
interactions. Hedin Hedin (1965) also used it to replace the bare Coulomb
interaction in the expansion of the Luttinger-Ward functional, Luttinger and
Ward (1960) resulting in the so-called GW approximation. The differences among
the three mentioned approaches should be described. In both DuBois’s and
Hedin’s formalism, they use the full polarization of the interacting system.
The difference between their formalism is in the electron propagators
employed. In DuBois’s approach, the free electron propagator is used, while in
Hedin’s method, similar to the Luttinger-Ward functional, Luttinger and Ward
(1960) the full electron propagator is employed. In our approach, it is the
polarization of the KS system that enters the calculation and the electron
propagator entering the diagrammatic calculation is the noninteracting KS
Green’s function.

###  IV.5 Homogeneous Electron Gas

Being the foundation for the LDA of the DFT, the homogeneous electron gas
(HEG) is also a simplified model for a metal or a plasma. The HEG system is an
interacting electron gas placed in a uniformly distributed positive background
chosen to ensure that the total system is neutral. Due to the translational
symmetry of the HEG system, the one-particle propagator functions only depend
on the coordinate difference between two variables (space-time points) instead
of on both variables. Consequently, each propagator (Green’s function) carries
a definite momentum even if the Coulomb interaction among electrons is fully
taken into consideration.

Investigation of the HEG system dates back to 1930s by Wigner,Wigner (1934)
who coined the term “correlation energy” to represent the ground state energy
per electron after subtracting away the average kinetic and exchange energy.
(The Hartree energy is cancelled exactly by the interaction between electrons
and the positive background, and by the Coulomb energy among the uniform
positive charges.) Apparently, after choosing the Rydberg (the negative of the
Hydrogen atom ground state energy) as the energy unit, one may express either
the correlation energy or the ground state energy of the HEG system in terms
of the dimensionless quantity $r_{s}\equiv r_{0}/r_{B}$, where ${4\pi\over
3}r_{0}^{3}=\frac{V}{N_{e}}$, $V$ being the volume of the system considered,
$N_{e}$ being the total number of electrons, and
$r_{B}=\frac{\hbar^{2}}{me^{2}}$ representing the Bohr radius. In fact,
efforts have been invested to express the correlation energy as a power series
in $r_{s}$ (and $r_{s}\ln r_{s}$ as well) at both the high density and low
density limits. Lacking real experimental results to compare with, however, it
is hard to assess how much improvement each increasing order of $r_{s}$ (or
$1/r_{s}$) can bring.

Since $r_{s}\propto 1/r_{B}\propto e^{2}$, a small $r_{s}$ expansion is most
naturally performed by treating the Coulomb interaction as a perturbation.
Equivalent to the method of computing the vacuum amplitude,Fetter and Walecka
(1971); Jones and March (1973); Negele and Orland (1988) the time-ordered
perturbative expansion developed by Goldstone Goldstone formalized the
Brueckner theory Brueckner and Levinson (1955) and made a direct connection to
diagrammatic expansion in calculation of the ground state energy. The electron
propagator in diagrams under this formalism is that of the noninteracting
electrons. An alternative formalism to compute the ground state energy (or the
grand potential) is to utilize the full propagator (i.e., the one with self-
energy included) in diagrammatic calculation. Under this alternative
formalism, only two-particle irreducible diagrams contribute. Indeed, the work
of Luttinger and Ward Luttinger and Ward (1960) (as well as that of Klein
Klein (1961)) was exactly along this line. Since the self-energy is not known
a priori, this type of energetic expression is a functional of the self-energy
(or the full Green’s function), which must be determined via a stationary
condition.Luttinger and Ward (1960)

An equivalent of the Brueckner-Goldstone formalism Goldstone ; Brueckner and
Gammel (1958) was used by Gell-Mann and Brueckner Gell-Mann and Brueckner
(1957) to compute the correlation energy of the HEG system. Given the long-
range nature of the Coulomb interaction, it is not surprising that Gell-Mann
and Brueckner identified the occurrence of divergence as early as in the
second order of $e^{2}$-based perturbative calculation. To circumvent this
unphysical divergence, they summed an infinite subset of diagrams to arrive at
a contribution proportional to $\ln r_{s}$. To eliminate unphysical
divergence, DuBois DuBois (1959) replaced the bare Coulomb interaction by the
screened one. Starting with calculating the vacuum amplitude, he expressed the
ground state energy in terms of integration of the full electronic
polarization over the Coulomb coupling strength, using a version of the
Feynman-Hellmann theorem. While the full polarization is used, the electron
propagator entering DuBois’s formalism is the free electron propagator instead
of the full propagator. Unfortunately, DuBois made a mistake in extracting the
higher order terms of the ground state energy at the high density (small
$r_{s}$) limit. Although this error was later corrected by Carr and
Maradudin,Carr and Maradudin (1964) according to Hedin,Hedin (1965) this high
density expansion violated at moderate $r_{s}$ values Ferrell’s
condition,Ferrell (1958) which is based on the simple fact that the second
order perturbation contribution to the ground state energy is always negative.

Except for using the screened Coulomb interaction to replace the bare one and
working directly at zero temperature, Hedin’s approach is largely similar to
the finite-temperature formalism of Luttinger and Ward Luttinger and Ward
(1960) in the sense that the full electron propagator is used as the
fundamental variable. As a consequence, the higher order diagrams contributing
to the full polarization appear different in DuBois’s work and in Hedin’s
formalism. Instead of solving his own self-consistent equations, however,
Hedin approximated the full Green’s function within his formalism by the non-
interacting Green’s function to compute the energetics of the HEG system.
Nevertheless, it is still possible to maintain the exactness within Hedin’s
approach if two-particle reducible diagrams are properly included. We will
provide an example of such a two-particle reducible diagram that would not be
included in Hedin’s approximate calculation but should be incorporated for
theoretical soundness.

Another important issue in obtaining the ground state energy of a many-
electron system has to do with whether the finite temperature formalism
(setting ${V\to\infty}$ first followed by ${T\to 0}$) or the zero temperature
formalism (setting $T\to 0$ and then followed by $V\to\infty$) is used. In
terms of diagrammatic expansion, there exist diagrams (termed anomalous
diagrams by Luttinger et al.Kohn and Luttinger (1960); Luttinger and Ward
(1960)) that are present (giving finite contribution) within the finite
temperature formalism but are absent (giving zero contribution) within the
zero temperature formalism. To be specific, a diagram is anomalous if within
it there exist two electron propagators linking two different times (say
$\tau_{1}$, $\tau_{2}$ and of course $\tau_{1}\neq\tau_{2}$) in the opposite
orders: ${\mathcal{G}}({\bf x},\tau_{1},{\bf
y},\tau_{2}){\mathcal{G}}({\mathbf{w}},\tau_{2},{\mathbf{z}},\tau_{1})$.

Zero temperature methods typically rely on the Gell-Mann and Low theorem,
Gell-Mann and Low (1951) which assures that adiabatic transformation of the
noninteracting ground state by gradually switching on the interaction leads to
an eigenstate of the interacting system.Fetter and Walecka (1971); Negele and
Orland (1988) The energetic difference of this eigenstate and the
noninteracting ground state is what the zero temperature formalism obtains.
This approach therefore makes sense only if the adiabatically transformed
state is the ground state of the interacting system, which occurs only if the
noninteracting Fermi surface is identical to that of the interacting system.
Jones and March (1973) For the HEG system, the perturbed Fermi surface remains
spherical (identical to that of the unperturbed one) since the Coulomb
interaction respects spherical symmetry and the background positive charge
distribution also has spherical symmetry. Consequently, for the HEG system,
the anomalous diagrams should end up giving no contribution. Indeed, Luttinger
et al. Kohn and Luttinger (1960); Luttinger and Ward (1960) illustrated how
the contributions from the anomalous diagrams in this case are cancelled by
the chemical potential shift. Luttinger and Ward Luttinger and Ward (1960)
also showed how one may avoid anomalous contributions from appearing by
expressing the grand potential as a sum of all possible linked diagrams
(including two-particle reducible ones). The key step there is to subtract
from each self-energy part a number, which is given by that self-energy part
evaluated at the Fermi surface. As we will illustrate later, under our finite
temperature formalism it is not necessary to implement such an elaborate
subtraction scheme because each two-particle reducible diagram is
automatically accompanied by another appropriate diagrammatic subtraction.

What sets our self-consistent equation (78) apart from that of Luttinger and
Ward and of Hedin is the variable to be solved for. It is the KS potential,
instead of the physical Green’s function, that enters our exact, self-
consistent equation. In general, one needs to first solve $\tilde{J}_{0}$
self-consistently using eq. (78) prior to the evaluation of the grand
potential (or the ground state energy). When limiting to a homogeneous system
with constant electron density, however, $\tilde{J}_{0}={\rm const.}$ becomes
the only possibility. Whether or not such a choice satisfies the self-
consistent equation (78) depends on whether a uniform electron density truly
represent the lowest energy configuration. For the present purpose of
considering a high density HEG, we assume a constant $\tilde{J}_{0}$ to
proceed. Because a constant $\tilde{J}_{0}$ can be easily absorbed into the
chemical potential, the resulting KS Green’s function carries a corresponding
energy that consists of kinetic energy only. That is, rather than being an
approximation as in Hedin’s case, the single-particle Green’s function
carrying only kinetic energy represents exactly the self-consistent KS Green’s
function under our formalism. Consequently, for the HEG system, the grand
potential (or the ground state energy) may be calculated using eq. (80). In
the remaining part of this section, we will first show how our formalism
naturally avoids divergence and how one may use it to obtain the ground state
energy of the HEG as a series in $r_{s}$ (and $\ln r_{s}$). We will then
illustrate with an explicit example how the anomalous contributions are
cancelled within our formalism, followed by a brief description of diagrams
that would be missed within Hedin’s approximation Hedin (1965) in computing
the ground state energy of the HEG system. To provide an easier comparison
with existing results, we restore in this section the electron spins that have
been suppressed thus far to simplify the exposition.

In our definition of the energy functional $E_{\upsilon}$, see eq. (4), the
$-\mu N_{e}$ term is included but the interaction among background charges is
not included. Since most zero temperature formalism does not include the $-\mu
N_{e}$ term and for the HEG system the interaction between background charges
is also included, the ground state energy in the literature will correspond to

$\lim_{\beta\to\infty}\left[\frac{1}{\beta}\Gamma[n]+\mu
N_{e}+\frac{1}{2}n_{\rm bg}{\cdot}\,U{\cdot}n_{\rm bg}\right]$

in our formalism. Eq. (104) can then be used to arrive at

$\displaystyle E_{g}$ $\displaystyle=$
$\displaystyle\lim_{\beta\to\infty}\left[\frac{1}{\beta}\Gamma[n]+\mu
N_{e}+\frac{1}{2}n_{\rm bg}{\cdot}\,U{\cdot}n_{\rm bg}\right]$ (114)
$\displaystyle=$
$\displaystyle\sum_{\alpha=1}^{N_{e}}\varepsilon_{\alpha}-J_{0}{\cdot}n+\frac{1}{2}n{\cdot}\,U{\cdot}n+\frac{1}{2}n_{\rm
bg}{\cdot}\,U{\cdot}n_{\rm
bg}+\lim_{\beta\to\infty}\left[\frac{1}{2\beta}\text{Tr}\ln\left(\tilde{\cal
D}_{0}^{-1}{\scriptstyle\circ}\,U\right)+\frac{1}{\beta}\sum_{i=2}^{\infty}\Gamma_{i}[n]\right]$
$\displaystyle=$
$\displaystyle\sum_{\alpha=1}^{N_{e}}(\varepsilon_{\alpha}-\tilde{J}_{0})+\lim_{\beta\to\infty}\left[\frac{1}{2\beta}\text{Tr}\ln\left(\tilde{\cal
D}_{0}^{-1}{\scriptstyle\circ}\,U\right)+\frac{1}{\beta}\sum_{i=2}^{\infty}\Gamma_{i}[n]\right]\;,$

where $J_{0}=\tilde{J}_{0}+U{\cdot}n$ is used and $n_{\rm bg}=n$ for the HEG
system is also employed. Note that the state label
$\alpha=({\mathbf{p}},\sigma)$ includes both momentum and spin. For the HEG
system, $\upsilon_{\rm ion}+U{\cdot}n=0$ and thus eq. (94) leads to
$\varepsilon_{\alpha}={\mathbf{p}}^{2}/2m+\tilde{J}_{0}$. Therefore, the
ground state energy of the HEG system may be written as

$E_{g}=2\sum_{{\mathbf{p}}}^{|{\mathbf{p}}|\leq
p_{F}}\frac{{\mathbf{p}}^{2}}{2m}+\lim_{\beta\to\infty}\frac{1}{\beta}\left[\frac{1}{2}\text{Tr}\ln\left(\tilde{\cal
D}_{0}^{-1}{\scriptstyle\circ}\,U\right)+\sum_{i=2}^{\infty}\Gamma_{i}[n]\right]\;,$
(115)

where $p_{F}$ indicates the Fermi momentum and the factor of $2$ in the first
term of the right hand side comes from noting that there are two spin states
associated with each momentum. Furthermore, the KS Green’s function in this
case may be written as

${\mathcal{G}}_{0}(x,\sigma;y,\sigma^{\prime})=\frac{1}{\beta
V}\sum_{{\mathbf{p}},n}\frac{e^{-i\omega_{n}(\tau_{x}-\tau_{y})}e^{-i{\mathbf{p}}\cdot({\bf
x}-{\bf
y})}}{-i\omega_{n}+\varepsilon_{\mathbf{p}}+\tilde{J}_{0}-\mu}\,\delta_{\sigma,\sigma^{\prime}}\,,$
(116)

where $\varepsilon_{\mathbf{p}}={\mathbf{p}}^{2}/2m$ is the kinetic energy of
an electron carrying momentum ${\mathbf{p}}$. By absorbing the constant
$\tilde{J}_{0}$ into the chemical potential in (116), one sees that the KS
Green’s function in the HEG system is indeed the free electron propagator and
at the zero temperature, the new chemical potential (the original one
subtracted by $\tilde{J}_{0}$) simply becomes $p_{F}^{2}/2m$.

The diagrammatic expression of our first correction term
$\Gamma_{1}[n]=\frac{1}{2}\text{Tr}\ln(\tilde{\mathcal{D}_{0}}^{-1}{\scriptstyle\circ}\,U)=\frac{1}{2}\text{Tr}\ln\left[{\mathbf{I}}-D_{0}{\scriptstyle\circ}\,U\right]$
is shown in Fig. 1. It can be seen that, upon being divided by the inverse
temperature $\beta$, our $\Gamma_{1}[n]$ contains the exchange energy
$\epsilon_{x}$ and all the ring-like correlation energy $\epsilon^{\prime}$
discussed by Gell-Mann and Brueckner Gell-Mann and Brueckner (1957) via the
relation $\Gamma_{1}[n]/\beta=N_{e}(\epsilon_{x}+\epsilon^{\prime})$, with
$N_{e}$ being the total number of electrons. Since both $D_{0}$ and $U$ are
diagonal in momentum space for the HEG system, we may write

$\frac{1}{\beta}\Gamma_{1}[n]=\frac{1}{2\beta}\sum_{{\mathbf{q}},\nu_{n}}\ln\left[1-D_{0}(q,\nu_{n})\,U(q)\right]=\frac{V}{2\beta}\int\frac{d{\mathbf{q}}}{(2\pi)^{3}}\sum_{\nu_{n}}\ln\left[1-D_{0}(q,\nu_{n})\,U(q)\right]\;.$
(117)

By comparing this equation to a derived result (eq. (30.16) of reference
Fetter and Walecka, 1971), it is also evident that our $\Gamma_{1}[n]/(\beta
N_{e})$ indeed gives $\epsilon_{x}+\epsilon^{\prime}$ of reference Gell-Mann
and Brueckner, 1957 as $\beta\to\infty$. In Fig. 1, except for the first
diagram, all other diagrams when evaluated individually exhibit infrared
divergence due to the piling up of $1/{\mathbf{q}}^{2}$ propagators.Gell-Mann
and Brueckner (1957) To obtain the leading contribution, one must pay
particular attention to the small $|{\mathbf{q}}|$ region.

The polarization $D_{0}\equiv\delta n(x)/\delta J_{0}(y)$ is defined in eq.
(71). When the electron spins are included,
$n(x)=-\sum_{\sigma}{\mathcal{G}}_{0}(x,\sigma;x,\sigma)$. Because the Coulomb
interaction does not flip spins, the KS propagator is diagonal in the electron
spins and thus

$D_{0}(x,y)=\sum_{\sigma}{\mathcal{G}}_{0}(x,\sigma;y,\sigma){\mathcal{G}}_{0}(y,\sigma;x,\sigma)=\,2\;{\mathcal{G}}_{0}(x,y){\mathcal{G}}_{0}(y,x)\;.$
(118)

We therefore have

$\displaystyle D_{0}({\mathbf{q}},\nu_{n})$ $\displaystyle\equiv$
$\displaystyle 2\int d({\bf x}-{\bf
y})d(\tau_{x}-\tau_{y})e^{i{\mathbf{q}}.({\bf x}-{\bf
y})}d^{i\nu_{n}(\tau_{x}-\tau_{y})}\left[{\mathcal{G}}_{0}(x,y){\mathcal{G}}_{0}(y,x)\right]$
$\displaystyle=$
$\displaystyle\frac{2}{\beta}\sum_{n^{\prime}}\int\frac{d{\mathbf{p}}}{(2\pi)^{3}}\frac{1}{(-i\omega_{n^{\prime}}-i\nu_{n}+\varepsilon_{{\mathbf{p}}+{\mathbf{q}}}-\mu)(-i\omega_{n^{\prime}}+\varepsilon_{\mathbf{p}}-\mu)}$
$\displaystyle=$ $\displaystyle
2\int\frac{d{\mathbf{p}}}{(2\pi)^{3}}\frac{n_{{\mathbf{p}}+{\mathbf{q}}}-n_{{\mathbf{p}}}}{-i\nu_{n}+\varepsilon_{{\mathbf{p}}+{\mathbf{q}}}-\varepsilon_{{\mathbf{p}}}}$
$\displaystyle=$
$\displaystyle-2\int\frac{d{\mathbf{p}}}{(2\pi)^{3}}\frac{n_{\mathbf{p}}(1-n_{{\mathbf{p}}+{\mathbf{q}}})-n_{{\mathbf{p}}+{\mathbf{q}}}(1-n_{{\mathbf{p}}})}{-i\nu_{n}+\varepsilon_{{\mathbf{p}}+{\mathbf{q}}}-\varepsilon_{{\mathbf{p}}}}$
$\displaystyle=$
$\displaystyle-2\int\frac{d{\mathbf{p}}}{(2\pi)^{3}}n_{\mathbf{p}}(1-n_{{\mathbf{p}}+{\mathbf{q}}})\left[\frac{1}{-i\nu_{n}+\varepsilon_{{\mathbf{p}}+{\mathbf{q}}}-\varepsilon_{{\mathbf{p}}}}+\frac{1}{i\nu_{n}+\varepsilon_{{\mathbf{p}}+{\mathbf{q}}}-\varepsilon_{{\mathbf{p}}}}\right]$
$\displaystyle=$
$\displaystyle-2\int\frac{d{\mathbf{p}}}{(2\pi)^{3}}\;n_{\mathbf{p}}\left[\frac{1}{-i\nu_{n}+\varepsilon_{{\mathbf{p}}+{\mathbf{q}}}-\varepsilon_{{\mathbf{p}}}}+\frac{1}{i\nu_{n}+\varepsilon_{{\mathbf{p}}+{\mathbf{q}}}-\varepsilon_{{\mathbf{p}}}}\right]$
$\displaystyle=$
$\displaystyle-4\int\frac{d{\mathbf{p}}}{(2\pi)^{3}}\frac{n_{\mathbf{p}}(1-n_{{\mathbf{p}}+{\mathbf{q}}})(\varepsilon_{{\mathbf{p}}+{\mathbf{q}}}-\varepsilon_{\mathbf{p}})}{(\varepsilon_{{\mathbf{p}}+{\mathbf{q}}}-\varepsilon_{\mathbf{p}})^{2}+\nu_{n}^{2}}=-4\int\frac{d{\mathbf{p}}}{(2\pi)^{3}}\frac{n_{\mathbf{p}}(\varepsilon_{{\mathbf{p}}+{\mathbf{q}}}-\varepsilon_{\mathbf{p}})}{(\varepsilon_{{\mathbf{p}}+{\mathbf{q}}}-\varepsilon_{\mathbf{p}})^{2}+\nu_{n}^{2}}$

where the third line of the above equation is obtained using the technique
described in section III.5, the fourth line is obtained by rewriting the
numerator of the third line as
$-n_{\mathbf{p}}(1-n_{{\mathbf{p}}+{\mathbf{q}}})+n_{{\mathbf{p}}+{\mathbf{q}}}(1-n_{\mathbf{p}})$,
the fifth line is obtained by a change of variable
(${\mathbf{p}}+{\mathbf{q}}\to-{\mathbf{p}}$) in the second half and assuming
the spherical property of $n_{\mathbf{p}}$, and the sixth line comes from the
fact that $n_{\mathbf{p}}n_{{\mathbf{p}}+{\mathbf{q}}}$ is symmetric with
respect to $({\mathbf{p}}+{\mathbf{q}})\leftrightarrow{\mathbf{p}}$ while the
quantity inside the square bracket is antisymmetric with respect to
$({\mathbf{p}}+{\mathbf{q}})\leftrightarrow{\mathbf{p}}$. Since
$\varepsilon_{{\mathbf{p}}}\propto|{\mathbf{p}}|^{2}$ is a monotonic
increasing function of $|{\mathbf{p}}|$, while both
$n_{{\mathbf{p}}}=1/(e^{\beta(\varepsilon_{\mathbf{p}}-\mu)}+1)$ and
$1-n_{\mathbf{p}}$ are nonnegative, we see that $D_{0}({\mathbf{q}},\nu_{n})$
is a negative definite quantity.

It turns out that for the HEG system one can further simplify the expression
of $D_{0}(q,\nu)$ by introducing a new variable $u$ via $\nu\equiv
u|{\mathbf{q}}|/m$. We then have (by choosing $\hat{q}$ as the $\hat{z}$
direction in ${\mathbf{p}}$ and using the fact that
$n_{\mathbf{p}}=n(|{\mathbf{p}}|)=n(p)$)

$\displaystyle D_{0}(q,\frac{uq}{m})$ $\displaystyle=$
$\displaystyle-2\int\frac{d{\mathbf{p}}}{(2\pi)^{3}}\;n(p)\left[\frac{1}{-iuq/m+\varepsilon_{{\mathbf{p}}+{\mathbf{q}}}-\varepsilon_{{\mathbf{p}}}}+\frac{1}{iuq/m+\varepsilon_{{\mathbf{p}}+{\mathbf{q}}}-\varepsilon_{{\mathbf{p}}}}\right]$
(120) $\displaystyle=$
$\displaystyle-2\frac{m}{q}\int_{0}^{\infty}\frac{n(p)\,p^{2}\,dp}{(2\pi)^{2}}\;\int_{-1}^{1}d\cos\theta\left[\frac{1}{-iu+q/2+p\cos\theta}+\frac{1}{iu+q/2+p\cos\theta}\right]$
$\displaystyle=$
$\displaystyle-2\frac{m}{q}\int_{0}^{\infty}\frac{n(p)\,dp}{(2\pi)^{2}}\;p\left[\ln\left(\frac{q}{2}+p+iu\right)+\ln\left(\frac{q}{2}+p-iu\right)\right.$
$\displaystyle\hskip
90.0pt\left.-\ln\left(\frac{q}{2}-p+iu\right)-\ln\left(\frac{q}{2}-p-iu\right)\right]$
$\displaystyle=$
$\displaystyle\frac{2m}{(2\pi)^{2}}\int_{0}^{\infty}\frac{\partial
n(p)}{\partial
p}\,dp\,\left[p+\frac{1}{2q}(p^{2}+u^{2}-\frac{q^{2}}{4})\ln\frac{(p+\frac{q}{2})^{2}+u^{2}}{(p-\frac{q}{2})^{2}+u^{2}}\right.$
$\displaystyle\hskip
115.0pt\left.-u\tan^{-1}\frac{\frac{q}{2}+p}{u}+u\tan^{-1}\frac{\frac{q}{2}-p}{u}\right]\;.$

As $\beta\to\infty$, $n(p)=\theta(p_{F}-p)$ and therefore $\partial
n(p)/\partial p=-\delta(p-p_{F})$, leading to

$\displaystyle D_{0}(q,\frac{uq}{m})$
$\displaystyle\xrightarrow[\beta\to\infty]{}$
$\displaystyle-\frac{2m}{(2\pi)^{2}}\left[p_{F}+\frac{1}{2q}(p_{F}^{2}+u^{2}-\frac{q^{2}}{4})\ln\frac{(p_{F}+\frac{q}{2})^{2}+u^{2}}{(p_{F}-\frac{q}{2})^{2}+u^{2}}\right.$
(121) $\displaystyle\hskip
75.0pt\left.-u\tan^{-1}\frac{\frac{q}{2}+p_{F}}{u}+u\tan^{-1}\frac{\frac{q}{2}-p_{F}}{u}\right]\;.$

The exact expression (121) allows one to extract the limits of $q\to\infty$
and $q\to 0$, both of which are important for determining the convergence
properties of the energy expansion. We have

$\displaystyle D_{0}(q\gg 1,\frac{uq}{m})$
$\displaystyle\xrightarrow[\beta\to\infty]{}$
$\displaystyle-\frac{4m\,p_{F}^{3}}{3\pi^{2}}\frac{1}{q^{2}+4u^{2}}+{\cal
O}\left((q^{2}+4u^{2})^{-2}\right)\;,$ (122) $\displaystyle D_{0}(q\ll
1,\frac{uq}{m})$ $\displaystyle\xrightarrow[\beta\to\infty]{}$
$\displaystyle-\frac{m}{\pi^{2}}\left[R_{0}(u)+R_{1}(u)\,q^{2}+R_{2}(u)\,q^{4}\right]+{\cal
O}\left(q^{6}\right)\;,$ (123)

where

$\displaystyle R_{0}(u)$ $\displaystyle\equiv$ $\displaystyle
p_{F}-u\tan^{-1}\frac{p_{F}}{u}\;,$ (124) $\displaystyle R_{1}(u)$
$\displaystyle=$ $\displaystyle-\frac{p_{F}^{3}}{12(p_{F}^{2}+u^{2})^{2}}\;,$
(125) $\displaystyle R_{2}(u)$ $\displaystyle=$
$\displaystyle-\frac{p_{F}^{3}(p_{F}^{2}-5u^{2})}{240(p_{F}^{2}+u^{2})^{4}}\;.$
(126)

With these asymptotic behaviors, we see from eq. (117) that $\Gamma_{1}[n]$ is
finite. Methods for extracting coefficients associated with the $e^{2}$-based
expansion for $\epsilon^{\prime}$ can be found in references Gell-Mann and
Brueckner, 1957 and Carr and Maradudin, 1964. The
$\epsilon^{\prime}=\lim_{\beta\to\infty}\Gamma_{1}[n]/(N_{e}\,\beta)-\epsilon_{x}$
part contains in general $r_{s}^{n\geq 0}\ln r_{s}$ and $r_{s}^{(n\geq 0)}$
terms when the energy is expressed in Rydbergs and expanded in power of
$r_{s}$ (or $e^{2}$).

We now turn our attention to $\Gamma_{2}[n]$, whose diagrammatic expression is
given in Figure 3. As mentioned earlier, within our formalism, $\Gamma_{i\geq
2}[n]$ correspond to diagrams containing only $\tilde{\cal D}_{0}$ propagator
and the KS electron propagator. Since the $\tilde{\cal D}_{0}$ propagator
retains a finite value even when its associated momentum approaches zero, each
of the diagrams corresponding to $\Gamma_{i\geq 2}[n]$ takes a finite value.
The absence of divergence no longer holds when one perform high density (
small $r_{s}$) expansion, as we will illustrate explicitly later. To make easy
comparison with existing HEG studies, one employs Dyson’s equation,

$\displaystyle\tilde{\cal D}_{0}$ $\displaystyle=$ $\displaystyle
U+U{\scriptstyle\circ}D_{0}{\scriptstyle\circ}\tilde{\cal D}_{0}$ (127)
\begin{picture}(20.0,30.0)(0.0,-3.0)\put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \end{picture} $\displaystyle=$
$\displaystyle\;\;\begin{picture}(20.0,30.0)(0.0,-3.0) \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){}
\end{picture}\;\;+\;\begin{picture}(60.0,30.0)(0.0,-3.0) \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \end{picture}\;\;$

to decompose the propagator $\tilde{\cal D}_{0}$ (of order $e^{2}$ and higher)
into the sum of a bare Coulomb propagator $U$ (of order $e^{2}$) and
$U{\scriptstyle\circ}D_{0}{\scriptstyle\circ}\tilde{\cal D}_{0}$ (of order
$e^{4}$ and higher). The reason that one should not expand further and write
$\tilde{\cal
D}_{0}=U+U{\scriptstyle\circ}D_{0}{\scriptstyle\circ}\,U+U{\scriptstyle\circ}D_{0}{\scriptstyle\circ}\,U{\scriptstyle\circ}D_{0}{\scriptstyle\circ}\tilde{\cal
D}_{0}$ is because the $U{\scriptstyle\circ}D_{0}{\scriptstyle\circ}\,U$ term
causes infrared divergence due to the momentum integral $\int
d{\mathbf{q}}/q^{4}$. In fact, this is exactly what causes (in the second
diagram of Figure 1) the infrared divergence, motivating the summation of the
ring diagrams.

To illustrate the main points, let us begin by examining the first three
diagrams of $\Gamma_{2}[n]$ and employ the decomposition rule mentioned above:

$\displaystyle\frac{1}{4}\;\;\;\;\;\;\;\;\;\;\;\;\;\begin{picture}(80.0,20.0)(0.0,-4.0)\put(-16.0,15.0){\makebox(0.0,0.0)[]{${\mathbf{p}}_{1}+{\mathbf{q}}$}}
\put(45.0,15.0){\makebox(0.0,0.0)[]{${\mathbf{p}}_{1}$}}
\put(-8.0,-15.0){\makebox(0.0,0.0)[]{$-{\mathbf{p}}_{2}$}}
\put(58.0,-15.0){\makebox(0.0,0.0)[]{$-{\mathbf{p}}_{2}-{\mathbf{q}}$}}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \end{picture}\;\;$ $\displaystyle=$
$\displaystyle\;\;\frac{1}{4}\;\begin{picture}(40.0,20.0)(0.0,-4.0)
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \raise-34.0pt\hbox to0.0pt{\kern
20.0pt\makebox(0.0,0.0)[cc]{$(a)$}\hss}
\ignorespaces\end{picture}\;\;+\frac{1}{2}\;\begin{picture}(50.0,25.0)(0.0,-4.0)
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\raise-34.0pt\hbox to0.0pt{\kern 25.0pt\makebox(0.0,0.0)[cc]{$(b)$}\hss}
\ignorespaces\end{picture}\;\;+\frac{1}{4}\;\begin{picture}(50.0,25.0)(0.0,-4.0)
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \raise-34.0pt\hbox to0.0pt{\kern
25.0pt\makebox(0.0,0.0)[cc]{$(c)$}\hss} \ignorespaces\end{picture}$ (128)
$\displaystyle\frac{1}{2}\;\;\;\;\;\;\;\;\;\;\;\;\;\begin{picture}(80.0,60.0)(0.0,-4.0)\put(-17.5,19.0){\makebox(0.0,0.0)[]{$-{\mathbf{p}}_{1}-{\mathbf{q}}$}}
\put(50.0,15.0){\makebox(0.0,0.0)[]{$-{\mathbf{p}}_{1}$}}
\put(-16.0,-19.0){\makebox(0.0,0.0)[]{$-{\mathbf{p}}_{2}-{\mathbf{q}}$}}
\put(60.0,-15.0){\makebox(0.0,0.0)[]{$-{\mathbf{p}}_{1}-{\mathbf{q}}$}}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \end{picture}\;\;$ $\displaystyle=$
$\displaystyle\;\;\frac{1}{2}\;\begin{picture}(40.0,20.0)(0.0,-4.0)
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \raise-34.0pt\hbox to0.0pt{\kern
20.0pt\makebox(0.0,0.0)[cc]{$(a)$}\hss}
\ignorespaces\end{picture}\;\;+\;\;\;\;\begin{picture}(50.0,25.0)(0.0,-4.0)
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\raise-34.0pt\hbox to0.0pt{\kern 25.0pt\makebox(0.0,0.0)[cc]{$(b)$}\hss}
\ignorespaces\end{picture}\;\;+\frac{1}{2}\;\begin{picture}(50.0,25.0)(0.0,-4.0)
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \raise-34.0pt\hbox to0.0pt{\kern
25.0pt\makebox(0.0,0.0)[cc]{$(c)$}\hss} \ignorespaces\end{picture}$ (129)
$\displaystyle-\frac{1}{2}\;\;\;\;\;\;\;\begin{picture}(108.0,60.0)(0.0,-4.0)\put(-5.5,15.0){\makebox(0.0,0.0)[]{${\mathbf{p}}_{2}$}}
\put(-6.5,8.0){\makebox(0.0,0.0)[]{$+$}}
\put(-6.5,2.0){\makebox(0.0,0.0)[]{${\mathbf{q}}^{\prime}$}}
\put(43.0,18.0){\makebox(0.0,0.0)[]{${\mathbf{p}}_{2}$}}
\put(43.0,-18.0){\makebox(0.0,0.0)[]{${\mathbf{p}}_{2}$}}
\put(59.0,18.0){\makebox(0.0,0.0)[]{${\mathbf{p}}_{1}$}}
\put(60.0,-18.0){\makebox(0.0,0.0)[]{${\mathbf{p}}_{1}$}}
\put(106.0,15.0){\makebox(0.0,0.0)[]{${\mathbf{p}}_{1}$}}
\put(107.0,8.0){\makebox(0.0,0.0)[]{$+$}}
\put(107.0,1.0){\makebox(0.0,0.0)[]{${\mathbf{q}}$}} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \end{picture}\;\;$ $\displaystyle=$
$\displaystyle\;\;-\frac{1}{2}\;\begin{picture}(100.0,60.0)(0.0,-4.0)
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \raise-20.0pt\hbox
to0.0pt{\kern 50.0pt\makebox(0.0,0.0)[cc]{$(a)$}\hss}
\ignorespaces\end{picture}\;\;-\;\;\begin{picture}(100.0,60.0)(0.0,-4.0)
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \raise-20.0pt\hbox to0.0pt{\kern
50.0pt\makebox(0.0,0.0)[cc]{$(b)$}\hss} \ignorespaces\end{picture}\;$ (130)
$\displaystyle\;\;-\frac{1}{2}\;\;\begin{picture}(100.0,50.0)(0.0,-4.0)
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \raise-20.0pt\hbox
to0.0pt{\kern 50.0pt\makebox(0.0,0.0)[cc]{$(c)$}\hss}
\ignorespaces\end{picture}$

Let us start with the diagrams to the left of the equal signs in (128-130). We
will for now keep the bosonic propagators general, that is, allowing them to
carry frequencies. After that, we will discuss the (a) diagrams in (128-130),
followed by the (b) diagrams and then the (c) diagrams.

In (128), let the vertical boson propagator, denoted by $B_{1}$, carry
momentum ${\mathbf{q}}$ (upward) and frequency $\nu$. Let the horizontal boson
propagator, denoted by $B_{2}$, carry momentum
${\mathbf{q}}^{\prime}\equiv-({\mathbf{p}}_{1}+{\mathbf{p}}_{2}+{\mathbf{q}})$
(leftward) and frequency $\nu^{\prime}$. Using techniques described in section
III.5 for frequency summation, one obtains the following generic result

$\displaystyle\frac{\beta
V}{4}\frac{2}{\beta^{2}}\sum_{\nu,\nu^{\prime}}\int\frac{d{\mathbf{p}}_{1}}{(2\pi)^{3}}\frac{d{\mathbf{p}}_{2}}{(2\pi)^{3}}\frac{d{\mathbf{q}}}{(2\pi)^{3}}\frac{B_{1}({\mathbf{q}},\nu)}{-i\nu+\varepsilon_{{\mathbf{p}}_{1}+{\mathbf{q}}}-\varepsilon_{{\mathbf{p}}_{1}}}\frac{B_{2}({\mathbf{q}}^{\prime},\nu^{\prime})}{i\nu+\varepsilon_{{\mathbf{p}}_{2}+{\mathbf{q}}}-\varepsilon_{{\mathbf{p}}_{2}}}\left\\{\frac{n_{{\mathbf{p}}_{2}}-n_{{\mathbf{p}}_{1}}}{\varepsilon_{{\mathbf{p}}_{2}}-\varepsilon_{{\mathbf{p}}_{1}}-i(\nu+\nu^{\prime})}\right.$
$\displaystyle\left.\hskip
15.0pt+\frac{n_{{\mathbf{p}}_{1}}-n_{{\mathbf{p}}_{1}+{\mathbf{q}}^{\prime}}}{\varepsilon_{{\mathbf{p}}_{1}+{\mathbf{q}}^{\prime}}-\varepsilon_{{\mathbf{p}}_{1}}-i\nu^{\prime}}+\frac{n_{{\mathbf{p}}_{2}}-n_{{\mathbf{p}}_{2}+{\mathbf{q}}^{\prime}}}{\varepsilon_{{\mathbf{p}}_{2}+{\mathbf{q}}^{\prime}}-\varepsilon_{{\mathbf{p}}_{2}}+i\nu^{\prime}}+\frac{n_{{\mathbf{p}}_{2}+{\mathbf{q}}}-n_{{\mathbf{p}}_{1}+{\mathbf{q}}}}{\varepsilon_{{\mathbf{p}}_{2}+{\mathbf{q}}}-\varepsilon_{{\mathbf{p}}_{1}+{\mathbf{q}}}+i(\nu-\nu^{\prime})}\right\\}\;,$
(131)

or equivalently,

$\displaystyle\frac{\beta
V}{4}\frac{2}{\beta^{2}}\sum_{\nu,\nu^{\prime}}\int\frac{d{\mathbf{p}}_{1}}{(2\pi)^{3}}\frac{d{\mathbf{p}}_{2}}{(2\pi)^{3}}\frac{d{\mathbf{q}}}{(2\pi)^{3}}\frac{B_{1}({\mathbf{q}},\nu)}{-i\nu^{\prime}+\varepsilon_{{\mathbf{p}}_{1}+{\mathbf{q}}^{\prime}}-\varepsilon_{{\mathbf{p}}_{1}}}\frac{B_{2}({\mathbf{q}}^{\prime},\nu^{\prime})}{i\nu^{\prime}+\varepsilon_{{\mathbf{p}}_{2}+{\mathbf{q}}^{\prime}}-\varepsilon_{{\mathbf{p}}_{2}}}\left\\{\frac{n_{{\mathbf{p}}_{2}}-n_{{\mathbf{p}}_{1}}}{\varepsilon_{{\mathbf{p}}_{2}}-\varepsilon_{{\mathbf{p}}_{1}}-i(\nu+\nu^{\prime})}\right.$
$\displaystyle\left.\hskip
15.0pt+\frac{n_{{\mathbf{p}}_{1}}-n_{{\mathbf{p}}_{1}+{\mathbf{q}}}}{\varepsilon_{{\mathbf{p}}_{1}+{\mathbf{q}}}-\varepsilon_{{\mathbf{p}}_{1}}-i\nu}+\frac{n_{{\mathbf{p}}_{2}}-n_{{\mathbf{p}}_{2}+{\mathbf{q}}}}{\varepsilon_{{\mathbf{p}}_{2}+{\mathbf{q}}}-\varepsilon_{{\mathbf{p}}_{2}}+i\nu}+\frac{n_{{\mathbf{p}}_{2}+{\mathbf{q}}}-n_{{\mathbf{p}}_{1}+{\mathbf{q}}}}{\varepsilon_{{\mathbf{p}}_{2}+{\mathbf{q}}}-\varepsilon_{{\mathbf{p}}_{1}+{\mathbf{q}}}+i(\nu-\nu^{\prime})}\right\\}\;.$
(132)

Note that the factor $2$ associated with $\frac{2}{\beta^{2}}$ comes from the
two possible spin states for an electron.

In (129), let the boson propagator on top, denoted by $B_{1}$, carry momentum
${\mathbf{q}}$ (towards lower-right) and frequency $\nu$. Let the boson
propagator at the bottom, denoted by $B_{2}$, carry momentum
${\mathbf{q}}^{\prime}\equiv({\mathbf{p}}_{1}-{\mathbf{p}}_{2})$ (towards
upper-left) and frequency $\nu^{\prime}$. Using techniques described in
section III.5 for frequency summation, one obtains the following generic
result

$\displaystyle\frac{\beta
V}{2}\frac{2}{\beta^{2}}\sum_{\nu,\nu^{\prime}}\int\frac{d{\mathbf{p}}_{1}}{(2\pi)^{3}}\frac{d{\mathbf{p}}_{2}}{(2\pi)^{3}}\frac{d{\mathbf{q}}}{(2\pi)^{3}}\frac{B_{1}({\mathbf{q}},\nu)}{-i\nu+\varepsilon_{{\mathbf{p}}_{1}+{\mathbf{q}}}-\varepsilon_{{\mathbf{p}}_{1}}}\frac{B_{2}({\mathbf{q}}^{\prime},\nu^{\prime})}{-i\nu^{\prime}+\varepsilon_{{\mathbf{p}}_{1}+{\mathbf{q}}}-\varepsilon_{{\mathbf{p}}_{2}+{\mathbf{q}}}}\left\\{-\beta
n_{{\mathbf{p}}_{1}+{\mathbf{q}}}(1-n_{{\mathbf{p}}_{1}+{\mathbf{q}}})\right.$
$\displaystyle\left.\hskip
15.0pt-\frac{n_{{\mathbf{p}}_{1}+{\mathbf{q}}}-n_{{\mathbf{p}}_{2}+{\mathbf{q}}}}{-i\nu^{\prime}+\varepsilon_{{\mathbf{p}}_{1}+{\mathbf{q}}}-\varepsilon_{{\mathbf{p}}_{2}+{\mathbf{q}}}}+\frac{n_{{\mathbf{p}}_{1}}-n_{{\mathbf{p}}_{1}+{\mathbf{q}}}}{-i\nu+\varepsilon_{{\mathbf{p}}_{1}+{\mathbf{q}}}-\varepsilon_{{\mathbf{p}}_{1}}}+\frac{n_{{\mathbf{p}}_{1}}-n_{{\mathbf{p}}_{2}+{\mathbf{q}}}}{i(\nu-\nu^{\prime})+\varepsilon_{{\mathbf{p}}_{1}}-\varepsilon_{{\mathbf{p}}_{2}+{\mathbf{q}}}}\right\\}\;.$
(133)

In (130), let the boson propagator at the right hand side, denoted by $B_{1}$,
carry momentum ${\mathbf{q}}$ (downwards) and frequency $\nu$. Let the boson
propagator at the left hand side, denoted by $B_{2}$, carry momentum
${\mathbf{q}}^{\prime}$ (also downwards) and frequency $\nu^{\prime}$. Using
techniques described in section III.5 for frequency summation, one obtains the
following generic result

$\displaystyle-\frac{\beta
V}{2}\frac{2^{2}}{\beta^{2}}\frac{-1}{2\beta\int\frac{d{\mathbf{p}}}{(2\pi)^{3}}n_{{\mathbf{p}}}(1-n_{{\mathbf{p}}})}\sum_{\nu,\nu^{\prime}}\int\frac{d{\mathbf{p}}_{1}}{(2\pi)^{3}}\frac{d{\mathbf{p}}_{2}}{(2\pi)^{3}}\frac{d{\mathbf{q}}}{(2\pi)^{3}}\frac{d{\mathbf{q}}^{\prime}}{(2\pi)^{3}}\frac{B_{1}({\mathbf{q}},\nu)}{-i\nu+\varepsilon_{{\mathbf{p}}_{1}+{\mathbf{q}}}-\varepsilon_{{\mathbf{p}}_{1}}}\frac{B_{2}({\mathbf{q}}^{\prime},\nu^{\prime})}{-i\nu^{\prime}+\varepsilon_{{\mathbf{p}}_{2}+{\mathbf{q}}^{\prime}}-\varepsilon_{{\mathbf{p}}_{2}}}$
$\displaystyle\hskip
15.0pt\left\\{\beta^{2}n_{{\mathbf{p}}_{1}}(1-n_{{\mathbf{p}}_{1}})n_{{\mathbf{p}}_{2}}(1-n_{{\mathbf{p}}_{2}})+\beta
n_{{\mathbf{p}}_{1}}(1-n_{{\mathbf{p}}_{1}})\frac{n_{{\mathbf{p}}_{2}+{\mathbf{q}}^{\prime}}-n_{{\mathbf{p}}_{2}}}{-i\nu^{\prime}+\varepsilon_{{\mathbf{p}}_{2}+{\mathbf{q}}^{\prime}}-\varepsilon_{{\mathbf{p}}_{2}}}\right.$
$\displaystyle\hskip 20.0pt\left.+\beta
n_{{\mathbf{p}}_{2}}(1-n_{{\mathbf{p}}_{2}})\frac{n_{{\mathbf{p}}_{1}+{\mathbf{q}}}-n_{{\mathbf{p}}_{1}}}{-i\nu+\varepsilon_{{\mathbf{p}}_{1}+{\mathbf{q}}}-\varepsilon_{{\mathbf{p}}_{1}}}+\frac{n_{{\mathbf{p}}_{2}+{\mathbf{q}}^{\prime}}-n_{{\mathbf{p}}_{2}}}{-i\nu^{\prime}+\varepsilon_{{\mathbf{p}}_{2}+{\mathbf{q}}^{\prime}}-\varepsilon_{{\mathbf{p}}_{2}}}\frac{n_{{\mathbf{p}}_{1}+{\mathbf{q}}}-n_{{\mathbf{p}}_{1}}}{-i\nu+\varepsilon_{{\mathbf{p}}_{1}+{\mathbf{q}}}-\varepsilon_{{\mathbf{p}}_{1}}}\right\\}\;,$
(134)

where the factor
$-1/\left\\{2\beta\int\frac{d{\mathbf{p}}}{(2\pi)^{3}}n_{\mathbf{p}}(1-n_{\mathbf{p}})\right\\}=D_{0}^{-1}({\mathbf{q}}=0,\nu=0)$
arises from the fact that in a translationally invariant system, momentum
conservation at each vertex demands that the inverse density correlator must
carry zero momentum and frequency.

For the diagram (a) on the right hand side of eq. (128),
$B_{1}({\mathbf{q}},\nu)=4\pi e^{2}/{\mathbf{q}}^{2}$ and
$B_{2}({\mathbf{q}}^{\prime},\nu^{\prime})=4\pi
e^{2}/({\mathbf{p}}_{1}+{\mathbf{p}}_{2}+{\mathbf{q}})^{2}$ are both frequency
independent. Thus, one may sum over both $\nu$ and $\nu^{\prime}$. And upon
doing so, we obtain

$\frac{1}{4\beta}\;\begin{picture}(40.0,20.0)(0.0,-4.0) \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){}
\end{picture}\;=\;{V}\int\frac{d{\mathbf{p}}_{1}}{(2\pi)^{3}}\frac{d{\mathbf{p}}_{2}}{(2\pi)^{3}}\frac{d{\mathbf{q}}}{(2\pi)^{3}}\frac{(4\pi
e^{2})^{2}}{{\mathbf{q}}^{2}({\mathbf{p}}_{1}+{\mathbf{p}}_{2}+{\mathbf{q}})^{2}}\frac{n_{{\mathbf{p}}_{1}}n_{{\mathbf{p}}_{2}}(1-n_{{\mathbf{p}}_{1}+{\mathbf{q}}})(1-n_{{\mathbf{p}}_{2}+{\mathbf{q}}})}{\varepsilon_{{\mathbf{p}}_{1}+{\mathbf{q}}}+\varepsilon_{{\mathbf{p}}_{2}+{\mathbf{q}}}-\varepsilon_{{\mathbf{p}}_{1}}-\varepsilon_{{\mathbf{p}}_{2}}}\;.\vspace*{8pt}$

If we divide the quantity above by $N_{e}$ and then take the zero temperature
limit, it gives rise to (using $p_{F}$ as the unit for momentum)

$\displaystyle\frac{mp_{F}^{3}V}{N_{e}(2\pi)^{9}}(4\pi
e^{2})^{2}\int_{|{\mathbf{p}}_{i}|<1,{\mathbf{p}}_{i}+{\mathbf{q}}|>1}\;d{\mathbf{q}}\;d{\mathbf{p}}_{1}\;d{\mathbf{p}}_{2}\frac{1}{{\mathbf{q}}^{2}({\mathbf{p}}_{1}+{\mathbf{p}}_{2}+{\mathbf{q}})^{2}}\frac{1}{{\mathbf{q}}^{2}+{\mathbf{q}}\cdot({\mathbf{p}}_{1}+{\mathbf{p}}_{2})}$
$\displaystyle=\left[\frac{3}{16\pi^{5}}\int_{|{\mathbf{p}}_{i}|<1,{\mathbf{p}}_{i}+{\mathbf{q}}|>1}\;d{\mathbf{q}}\;d{\mathbf{p}}_{1}\;d{\mathbf{p}}_{2}\frac{1}{{\mathbf{q}}^{2}({\mathbf{p}}_{1}+{\mathbf{p}}_{2}+{\mathbf{q}})^{2}}\frac{1}{{\mathbf{q}}^{2}+{\mathbf{q}}\cdot({\mathbf{p}}_{1}+{\mathbf{p}}_{2})}\right]\;\;{\rm
Rydberg},$

which is exactly the $\epsilon_{b}^{(2)}$ term of Gell-Mann and
Brueckner.Gell-Mann and Brueckner (1957)

Diagram (a) on the right hand side of eq. (129) is one of the anomalous
diagramsKohn and Luttinger (1960); Luttinger and Ward (1960) that give rise to
finite contribution as the $T\to 0$ limit is taken within finite temperature
formalism but are absent within zero temperature formalism. For this diagram,
$B_{1}({\mathbf{q}},\nu)=4\pi e^{2}/{\mathbf{q}}^{2}$ and
$B_{2}({\mathbf{q}}^{\prime},\nu^{\prime})=4\pi
e^{2}/({\mathbf{p}}_{1}-{\mathbf{p}}_{2})^{2}$ are both frequency independent.
Thus, one may evaluate its zero temperature contribution by summing over both
$\nu$ and $\nu^{\prime}$. Upon doing so, one obtains from diagram (a) of (129)
the following anomalous contribution

$\displaystyle\frac{1}{2\beta}\;\begin{picture}(40.0,20.0)(0.0,-4.0)
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \end{picture}\;\;$ $\displaystyle=$
$\displaystyle-(V)\int\frac{d{\mathbf{p}}_{1}\,d{\mathbf{p}}_{2}\,d{\mathbf{q}}}{(2\pi)^{9}}\frac{4\pi
e^{2}}{{\mathbf{q}}^{2}}\frac{4\pi
e^{2}}{({\mathbf{p}}_{1}-{\mathbf{p}}_{2})^{2}}\,\beta
n_{{\mathbf{p}}_{1}+{\mathbf{q}}}(1-n_{{\mathbf{p}}_{1}+{\mathbf{q}}})n_{{\mathbf{p}}_{2}+{\mathbf{q}}}n_{{\mathbf{p}}_{1}}$
(135) $\displaystyle=$
$\displaystyle-(V)\int\frac{d{\mathbf{p}}_{1}\,d{\mathbf{q}}\,d{\mathbf{q}}^{\prime}}{(2\pi)^{9}}\frac{4\pi
e^{2}}{{\mathbf{q}}^{2}}\frac{4\pi e^{2}}{({\mathbf{q}}^{\prime})^{2}}\,\beta
n_{{\mathbf{p}}_{1}}(1-n_{{\mathbf{p}}_{1}})n_{{\mathbf{p}}_{1}+{\mathbf{q}}^{\prime}}n_{{\mathbf{p}}_{1}+{\mathbf{q}}}\;,$

where the last expression is obtained via the following change of variables:
${\mathbf{p}}_{1(2)}\to-{\mathbf{p}}_{1(2)}-{\mathbf{q}}$, eliminating
${\mathbf{p}}_{2}$ by
${\mathbf{q}}^{\prime}\equiv{\mathbf{p}}_{1}-{\mathbf{p}}_{2}$, and
${\mathbf{q}}^{\prime}\to-{\mathbf{q}}^{\prime}$. This expression can be
further simplified via the following definition

$f({\mathbf{p}})\equiv\int\frac{d{\mathbf{q}}}{(2\pi)^{3}}\frac{1}{{\mathbf{q}}^{2}}n_{{\mathbf{p}}+{\mathbf{q}}}\;,$

leading to

$\displaystyle\frac{1}{2\beta}\;\begin{picture}(40.0,20.0)(0.0,-4.0)
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \end{picture}\;\;$ $\displaystyle=$
$\displaystyle-V(4\pi e^{2})^{2}\int\frac{d{\mathbf{p}}_{1}}{(2\pi)^{3}}\beta
n_{{\mathbf{p}}_{1}}(1-n_{{\mathbf{p}}_{1}})f^{2}({\mathbf{p}}_{1})$ (136)
$\displaystyle=$ $\displaystyle-V(4\pi
e^{2})^{2}\int\frac{d{\mathbf{p}}_{1}}{(2\pi)^{3}}\frac{\partial
n_{{\mathbf{p}}_{1}}}{\partial\mu}f^{2}({\mathbf{p}}_{1})\;.\vspace*{18pt}$

As mentioned earlier, each two-particle reducible diagram within our formalism
is accompanied by a corresponding diagram that will eliminate anomalous
contribution when applicable. Diagrams in eq. (130), appearing only under the
effective action formalism, are such diagrams. They are neither present in the
zero temperature formalism of Goldstone Goldstone and Brueckner Brueckner and
Levinson (1955) nor in the finite temperature formalism of Luttinger and
Ward.Luttinger and Ward (1960) Diagram (a) on the right hand side of (130)
corresponds to the case $B_{1}({\mathbf{q}},\nu)=4\pi e^{2}/{\mathbf{q}}^{2}$
and $B_{2}({\mathbf{q}}^{\prime},\nu^{\prime})=4\pi
e^{2}/({\mathbf{q}}^{\prime})^{2}$, both being frequency independent. Thus,
upon summing over both $\nu$ and $\nu^{\prime}$, one obtains

$\displaystyle-\frac{1}{2\beta}\;\begin{picture}(100.0,40.0)(0.0,-4.0)
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \end{picture}$
$\displaystyle=$ $\displaystyle\frac{2(-1)^{2}V(4\pi
e^{2})^{2}}{2\beta\int\frac{d{\mathbf{p}}}{(2\pi)^{3}}n_{\mathbf{p}}(1-n_{\mathbf{p}})}\int\frac{d{\mathbf{p}}_{1}\,d{\mathbf{p}}_{2}\,d{\mathbf{q}}\,d{\mathbf{q}}^{\prime}}{(2\pi)^{12}}\frac{4\pi
e^{2}}{q^{2}}\frac{4\pi e^{2}}{{q^{\prime}}^{2}}\times$ (137)
$\displaystyle\times\beta^{2}n_{{\mathbf{p}}_{1}}(1-n_{{\mathbf{p}}_{1}})n_{{\mathbf{p}}_{2}}(1-n_{{\mathbf{p}}_{2}})n_{{\mathbf{p}}_{1}+{\mathbf{q}}}n_{{\mathbf{p}}_{2}+{\mathbf{q}}^{\prime}}$
$\displaystyle=$ $\displaystyle\frac{V(4\pi
e^{2})^{2}}{\int\frac{d{\mathbf{p}}}{(2\pi)^{3}}\frac{\partial
n_{\mathbf{p}}}{\partial\mu}}\left[\int\frac{d{\mathbf{p}}}{(2\pi)^{3}}\frac{\partial
n_{\mathbf{p}}}{\partial\mu}f({\mathbf{p}})\right]^{2}\;.$ (138)

When combined with eq. (136), one obtains

$\frac{1}{2\beta}\;\begin{picture}(40.0,20.0)(0.0,-4.0) \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){}
\end{picture}\;\;-\frac{1}{2\beta}\;\begin{picture}(100.0,30.0)(0.0,-4.0)
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \end{picture}=-V(4\pi
e^{2})^{2}\left(\int\frac{d{\mathbf{p}}}{(2\pi)^{3}}\frac{\partial
n_{\mathbf{p}}}{\partial\mu}\right)\overline{\left(f-\overline{f}\right)^{2}}\;,\vspace*{10pt}$
(139)

where the overline symbol is defined as
$\overline{f}\equiv\int\frac{d{\mathbf{p}}}{(2\pi)^{3}}\frac{\partial
n_{\mathbf{p}}}{\partial\mu}f({\mathbf{p}})$. Apparently, the expression (139)
is in general negative unless $f=\overline{f}$. We will show that this happens
only at zero temperature and only if the system has spherical symmetry.

In order to have $f({\mathbf{p}})=\overline{f}$, ${\mathbf{p}}$ can only have
support at a constant $f({\mathbf{p}})$ surface. This is achieved when $T\to
0$, where $n_{\mathbf{p}}\to\theta(\mu-\varepsilon_{\mathbf{p}})$ and
$\partial n_{\mathbf{p}}/\partial\mu=\delta(\mu-\varepsilon_{\mathbf{p}})$
forcing ${\mathbf{p}}$ to lie on a constant $\varepsilon_{\mathbf{p}}$
surface. Furthermore, $f({\mathbf{p}})=\overline{f}$ demands that
$f({\mathbf{p}})$ depends only on the magnitude of ${\mathbf{p}}$, i.e.,
$f({\mathbf{p}})=f(p)$. This can only be achieved if $n_{\mathbf{p}}$ depends
only on $|{\mathbf{p}}|=p$. For the HEG system,
$\varepsilon_{\mathbf{p}}={\mathbf{p}}^{2}/2m$, thus $n_{\mathbf{p}}=n(p)$ and
$\mu=p_{F}^{2}/2m$. We therefore have $\partial
n_{\mathbf{p}}/\partial\mu=m\delta(p-p_{F})/p_{F}$, fixing the length of
${\mathbf{p}}$. Now
$n_{{\mathbf{p}}+{\mathbf{q}}}=\theta(\mu-({\mathbf{p}}+{\mathbf{q}})^{2}/2m)|_{|{\mathbf{p}}|=p_{F}}=\theta(-\cos\vartheta-\frac{q}{2p_{F}})$,
with $\vartheta$ being the angle between ${\mathbf{q}}$ and ${\mathbf{p}}$.
Thus, in the integral defining $f({\mathbf{p}})$, although $\hat{\mathbf{p}}$
defines the $\hat{z}$ direction of vector ${\mathbf{q}}$ the integral is
independent of the choice of $\hat{\mathbf{p}}$. Therefore, the anomalous
contribution (136) may be written as (when $T\to 0$)

$-(V)\left[\int\frac{d{\mathbf{p}}_{1}}{(2\pi)^{3}}\frac{m}{p_{F}}\delta(p_{1}-p_{F})\right]\left[\int_{-1}^{1}dx\int_{0}^{\infty}\frac{q^{2}\,dq}{(2\pi)^{2}}\frac{4\pi
e^{2}}{q^{2}}\;\theta\big{(}-x-\frac{q}{2p_{F}}\big{)}\right]^{2}\;,$ (140)

while the corresponding subtraction term (138) may be written as

$\displaystyle\frac{V}{\int\frac{d{\mathbf{p}}_{1}}{(2\pi)^{3}}\frac{m}{p_{F}}\delta(p_{1}-p_{F})}\left[\int\frac{d{\mathbf{p}}_{1}}{(2\pi)^{3}}\frac{m}{p_{F}}\delta(p_{1}-p_{F})\right]^{2}\left[\int_{-1}^{1}dx\int_{0}^{\infty}\frac{q^{2}\,dq}{(2\pi)^{2}}\frac{4\pi
e^{2}}{q^{2}}\;\theta\big{(}-x-\frac{q}{2p_{F}}\big{)}\right]^{2}\;$
$\displaystyle=V\left[\int\frac{d{\mathbf{p}}_{1}}{(2\pi)^{3}}\frac{m}{p_{F}}\delta(p_{1}-p_{F})\right]\left[\int_{-1}^{1}dx\int_{0}^{\infty}\frac{q^{2}\,dq}{(2\pi)^{2}}\frac{4\pi
e^{2}}{q^{2}}\;\theta\big{(}-x-\frac{q}{2p_{F}}\big{)}\right]^{2}\;,$ (141)

cancelling exactly the anomalous contribution (140) in the HEG case.

Within the framework of Luttinger et al.,Kohn and Luttinger (1960); Luttinger
and Ward (1960) the anomalous contribution is cancelled by the chemical
potential shift, the difference between $\mu(T\to 0)$, under finite
temperature formalism, and $\mu(T=0)=p_{F}^{2}/2m$, under zero temperature
formalism. Within the current finite temperature framework, however, the
$\mu(T\to 0)$ is identical to $\mu(T=0)$ and the cancellation of anomalous
contribution is explicit. As we will show later, however, diagram (b) of (129)
can’t be cancelled by diagram (b) of (130). This is because diagram (b) of
(129) is not an anomalous diagram.

Before moving onto (b) diagrams in (128-130), let us remark that diagrams
within our formalism, i.e., diagrams on the left hand side of (128-130) with
$B_{1(2)}\to\tilde{\cal D}_{0}$ yield only finite contribution. The
decomposition made in (128-130), however, may introduce divergence as we will
illustrate using the (b) diagrams that correspond to the main results of
DuBois.DuBois (1959)

For the (b) diagram of (128), $B_{2}({\mathbf{q}}^{\prime},\nu^{\prime})=4\pi
e^{2}/({\mathbf{p}}_{1}+{\mathbf{p}}_{2}+{\mathbf{q}})^{2}$ while
$B_{1}({\mathbf{q}},\nu)=(4\pi
e^{2})^{2}D_{0}({\mathbf{q}},\nu)/\left\\{{\mathbf{q}}^{2}[{\mathbf{q}}^{2}-4\pi
e^{2}D_{0}({\mathbf{q}},\nu)]\right\\}$. Therefore, one may sum over
$\nu^{\prime}$ to simplify the expression. Using again the methods described
in section III.5 for frequency summation, one obtains

$\displaystyle\frac{1}{2\beta}\;\;\begin{picture}(50.0,25.0)(0.0,-4.0)
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\end{picture}\;$ $\displaystyle=$ $\displaystyle\;2\frac{V}{4}\frac{(4\pi
e^{2})^{3}}{(2\pi)^{9}}\frac{2}{\beta}\sum_{\nu}\int
d{\mathbf{p}}_{1}d{\mathbf{p}}_{2}d{\mathbf{q}}\frac{D_{0}({\mathbf{q}},\nu)}{{\mathbf{q}}^{2}({\mathbf{q}}^{2}-4\pi
e^{2}D_{0}({\mathbf{q}},\nu))}\frac{1}{({\mathbf{p}}_{1}+{\mathbf{p}}_{2}+{\mathbf{q}})^{2}}\times$
(142) $\displaystyle\hskip
60.0pt\times\frac{n_{{\mathbf{p}}_{1}+{\mathbf{q}}}-n_{{\mathbf{p}}_{1}}}{-i\nu+\varepsilon_{{\mathbf{p}}_{1}+{\mathbf{q}}}-\varepsilon_{{\mathbf{p}}_{1}}}\,\frac{n_{{\mathbf{p}}_{2}+{\mathbf{q}}}-n_{{\mathbf{p}}_{2}}}{i\nu+\varepsilon_{{\mathbf{p}}_{2}+{\mathbf{q}}}-\varepsilon_{{\mathbf{p}}_{2}}}$
$\displaystyle=$ $\displaystyle\;\frac{V}{\beta}\frac{(4\pi
e^{2})^{3}}{(2\pi)^{9}}\sum_{\nu}\int
d{\mathbf{p}}_{1}d{\mathbf{p}}_{2}d{\mathbf{q}}\frac{D_{0}({\mathbf{q}},\nu)}{{\mathbf{q}}^{2}({\mathbf{q}}^{2}-4\pi
e^{2}D_{0}({\mathbf{q}},\nu))}\frac{1}{({\mathbf{p}}_{1}-{\mathbf{p}}_{2})^{2}}\times$
$\displaystyle\hskip
60.0pt\times\frac{n_{{\mathbf{p}}_{1}+{\mathbf{q}}}-n_{{\mathbf{p}}_{1}}}{-i\nu+\varepsilon_{{\mathbf{p}}_{1}+{\mathbf{q}}}-\varepsilon_{{\mathbf{p}}_{1}}}\,\frac{n_{{\mathbf{p}}_{2}+{\mathbf{q}}}-n_{{\mathbf{p}}_{2}}}{-i\nu+\varepsilon_{{\mathbf{p}}_{2}+{\mathbf{q}}}-\varepsilon_{{\mathbf{p}}_{2}}}\;,$

where the last expression is obtained by changing variables:
${\mathbf{p}}_{2}+{\mathbf{q}}\to-{\mathbf{p}}_{2}^{\prime}$ followed by
${\mathbf{p}}_{2}^{\prime}\to{\mathbf{p}}_{2}$. For the high density
expansion, where $e^{2}$ is treated as a small parameter, the major
contribution in (142) comes from the region ${\mathbf{q}}\to 0$. One thus
writes

$n_{{\mathbf{p}}+{\mathbf{q}}}-n_{\mathbf{p}}\xrightarrow[|{\mathbf{q}}|\to
0]{}(\varepsilon_{{\mathbf{p}}+q}-\varepsilon_{{\mathbf{p}}})\frac{\partial
n_{{\mathbf{p}}}}{\partial\varepsilon_{{\mathbf{p}}}}=-\beta
n_{{\mathbf{p}}}(1-n_{{\mathbf{p}}})(\varepsilon_{{\mathbf{p}}+{\mathbf{q}}}-\varepsilon_{{\mathbf{p}}})\xrightarrow[T\to
0]{}-q\cos\vartheta\,\delta(p-p_{F})\;,$

where $\vartheta$ is the angle between ${\mathbf{p}}$ and ${\mathbf{q}}$. As
$T\to 0$,
$\frac{1}{\beta}\sum_{n}F(\nu_{n})\to\int_{-\infty}^{\infty}\frac{d\nu}{2\pi}F(\nu)$
if $F(\nu)$ does not have pole strength greater than one. Making a change of
variable $\nu\equiv u|{\mathbf{q}}|/m$ and treating $q$ as a small quantity,
the major contribution of the (b) diagram of (128) is given by

$\displaystyle\frac{1}{2\beta}\;\;\begin{picture}(50.0,25.0)(0.0,-4.0)
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \end{picture}$
$\displaystyle\approx$ $\displaystyle\frac{V}{m}\frac{(4\pi
e^{2})^{3}}{(2\pi)^{9}}\int_{-\infty}^{\infty}\frac{du}{2\pi}\int
d{\mathbf{p}}_{1}\,d{\mathbf{p}}_{2}\frac{m\delta(p_{1}-p_{F})\cos\vartheta_{1}}{-iu+p_{F}\cos\vartheta_{1}}\frac{m\delta(p_{2}-p_{F})\cos\vartheta_{2}}{iu+p_{F}\cos\vartheta_{2}}\times$
(143) $\displaystyle\hskip
10.0pt\frac{1}{({\mathbf{p}}_{1}-{\mathbf{p}}_{2})^{2}}\int_{0}^{q_{c}}4\pi
qdq\frac{D_{0}(q,uq/m)}{q^{2}-4\pi e^{2}D_{0}(q,uq/m)}$ $\displaystyle\approx$
$\displaystyle Vm\,\frac{(4\pi
e^{2})^{3}}{(2\pi)^{9}}\int_{-\infty}^{\infty}\frac{du}{2\pi}\int_{-1}^{1}\frac{\cos\vartheta_{1}\,d\cos\vartheta_{1}}{-iu+p_{F}\cos\vartheta_{1}}\,\frac{\cos\vartheta_{2}\,d\cos\vartheta_{2}}{-iu+p_{F}\cos\vartheta_{2}}\times$
$\displaystyle\hskip
10.0pt\frac{(2\pi)^{2}p_{F}^{2}/2}{|\cos\vartheta_{1}-\cos\vartheta_{2}|}\int
4\pi qdq\frac{D_{0}(q,uq/m)}{q^{2}-4\pi e^{2}D_{0}(q,uq/m)}\;.$

Apparently, the $1/|\cos\vartheta_{1}-\cos\vartheta_{2}|$ factor in the
integrand causes an undesirable divergence in the $\vartheta$ integral that is
absent if one does not decompose the diagrams using (127). A finite result
emerges only when one combines (143) with the (b) diagram of (129), which we
soon turn to.

Before elaborating on the evaluation of the (b) diagram of (129), we wish to
point out that this is the type of diagram (two-particle reducible ones) that
was missed in Hedin’s approximation Hedin (1965) but should have been included
for the exactness of the theory. For this diagram,
$B_{2}({\mathbf{q}}^{\prime},\nu^{\prime})=4\pi
e^{2}/({\mathbf{p}}_{1}-{\mathbf{p}}_{2})^{2}$ while
$B_{1}({\mathbf{q}},\nu)=(4\pi
e^{2})^{2}D_{0}({\mathbf{q}},\nu)/\left\\{{\mathbf{q}}^{2}[{\mathbf{q}}^{2}-4\pi
e^{2}D_{0}({\mathbf{q}},\nu)]\right\\}$. Summing over $\nu^{\prime}$ via
methods described in section III.5, one obtains

$\displaystyle\frac{1}{\beta}\;\begin{picture}(50.0,25.0)(0.0,-4.0)
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \end{picture}$
$\displaystyle=$ $\displaystyle V\frac{(4\pi
e^{2})^{3}}{(2\pi)^{9}}\frac{2}{\beta}\sum_{\nu}\int
d{\mathbf{p}}_{1}d{\mathbf{p}}_{2}\,d{\mathbf{q}}\frac{1}{({\mathbf{p}}_{1}-{\mathbf{p}}_{2})^{2}}\frac{D_{0}(q,\nu)}{q^{2}(q^{2}-4\pi
e^{2}D_{0}(q,\nu))}\times$ $\displaystyle\hskip
10.0pt\left[-\frac{n_{{\mathbf{p}}_{2}+{\mathbf{q}}}(n_{{\mathbf{p}}_{1}+{\mathbf{q}}}-n_{{\mathbf{p}}_{1}})}{(-i\nu+\varepsilon_{{\mathbf{p}}_{1}+{\mathbf{q}}}-\varepsilon_{{\mathbf{p}}_{1}})^{2}}-\frac{\beta
n_{{\mathbf{p}}_{1}}(1-n_{{\mathbf{p}}_{1}})n_{{\mathbf{p}}_{2}}}{-i\nu+\varepsilon_{{\mathbf{p}}_{1}}-\varepsilon_{{\mathbf{p}}_{1}+{\mathbf{q}}}}\right]$
$\displaystyle=$ $\displaystyle\frac{V}{\beta}\frac{(4\pi
e^{2})^{3}}{(2\pi)^{9}}\sum_{\nu}\int
d{\mathbf{p}}_{1}d{\mathbf{p}}_{2}\,d{\mathbf{q}}\frac{1}{({\mathbf{p}}_{1}-{\mathbf{p}}_{2})^{2}}\frac{D_{0}(q,\nu)}{q^{2}(q^{2}-4\pi
e^{2}D_{0}(q,\nu))}\times$ $\displaystyle\hskip
10.0pt\left[-\frac{(n_{{\mathbf{p}}_{2}+{\mathbf{q}}}-n_{{\mathbf{p}}_{2}})(n_{{\mathbf{p}}_{1}+{\mathbf{q}}}-n_{{\mathbf{p}}_{1}})}{(-i\nu+\varepsilon_{{\mathbf{p}}_{1}+{\mathbf{q}}}-\varepsilon_{{\mathbf{p}}_{1}})^{2}}+\frac{2\beta
n_{{\mathbf{p}}_{1}}(1-n_{{\mathbf{p}}_{1}})n_{{\mathbf{p}}_{2}}}{-i\nu+\varepsilon_{{\mathbf{p}}_{1}+{\mathbf{q}}}-\varepsilon_{{\mathbf{p}}_{1}}}\right]$
$\displaystyle\xrightarrow[T\to 0]{|{\mathbf{q}}|\to 0}$
$\displaystyle(-Vm)\,\frac{(4\pi
e^{2})^{3}}{(2\pi)^{9}}\int_{-\infty}^{\infty}\frac{du}{2\pi}\int_{-1}^{1}\frac{\cos\vartheta_{1}\,\cos\vartheta_{2}\,d\cos\vartheta_{1}\,d\cos\vartheta_{2}}{(-iu+p_{F}\cos\vartheta_{1})^{2}}\times$
$\displaystyle\hskip
30.0pt\frac{(2\pi)^{2}p_{F}^{2}/2}{|\cos\vartheta_{1}-\cos\vartheta_{2}|}\int
4\pi qdq\frac{D_{0}(q,uq/m)}{q^{2}-4\pi e^{2}D_{0}(q,uq/m)}$
$\displaystyle+(2V)\frac{(4\pi
e^{2})^{3}}{(2\pi)^{3}}\int_{-\infty}^{\infty}\frac{du}{2\pi}\int_{-1}^{1}\frac{4\pi
mp_{F}\,d\\!\cos\vartheta}{-iu+p_{F}\cos\vartheta}\int\frac{dq}{(2\pi)^{2}}\frac{D_{0}(q,uq/m)}{q^{2}-4\pi
e^{2}D_{0}(q,uq/m)}\times$ $\displaystyle\hskip
15.0pt\left[\int\frac{q^{\prime
2}dq^{\prime}d\\!\cos\vartheta^{\prime}}{(2\pi)^{2}}\frac{1}{({\mathbf{q}}^{\prime})^{2}}\theta(-\cos\vartheta^{\prime}-\frac{q^{\prime}}{2p_{F}})\right]\;$
$\displaystyle\equiv$ $\displaystyle(-Vm)L_{1}+(2V)L_{2}\;,$

where the first term inside the square brackets after the second equal sign is
obtained by adding to it an equivalent expression with
${\mathbf{p}}_{1}+{\mathbf{q}}\to-{\mathbf{p}}_{1}$,
${\mathbf{p}}_{2}+{\mathbf{q}}\to-{\mathbf{p}}_{2}$, $\nu\to-\nu$, and then
taking the average, while the second term inside the same square brackets
results from changing $\nu\to-\nu$. The $L_{2}$ part, where the dummy variable
of the integral is switched from ${\mathbf{p}}_{2}$ to
${\mathbf{q}}^{\prime}\equiv{\mathbf{p}}_{2}-{\mathbf{p}}_{1}$, can be
cancelled by one of the terms contributing to the (b) diagram of (130). The
$L_{1}$ part, however, upon taking the $T\to 0$ and $q\to 0$ limits, may be
combined with (143) to yield a finite expression

$\displaystyle\frac{Vmp_{F}^{3}}{2}\,\frac{(4\pi
e^{2})^{3}}{(2\pi)^{7}}\int_{-\infty}^{\infty}\frac{du}{2\pi}\int_{-1}^{1}\frac{\cos\vartheta_{1}\,d\cos\vartheta_{1}}{(-iu+p_{F}\cos\vartheta_{1})^{2}}\,\frac{\cos\vartheta_{2}\,d\cos\vartheta_{2}}{-iu+p_{F}\cos\vartheta_{2}}\times$
$\displaystyle\hskip
10.0pt\frac{\cos\vartheta_{1}-\cos\vartheta_{2}}{|\cos\vartheta_{1}-\cos\vartheta_{2}|}\int
2\pi\;d(q^{2})\frac{D_{0}(q,uq/m)}{q^{2}-4\pi e^{2}D_{0}(q,uq/m)}\;.$ (145)

Note that this expression, in agreement with Carr and Maradudin,Carr and
Maradudin (1964) carries a different sign when compared to the original result
of DuBois.DuBois (1959)

It was conjectured before Valiev and Fernando (1997) that cancellation of
diagrams of the following form always hold true for HEG

$\begin{picture}(30.0,15.0)(0.0,-4.0) \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){}
\end{picture}\;\;\;+\;\;\;\;\;\begin{picture}(80.0,15.0)(0.0,-4.0)
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \end{picture}\;=\;0\;.\vspace*{4pt}$

Here, black circles denote parts of the diagram that are connected to each
other via two propagators. If this is the case, then the contributions of the
(b) diagrams of (129) and (130) will cancel each other. We don’t expect this
to happen since the (b) diagram of (129) is not anomalous. To illustrate that
the (b) diagram of (130) does not eliminate the (b) diagram of (129), we now
proceed to evaluate the (b) diagram of (130).

Here, $B_{2}({\mathbf{q}}^{\prime},\nu^{\prime})=4\pi
e^{2}/({\mathbf{q}}^{\prime})^{2}$ while $B_{1}({\mathbf{q}},\nu)=(4\pi
e^{2})^{2}D_{0}({\mathbf{q}},\nu)/\left\\{{\mathbf{q}}^{2}[{\mathbf{q}}^{2}-4\pi
e^{2}D_{0}({\mathbf{q}},\nu)]\right\\}$. Summing over $\nu^{\prime}$ via
methods described in section III.5, one obtains

$\displaystyle-\frac{1}{\beta}\;\begin{picture}(90.0,40.0)(0.0,-4.0)
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \end{picture}$ $\displaystyle=$ $\displaystyle-(2V)\frac{(4\pi
e^{2})^{3}/(2\pi)^{12}}{2\beta\int\frac{d{\mathbf{p}}}{(2\pi)^{3}}n_{\mathbf{p}}(1-n_{\mathbf{p}})}\frac{2}{\beta}\sum_{\nu}\int
d{\mathbf{p}}_{1}d{\mathbf{p}}_{2}\,d{\mathbf{q}}\,d{\mathbf{q}}^{\prime}\frac{n_{{\mathbf{p}}_{2}+{\mathbf{q}}^{\prime}}}{({\mathbf{q}}^{\prime})^{2}}\times$
(146) $\displaystyle\frac{D_{0}(q,\nu)\beta
n_{{\mathbf{p}}_{2}}(1-n_{{\mathbf{p}}_{2}})}{q^{2}(q^{2}-4\pi
e^{2}D_{0}(q,\nu))}\left[\frac{\beta
n_{{\mathbf{p}}_{1}}(1-n_{{\mathbf{p}}_{1}})}{-i\nu+\varepsilon_{{\mathbf{p}}_{1}+{\mathbf{q}}}-\varepsilon_{{\mathbf{p}}_{1}}}+\frac{n_{{\mathbf{p}}_{1}+{\mathbf{q}}}-n_{{\mathbf{p}}_{1}}}{(-i\nu+\varepsilon_{{\mathbf{p}}_{1}+{\mathbf{q}}}-\varepsilon_{{\mathbf{p}}_{1}})^{2}}\right]$
$\displaystyle\xrightarrow[T\to 0]{|{\mathbf{q}}|\to 0}$
$\displaystyle\frac{-(2V)(4\pi
e^{2})^{3}/(2\pi)^{3}}{2\int\frac{d{\mathbf{p}}}{(2\pi)^{3}}\delta(\mu-\varepsilon_{\mathbf{p}})}\int\frac{du}{2\pi}\int
d{\mathbf{p}}_{1}\frac{d{\mathbf{q}}}{(2\pi)^{3}}\frac{q}{m}\frac{D_{0}(q,uq/m)}{q^{2}(q^{2}-4\pi
e^{2}D_{0}(q,uq/m))}\times$ $\displaystyle\hskip
5.0pt\left[2\int\frac{d{\mathbf{p}}_{2}\,d{\mathbf{q}}^{\prime}}{(2\pi)^{6}}\delta(\mu-\varepsilon_{{\mathbf{p}}_{2}})\frac{n_{{\mathbf{p}}_{2}+{\mathbf{q}}^{\prime}}}{({\mathbf{q}}^{\prime})^{2}}\right]\times$
$\displaystyle\hskip
5.0pt\left[\frac{\delta(\mu-\varepsilon_{{\mathbf{p}}_{1}})m/q}{-iu+p_{F}\cos\vartheta_{1}}+\frac{\delta(\mu-\varepsilon_{{\mathbf{p}}_{1}})m^{2}/q^{2}}{(-iu+p_{F}\cos\vartheta_{1})^{2}}\frac{qp_{F}\cos\vartheta_{1}}{m}\right]$
$\displaystyle=$ $\displaystyle\frac{-(2V)(4\pi
e^{2})^{3}/(2\pi)^{3}}{2\int\frac{d{\mathbf{p}}}{(2\pi)^{3}}\delta(\mu-\varepsilon_{\mathbf{p}})}\int\frac{du}{2\pi}\int
d{\mathbf{p}}_{1}\frac{d{\mathbf{q}}}{(2\pi)^{3}}\frac{D_{0}(q,uq/m)}{q^{2}(q^{2}-4\pi
e^{2}D_{0}(q,uq/m))}\times$ $\displaystyle\hskip
5.0pt\left[2\int\frac{d{\mathbf{p}}_{2}}{(2\pi)^{3}}\delta(\mu-\varepsilon_{{\mathbf{p}}_{2}})\right]\left[\int\frac{d{\mathbf{q}}^{\prime}}{(2\pi)^{3}}\frac{\theta(-\cos\vartheta^{\prime}-\frac{q^{\prime}}{2p_{F}})}{({\mathbf{q}}^{\prime})^{2}}\right]\times$
$\displaystyle\hskip
5.0pt\left[\frac{\delta(p_{1}-p_{F})m/p_{F}}{-iu+p_{F}\cos\vartheta_{1}}+\frac{\delta(p_{1}-p_{F})m\cos\vartheta_{1}}{(-iu+p_{F}\cos\vartheta_{1})^{2}}\right]$
$\displaystyle=$ $\displaystyle-(2V)\frac{(4\pi
e^{2})^{3}}{(2\pi)^{3}}\int_{-\infty}^{\infty}\frac{du}{2\pi}\int_{-1}^{1}\frac{4\pi
mp_{F}d\cos\vartheta_{1}}{-iu+p_{F}\cos\vartheta_{1}}\left(1+\frac{p_{F}\cos\vartheta_{1}}{-iu+p_{F}\cos\vartheta_{1}}\right)\times$
$\displaystyle\hskip
5.0pt\left[\int\frac{dq}{(2\pi)^{2}}\frac{D_{0}(q,uq/m)}{q^{2}-4\pi
e^{2}D_{0}(q,uq/m)}\right]\left[\int\frac{d{\mathbf{q}}^{\prime}}{(2\pi)^{3}}\frac{\theta(-\cos\vartheta^{\prime}-\frac{q^{\prime}}{2p_{F}})}{({\mathbf{q}}^{\prime})^{2}}\right]\;.$

The contribution of (146) can be easily divided in two by explicitly expanding
the two parts inside the round parentheses. It is obvious that the
contribution associated with $1$ cancels exactly the $(2V)L_{2}$ part of
(IV.5) while the contribution associated with
$\frac{p_{F}\cos\vartheta_{1}}{-iu+p_{F}\cos\vartheta_{1}}$ cannot cancel the
$(-Vm)L_{1}$ part of (IV.5).

For the (c) diagrams of (128-130), both $B_{1}$ and $B_{2}$ are of order
$(4\pi e^{2})^{2}$, leading to contributions of order $(e^{2})^{4}$ and
higher. The (c) diagrams thus already lead us beyond what was studied by
DuBois DuBois (1959) and by Carr and Maradudin.Carr and Maradudin (1964) The
last two diagrams of $\Gamma_{2}[n]$ (shown in Figure 3) gives rise to the
$E^{\prime}_{3}$ term of Carr and Maradudin Carr and Maradudin (1964) when the
$\tilde{\mathcal{D}}_{0}$ lines are each replaced by $U$, the first term in
the decomposition of $\tilde{\mathcal{D}}_{0}$. In principle, one may go on to
study terms of order $(e^{2})^{4}$; we will not, however, delve into this
endeavor since this is not the primary aim here. We would like to emphasize
the following points. First, unlike conventional $e^{2}$ based perturbation
theory, the formalism presented here naturally avoids divergence. This is
shown by the fact that each diagram in our formalism contains no singularity
while attempts to perform $e^{2}$ based purturbation necessarily require
further re-grouping, such as combining the (b) diagrams of (128) and (129), to
tame the divergence. Second, even if one were to pursue $e^{2}$ expansion,
using Dyson’s equation (127) within our formalism still makes the task
straightforward. In fact, as shown in this section, the
$\tilde{\Gamma}_{0}[n]+\Gamma_{1}[n]+\Gamma_{2}[n]$ part when applied to the
HEG already contain the celebrated results of Carr and Maradudin.Carr and
Maradudin (1964) Third, the removal of anomalous contributions that required
special attention within the formalism of Luttinger et al. Kohn and Luttinger
(1960); Luttinger and Ward (1960) becomes automatic under this formalism.

##  V Excitations

To obtain information regarding the excitations, one needs a time-dependent
probe, although one with infinitesimal amplitude is sufficient. Runge and
Gross Runge and Gross (1984) extended the correspondence between the external
probe potential and the ground state charge density of the DFT to the time-
dependent case. This relationship provides the foundation for studying
excitation energies under DFT. Fukuda et al. Fukuda et al. (1994) expressed
the excitation energy condition using effective action formalism, without
explicit connection to the Kohn-Sham formalism. By introducing time-dependent
Kohn-Sham orbitals while assuming time idependence of the orbital occupation
numbers, Casida Casida (1995b) derived via linear response theory a self-
consistent condition on the density matrix response that leads to
determination of excitation energies. Along a similar line, Petersilka et al.
Petersilka et al. (1996) proposed the so-called optimized effective potential
(expanded in time-dependent Kohn-Sham orbitals) to tackle the problem of
excitation energies. Since there exist formalisms to extract the excitation
energies of the system provided that the UDF is known, our result for
excitation energies should not be considered novel. The reasons for this
section are twofold. First, we would like to explicitly show that the
excitation energies can be obtained using the formalism of section III without
introducing time-dependent orbitals. Second, although it is possible to find
derivations of formulas Valiev and Fernando (1997) similar to those that will
be shown here, they do appear slightly different and thus a self-contained
exposition may be helpful.

Intuitively speaking, by varying the frequency of the probe, one seeks the
frequency/energy where the amplitude of the response function diverges.
Indeed, it is known that the spectral representation of the correlation
function has poles at the excitation energies of the system Kobe (1962).
Having obtained the effective action $\Gamma[n]$, we also note that the second
(functional) derivative of $\Gamma[n]$ with respect to the local electron
density is the inverse of the density-density correlation function. Therefore,
any pole associated with the correlation function becomes a root of the
effective action. It can be shown that upon analytic continuation the
correlation function, obtained using the imaginary time (finite temperature)
formalism, can be turned into the response function of the real time. Below we
briefly illustrate this point. Readers interested in more details can find an
extensive exposition in reference Fetter and Walecka, 1971.

From Eqs. (13) and (28), we know that

$\displaystyle n(x)$ $\displaystyle=$ $\displaystyle{\delta(\beta
W[J])\over\delta J(x)}\;,$ $\displaystyle{\rm and\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ }J(x)$
$\displaystyle=$ $\displaystyle-\frac{\delta\Gamma[n]}{\delta n(x)}\;.$

One then considers

$\delta(x-y)=\frac{\delta J(x)}{\delta J(y)}=-\int
dz\frac{\delta^{2}\Gamma[n]}{\delta n(x)\delta n(z)}\frac{\delta n(z)}{\delta
J(y)}=-\int dz\frac{\delta^{2}\Gamma[n]}{\delta n(x)\delta
n(z)}\frac{\delta^{2}(\beta W[J])}{\delta J(z)\delta J(y)}\;.$ (147)

Although a time-dependence of $J$ is introduced to probe the excitations, in
the end we will return to a time-independent source ($J(x)\to J({\bf x})$)
while computing the excitation energies. As we will show below, it is most
convenient to go to the zero temperature limit to compute the excitation
energies.

Note that

$\displaystyle\frac{\delta^{2}(\beta W[J])}{\delta J(x)\delta J(y)}$
$\displaystyle=$
$\displaystyle-\left[\hat{n}(x)\hat{n}(y)\right]_{T}+\left[\hat{n}(x)\right]_{T}\left[\hat{n}(y)\right]_{T}=-\left[\hat{n}(x)\hat{n}(y)\right]_{T}+n(x)n(y)$
(148) $\displaystyle=$
$\displaystyle-\left[(\hat{n}(x)-n(x))(\hat{n}(y)-n(y))\right]_{T}\equiv-\left[\tilde{n}(x)\,\tilde{n}(y)\right]_{T}$

where $\hat{n}(x)=\psi^{{\dagger}}(x)\psi(x)$ is the electron density operator
and the square bracket $[]_{T}$ indicates an (imaginary) time-ordered thermal
average such that

$\left[\hat{\mathcal{O}}_{1}(t_{1})\hat{\mathcal{O}}_{2}(t_{2})\right]_{T}=\frac{\text{Tr}\left[e^{-\beta
H[J]}T\left(\hat{\mathcal{O}}_{1}(t_{1})\hat{\mathcal{O}}_{2}(t_{2})\right)\right]}{\text{Tr}\left[e^{-\beta
H[J]}\right]}=\frac{\text{Tr}\left[e^{-\beta
H[J]}T\left(\hat{\mathcal{O}}_{1}(t_{1})\hat{\mathcal{O}}_{2}(t_{2})\right)\right]}{Z[J]}\;.$

Let us denote by $\\{|\ell\rangle_{ex}\\}_{\ell=0}^{\infty}$ the eigenstates
of the Hamiltonian $H[J]$ with the corresponding eigenenergies
$\\{{\mathcal{E}}_{\ell}\\}_{\ell=0}^{\infty}$. The spectral representation is
obtained by first writing the time-ordered product (assuming that operators
$\hat{\mathcal{O}}_{1}(t_{1})$ and $\hat{\mathcal{O}}_{2}(t_{2})$ are bosonic)
as

$\phantom{}{}_{ex}\\!\langle\ell|T\left(\hat{\mathcal{O}}_{1}(t_{1})\hat{\mathcal{O}}_{2}(t_{2})\right)|\ell\rangle_{ex}=\theta(t_{1}-t_{2})\,\phantom{}_{ex}\\!\langle\ell|\hat{\mathcal{O}}_{1}(t_{1})\hat{\mathcal{O}}_{2}(t_{2})|\ell\rangle_{ex}+\theta(t_{2}-t_{1})\,\phantom{}_{ex}\\!\langle\ell|\hat{\mathcal{O}}_{2}(t_{2})\hat{\mathcal{O}}_{1}(t_{1})|\ell\rangle_{ex}\;,$

then by inserting the identity operator
$\sum_{\ell\,^{\prime}}|\ell\,^{\prime}\rangle_{ex}\;\phantom{}{}_{ex}\\!\langle\ell\,^{\prime}|$
between the two operators, and finally by multiplying by
$e^{i\omega(t_{1}-t_{2})}$ and then integrating over the time variable
$t_{1}-t_{2}$. Proceeding in this way, one will then obtain information on
${\mathcal{E}}_{\ell\,^{\prime}}-{\mathcal{E}}_{\ell}$. Since knowing
$\Omega_{\ell\,^{\prime}}\equiv{\mathcal{E}}_{\ell\,^{\prime}}-{\mathcal{E}}_{0}$
also provides complete information on
${\mathcal{E}}_{\ell\,^{\prime}}-{\mathcal{E}}_{\ell}$, one may also focus on
$\ell=0$ by taking the limit $\beta\to\infty$. Let

$-W^{(2)}(x,y)\equiv-\frac{\delta^{2}(\beta W[J])}{\delta J(x)\delta J(y)}\;,$
(149)

and

$-W^{(2)}({\bf x},{\bf
y},i\nu_{n})\equiv\int_{0}^{\beta}d(\tau_{x}-\tau_{y})e^{i\nu_{n}(\tau_{x}-\tau_{y})}\;W^{(2)}(x,y)\equiv\int_{-\infty}^{\infty}\frac{d\omega^{\prime}}{2\pi}\frac{A(\omega^{\prime})}{i\nu_{n}-\omega^{\prime}}\;,$
(150)

with $\nu_{n}=2\pi n/\beta$ and

$A(\omega)=e^{\beta
W[J]}\sum_{\ell,m}e^{-\beta{\mathcal{E}}_{\ell}}\left(e^{-\beta\omega}-1\right)2\pi\delta(\omega+{\mathcal{E}}_{\ell}-{\mathcal{E}}_{m})\,\phantom{}_{ex}\\!\langle\ell|\tilde{n}({\bf
x})|m\rangle_{ex}\,\phantom{}{}_{ex}\\!\langle m|\tilde{n}({\bf
y})|\ell\rangle_{ex}\;.$ (151)

Note that $\tilde{n}({\bf x})$ measures the deviation from the thermally
averaged electronic density $n_{T}({\bf x})$. Since the expression (149) is
evaluated at static $J({\bf x})$, the Hamiltonian contains no time dependence.
We may thus write $\tilde{n}(x)=e^{H\tau_{x}}\tilde{n}({\bf
x})e^{-H\tau_{x}}$.

Let’s now express using real time the retarded correlation function
(${\mathcal{R}}$), also called the response function, and the advanced
correlation function (${\mathcal{A}}$) as follows (with
$n(x)=e^{it_{x}H}n({\bf x})e^{-it_{x}H}$)

$\displaystyle\left(\begin{array}[]{c}{\cal R}(x,y)\\\
{\mathcal{A}}(x,y)\end{array}\right)$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{r}-i\theta(t_{x}-t_{y})\\\
i\theta(t_{y}-t_{x})\end{array}\right)e^{\beta W[J]}\text{Tr}\left(e^{-\beta
H}\left[\tilde{n}(x),\tilde{n}(y)\right]\right)$ (156) $\displaystyle=$
$\displaystyle\left(\begin{array}[]{r}-i\theta(t_{x}-t_{y})\\\
i\theta(t_{y}-t_{x})\end{array}\right)\sum_{\ell,m}\,e^{\beta(W[J]-{\mathcal{E}}_{\ell})}\left[e^{i(t_{x}-t_{y})({\mathcal{E}}_{\ell}-{\mathcal{E}}_{m})}\,\phantom{}_{ex}\\!\langle\ell|\tilde{n}({\bf
x})|m\rangle_{ex}\,\phantom{}{}_{ex}\\!\langle m|\tilde{n}({\bf
y})|\ell\rangle_{ex}\right.$ (160) $\displaystyle\left.\hskip
90.0pt-e^{-i(t_{x}-t_{y})({\mathcal{E}}_{\ell}-{\mathcal{E}}_{m})}\phantom{}_{ex}\\!\langle\ell|\tilde{n}({\bf
y})|m\rangle_{ex}\,\phantom{}{}_{ex}\\!\langle m|\tilde{n}({\bf
x})|\ell\rangle_{ex}\right]\;.$

Taking the Fourier transform of the response function, we consider

$\displaystyle\left(\\!\\!\begin{array}[]{c}{\cal R}({\bf x},{\bf
y},\omega)\\\ {\mathcal{A}}({\bf x},{\bf y},\omega)\end{array}\\!\\!\right)$
$\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}d(t_{x}-t_{y})\;e^{i\omega(t_{x}-t_{y})}\
\left(\begin{array}[]{c}{\cal R}(x,y)\\\ {\mathcal{A}}(x,y)\end{array}\right)$
$\displaystyle=$ $\displaystyle e^{\beta
W[J]}\sum_{\ell,m}e^{-\beta{\mathcal{E}}_{\ell}}\frac{e^{\beta({\mathcal{E}}_{\ell}-{\mathcal{E}}_{m})}-1}{\omega+{\mathcal{E}}_{\ell}-{\mathcal{E}}_{m}\pm
i\eta}\,\phantom{}_{ex}\\!\langle\ell|\tilde{n}({\bf
x})|m\rangle_{ex}\,\phantom{}{}_{ex}\\!\langle m|\tilde{n}({\bf
y})|\ell\rangle_{ex}$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\frac{d\omega^{\prime}}{2\pi}\frac{A(\omega^{\prime})}{\omega-\omega^{\prime}\pm
i\eta}\;.$ (166)

Comparing Eqs. (150) and (166), one finds that substituting
$i\nu_{n}\to\omega+i\eta$ ($\omega-i\eta$) in the imaginary time-ordered
correlation function leads to the retarded (advanced) correlation function.
The validity of this analytic continuation was discussed by Baym and Mermin
Baym and Mermin (1961).

As $T\to 0$, $e^{-\beta
W[J]}=e^{-\beta{\mathcal{E}}_{0}}\left[1+{\mathcal{O}}(e^{-\beta({\mathcal{E}}_{1}-{\mathcal{E}}_{0})})\right]$,
and $e^{\beta
W[J]}=e^{\beta{\mathcal{E}}_{0}}\left[1-{\mathcal{O}}(e^{-\beta({\mathcal{E}}_{1}-{\mathcal{E}}_{0})})\right]$.
Under this limit, we may rewrite the spectral weight $A(\omega)$ as

$\displaystyle\lim_{\beta\to\infty}\frac{A(\omega)}{2\pi}$ $\displaystyle=$
$\displaystyle\sum_{\ell,m}\left(e^{-\beta({\mathcal{E}}_{m}-{\mathcal{E}}_{0})}-e^{-\beta({\mathcal{E}}_{\ell}-{\mathcal{E}}_{0})}\right)\delta(\omega+{\mathcal{E}}_{\ell}-{\mathcal{E}}_{m})\phantom{}_{ex}\\!\langle\ell|\tilde{n}({\bf
x})|m\rangle_{ex}\,\phantom{}{}_{ex}\\!\langle m|\tilde{n}({\bf
y})|\ell\rangle_{ex}$ (167)
$\displaystyle+{\mathcal{O}}(e^{-\beta({\mathcal{E}}_{1}-{\mathcal{E}}_{0})})$
$\displaystyle=$
$\displaystyle\sum_{\ell}\,\left[\,\delta(\omega+{\mathcal{E}}_{\ell}-{\mathcal{E}}_{0})\phantom{}_{ex}\\!\langle
0|\tilde{n}({\bf
y})|\ell\rangle_{ex}\,\phantom{}{}_{ex}\\!\langle\ell|\tilde{n}({\bf
x})|0\rangle_{ex}\right.$ $\displaystyle\hskip
10.0pt\left.-\delta(\omega-({\mathcal{E}}_{\ell}-{\mathcal{E}}_{0}))\phantom{}_{ex}\\!\langle
0|\tilde{n}({\bf
x})|\ell\rangle_{ex}\,\phantom{}{}_{ex}\\!\langle\ell|\tilde{n}({\bf
y})|0\rangle_{ex}\right]$ $\displaystyle\equiv$
$\displaystyle\sum_{\ell}\left[\delta(\omega+{\mathcal{E}}_{\ell}-{\mathcal{E}}_{0})n_{\ell,e}^{*}({\bf
y})n_{\ell,e}({\bf
x})-\delta(\omega-({\mathcal{E}}_{\ell}-{\mathcal{E}}_{0}))n_{\ell,e}^{*}({\bf
x})n_{\ell,e}({\bf y})\right]$

where $n_{\ell,e}({\bf y})\equiv\phantom{}_{ex}\\!\langle\ell|\tilde{n}({\bf
y})|0\rangle_{ex}$ .

When continued to the retarded correlation function (response function), the
density correlation function $W^{(2)}$ reads

$\lim_{\beta\to\infty}W^{(2)}({\bf x},{\bf
y},\omega)=\sum_{\ell}\left[\frac{n^{*}_{\ell,e}({\bf x})\,n_{\ell,e}({\bf
y})}{\omega-\Omega_{\ell}+i\eta}-\frac{n^{*}_{\ell,e}({\bf
y})\,n_{\ell,e}({\bf x})}{\omega+\Omega_{\ell}+i\eta}\right]\;,$ (168)

with

$\Omega_{\ell}\equiv{\mathcal{E}}_{\ell}-{\mathcal{E}}_{0}\;.$ (169)

Provided that the amplitudes $\phantom{}{}_{ex}\\!\langle
0|\tilde{n}|\ell\rangle_{ex}$ are nonzero, we see from Eq. (168) that
$\omega+i\eta=\pm\left({\mathcal{E}}_{l}-{\mathcal{E}}_{0}\right)$ are simple
poles of $W^{(2)}({\bf x},{\bf y},\omega)$. Furthermore, we also see that

$\left(W^{(2)}({\bf x},{\bf y},\omega)\right)^{*}=W^{(2)}({\bf x},{\bf
y},-\omega^{*})\;,$ (170)

and when $\omega$ is real

$\left(W^{(2)}({\bf x},{\bf y},\omega)\right)^{*}=W^{(2)}({\bf x},{\bf
y},-\omega)\;.$

Let us also define

$\Gamma^{(2)}(x,y)\equiv\frac{\delta^{2}\Gamma[n]}{\delta n(x)\delta n(y)}\;.$

Eq. (147) may thus be rewritten as

$-\delta(x-y)=\int dz\;\Gamma^{(2)}(x,z)\,W^{(2)}(z,y)\;.$ (171)

Since eventually, $\tau_{x}$ must agree with $\tau_{y}$ in the equation above,
if $\tau_{z}>\tau_{y}$, we must have $\tau_{z}>\tau_{x}$ as well. Similarly to
Eq. (150), the Fourier transform for $\Gamma^{(2)}$can be written as

$\Gamma^{(2)}({\bf x},{\bf
z},i\nu_{n})\equiv\int_{0}^{\beta}\\!\\!d(\tau_{z}-\tau_{x})\,e^{i\nu_{n}(\tau_{z}-\tau_{x})}\,\Gamma^{(2)}(x,z)\;.$
(172)

The inverse transform of (150) and (172) can be written as

$\displaystyle W^{(2)}(z,y)$ $\displaystyle=$
$\displaystyle\frac{1}{\beta}\sum_{n}e^{-i\nu_{n}(\tau_{z}-\tau_{y})}W^{(2)}({\bf
z},{\bf y},i\nu_{n})$ $\displaystyle\Gamma^{(2)}(x,z)$ $\displaystyle=$
$\displaystyle\frac{1}{\beta}\sum_{n}e^{-i\nu_{n}(\tau_{z}-\tau_{x})}\Gamma^{(2)}({\bf
x},{\bf z},i\nu_{n})$

and

$\int_{0}^{\beta}\\!\\!\\!\\!d\tau_{z}\\!\\!\int\\!\\!d{\bf
z}\;\Gamma^{(2)}(x,z)W^{(2)}(z,y)=\int d{\bf
z}\,\frac{1}{\beta}\sum_{n}\,e^{-i\nu_{n}(\tau_{x}-\tau_{y})}\Gamma^{(2)}({\bf
x},{\bf z},-i\nu_{n})\,W^{(2)}({\bf z},{\bf y},i\nu_{n})\;.$

Since
$\delta(\tau_{x}-\tau_{y})=\frac{1}{\beta}\sum_{n}e^{-i\nu_{n}(\tau_{x}-\tau_{y})}$,
one obtains

$-\delta({\bf x}-{\bf y})=\int\\!\\!d{\bf z}\;\Gamma^{(2)}({\bf x},{\bf
z},-i\nu_{n})W^{(2)}({\bf z},{\bf y},i\nu_{n})=\int\\!\\!d{\bf
z}\;W^{(2)}({\bf x},{\bf z},-i\nu_{n})\Gamma^{(2)}({\bf z},{\bf
y},i\nu_{n})\;,$ (173)

which under analytic continuation becomes

$-\delta({\bf x}-{\bf y})=\int\\!\\!d{\bf z}\;\Gamma^{(2)}({\bf x},{\bf
z},-\omega)W^{(2)}({\bf z},{\bf y},\omega)=\int\\!\\!d{\bf z}\;W^{(2)}({\bf
x},{\bf z},-\omega)\Gamma^{(2)}({\bf z},{\bf y},\omega)\;.$ (174)

From Eqs. (170) and (174), one has

$\left(\Gamma^{(2)}({\bf x},{\bf z},\omega)\right)^{*}=\Gamma^{(2)}({\bf
x},{\bf z},-\omega^{*})\;.$ (175)

Multiplying the LHS and the middle expression of Eq. (174) by
$(\omega\mp\Omega_{\ell}+i\eta)$ and then setting
$\omega\to\pm\,\Omega_{\ell}-i\eta$, we see that

$\displaystyle\int\\!\\!d{\bf z}\,\Gamma^{(2)}({\bf x},{\bf
z},-\omega\to-\Omega_{\ell}-i\eta)\,n^{*}_{\ell,e}({\bf z})$ $\displaystyle=$
$\displaystyle 0\;,$ (176) $\displaystyle{\rm and\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
}\int\\!\\!d{\bf z}\,\Gamma^{(2)}({\bf x},{\bf
z},-\omega\to+\Omega_{\ell}-i\eta)\,n_{\ell,e}({\bf z})$ $\displaystyle=$
$\displaystyle 0\;.\;,$ (177)

That is, $n_{\ell,e}^{(*)}({\bf y})$ become eigenvectors of $\Gamma^{(2)}({\bf
x},{\bf y},-\omega)$ with zero eigenvalues. The matter of finding excitation
energies and the corresponding electronic densities thus reduces to finding
for $\Gamma^{(2)}$ the eigenvectors with zero eigenvalue.Fukuda et al. (1994)
Since we are interested in obtaining the excitation energy under the physical
condition, $J({\bf x})\to 0$, this also means that the derivatives of $\Gamma$
above are evaluated at the ground state electronic density in the zero
temperature limit.

From the expressions (65) and (70), one sees that the effective action is
split into the free particle part $\Gamma_{0}$, the Hartree functional and the
exchange-correlational functional

$\Gamma[n]=\Gamma_{0}[n]+\frac{1}{2}n{\scriptstyle\circ}\,U{\scriptstyle\circ}n+\Gamma_{xc}[n]\equiv\Gamma_{0}[n]+\Gamma_{\rm
int}[n],$ (178)

where $\Gamma_{xc}[n]=\sum_{l=1}^{\infty}\Gamma_{l}[n]$. Letting $\Delta(x)$
be the eigenvector, the excitation condition becomes

$\int d{\bf z}\Gamma^{(2)}({\bf x},{\bf z})\Delta({\bf z})=0\;.$ (179)

After splitting the effective action into $\Gamma_{0}$ and $\Gamma_{\rm int}$
the eigenvalue equations (176-177) may be expressed as

$-\int\\!\\!d{\bf z}\,\Gamma_{0}^{(2)}({\bf y},{\bf z},-\omega)\Delta({\bf
z})=\int\\!\\!d{\bf z}\,\Gamma_{\rm int}^{(2)}({\bf y},{\bf
z},-\omega)\Delta({\bf z})\;.$ (180)

Multiplying both sides of (180) by $W_{0}^{(2)}({\bf x},{\bf y},\omega)$ and
then integrating over $d{\bf y}$, one obtains

$\Delta({\bf x})=\int d{\bf y}d{\bf z}\;W_{0}^{(2)}({\bf x},{\bf
y},\omega)\Gamma_{\rm int}^{(2)}({\bf y},{\bf z},-\omega)\Delta({\bf z})\;.$
(181)

Note that $W_{0}[J_{0}]$ describes our Kohn-Sham system, a constructed non-
interacting system that produces the same ground state electron density as
that of the physical system considered. From Eq. (168), one can write down
$W^{(2)}_{0}({\bf x},{\bf y},\omega)$ in terms of the excitation energies
associated with the Kohn-Sham non-interacting system:

$W^{(2)}_{0}({\bf x},{\bf y},\omega)=\sum_{\ell}\left[\frac{n^{*}_{\ell}({\bf
x})\,n_{\ell}({\bf y})}{\omega-\omega_{\ell}+i\eta}-\frac{n^{*}_{\ell}({\bf
y})\,n_{\ell}({\bf x})}{\omega+\omega_{\ell}+i\eta}\right]\;,$ (182)

where

$\omega_{\ell}\equiv E_{\ell}-E_{0}\;,$ (183)

and $E_{\ell}$ is the energy of $|\ell\rangle_{ks}$, the $\ell$th state of the
many-particle Kohn-Sham system with $J_{0}({\bf x})$ chosen to generate the
correct ground state electron density. Therefore, for any $\ell$, $E_{\ell}$
is simply the sum of single-particle energies $\varepsilon_{m}$. Note that in
Eq. (182), $n_{\ell}({\bf x})$ is defined by
$\phantom{}{}_{ks}\\!\langle\ell|\tilde{n}({\bf x})|0\rangle_{ks}$ except that
$|\ell\rangle_{ks}$ now describes the $\ell$th state of the Kohn-Sham system,
not the physical system considered.

We seek a general solution for $\Delta({\bf x})$ of the form

$\Delta({\bf x})=\sum_{\ell}\,\left[a_{\ell}\,n_{\ell}({\bf
x})+b_{\ell}\,n^{*}_{\ell}({\bf x})\right]\;.$ (184)

Evidently, if a frequency $\hat{\omega}$ leads to a solution $\\{b_{\ell}\\}$,
then $-\hat{\omega}^{*}$ should lead to a solution $\\{a_{\ell}\\}$, which
plays the role of $\\{b_{\ell}^{*}\\}$. Substituting Eqs. (182) and (184) into
Eq. (181), we find that

$\displaystyle\sum_{\ell}\left[a_{\ell}n_{\ell}({\bf
x})+b_{\ell}n^{*}_{\ell}({\bf x})\right]$
$\displaystyle=\sum_{\ell,\ell\,^{\prime}}\int\\!\\!d{\bf y}d{\bf
z}\left[\frac{n^{*}_{\ell}({\bf x})n_{\ell}({\bf
y})}{\omega-\omega_{\ell}+i\eta}-\frac{n_{\ell}({\bf x})n^{*}_{\ell}({\bf
y})}{\omega+\omega_{\ell}+i\eta}\right]\Gamma^{(2)}_{\rm int}({\bf y},{\bf
z},-\omega)\left[a_{\ell\,^{\prime}}n_{\ell\,^{\prime}}({\bf
z})+b_{\ell\,^{\prime}}n^{*}_{\ell\,^{\prime}}({\bf z})\right].$ (185)

Equating the coefficients associated with $n_{\ell}({\bf x})$ and
$n_{\ell}^{*}({\bf x})$, we find that

$\displaystyle a_{\ell}$ $\displaystyle=$
$\displaystyle-\sum_{\ell\,^{\prime}}\int d{\bf y}d{\bf
z}\frac{n_{\ell}^{*}({\bf y})}{\omega+\omega_{\ell}+i\eta}\Gamma^{(2)}_{\rm
int}({\bf y},{\bf z},-\omega)\left[a_{\ell\,^{\prime}}n_{\ell\,^{\prime}}({\bf
z})+b_{\ell\,^{\prime}}n^{*}_{\ell\,^{\prime}}({\bf z})\right]\;,$
$\displaystyle b_{\ell}$ $\displaystyle=$
$\displaystyle\sum_{\ell\,^{\prime}}\int d{\bf y}d{\bf z}\frac{n_{\ell}({\bf
y})}{\omega-\omega_{\ell}+i\eta}\Gamma^{(2)}_{\rm int}({\bf y},{\bf
z},-\omega)\left[a_{\ell\,^{\prime}}n_{\ell\,^{\prime}}({\bf
z})+b_{\ell\,^{\prime}}n^{*}_{\ell\,^{\prime}}({\bf z})\right]\;.$

Let us define

$\displaystyle Y_{\ell,\ell\,^{\prime}}(\omega)$ $\displaystyle=$
$\displaystyle\int d{\bf y}d{\bf z}\;n_{\ell}^{*}({\bf y})\Gamma^{(2)}_{\rm
int}({\bf y},{\bf z},-\omega)n_{\ell\,^{\prime}}({\bf
z})\;+\omega_{\ell}\,\delta_{\ell,\ell\,^{\prime}}\;,$ (186) $\displaystyle
K_{\ell,\ell\,^{\prime}}(\omega)$ $\displaystyle=$ $\displaystyle\int d{\bf
y}d{\bf z}\;n_{\ell}^{*}({\bf y})\Gamma^{(2)}_{\rm int}({\bf y},{\bf
z},-\omega)n^{*}_{\ell\,^{\prime}}({\bf z})\;,\;$ (187)

and we obtain the following matrix equation

$\left(\begin{array}[]{cc}{\mathbf{Y}}(\omega)&{\mathbf{K}}(\omega)\\\
{\mathbf{K}}^{*}(-\omega^{*})&{\mathbf{Y}}^{*}(-\omega^{*})\end{array}\right)\left(\begin{array}[]{c}A\\\
B\end{array}\right)=\left(\omega+i\eta\right)\left(\begin{array}[]{rr}-1&0\\\
0&1\end{array}\right)\left(\begin{array}[]{c}A\\\ B\end{array}\;,\right)$
(188)

where $(A)_{\ell}=a_{\ell}$ and $(B)_{\ell}=b_{\ell}$. Evidently, one seeks
$\hat{\omega}$ such that

$\det\left[\left(\begin{array}[]{cc}{\mathbf{Y}}(\hat{\omega})&{\mathbf{K}}(\hat{\omega})\\\
{\mathbf{K}}^{*}(-\hat{\omega}^{*})&{\mathbf{Y}}^{*}(-\hat{\omega}^{*})\end{array}\right)-\left(\hat{\omega}+i\eta\right)\left(\begin{array}[]{rr}-1&0\\\
0&1\end{array}\right)\right]=0\;.$ (189)

As mentioned earlier, one anticipates
$\left(A(\hat{\omega})\right)=\left(B^{*}(-\hat{\omega}^{*})\right)$. To see
this, we perform the change $\omega\Leftrightarrow-\omega^{*}$ in (188) and
find that one can then rearrange the resulting equation into

$\left(\begin{array}[]{cc}{\mathbf{Y}}(\omega)&{\mathbf{K}}(\omega)\\\
{\mathbf{K}}^{*}(-\omega^{*})&{\mathbf{Y}}^{*}(-\omega^{*})\end{array}\right)\left(\begin{array}[]{c}B^{*}(-\omega^{*})\\\
A^{*}(-\omega^{*})\end{array}\right)=\left(\omega+i\eta\right)\left(\begin{array}[]{rr}-1&0\\\
0&1\end{array}\right)\left(\begin{array}[]{c}B^{*}(-\omega^{*})\\\
A^{*}(-\omega^{*})\end{array}\right)\;,$ (190)

which is identical to Eq. (188) except with
$\left(B^{*}(-\hat{\omega}^{*})\right)$ playing the role of
$\left(A(\hat{\omega})\right)$ and $\left(A^{*}(-\hat{\omega}^{*})\right)$
playing the role of $\left(B(\hat{\omega})\right)$ .

We now compare Eq. (188) with similar existing results. In references
Bauernschmitt and Ahlrichs, 1996 and Valiev and Fernando, 1997, equations
similar to (188) were obtained, and those will be identical to Eq. (188)
provided that ${\mathbf{K}}^{*}(-\omega^{*})={\mathbf{K}}^{*}(\omega)$. This
will happen if $\Gamma({\bf x},{\bf y},-\omega)=\Gamma({\bf x},{\bf
y},\omega)$ for real $\omega$.

##  VI Saddle-Point as an Alternative Formalism

Below, we will obtain the effective action using a classical variable
$i\varphi_{c}$ that corresponds to the saddle-point of the auxiliary field
path integral. At the physical condition $J=0$, $i\varphi_{c}$ is interpreted
as the electron density of a self-consistent Hartree solution. Our Hartee
problem is not of the conventional type, but rather similar to what Kohn
described in his Nobel lecture. Kohn (1999) In the conventional Hartree
calculation, the wave functions obtained may not be orthogonal to each other
due to the fact that each particle’s wave function is solved with a different
potential. Slater (1930) In the method below and mentioned by Kohn Kohn
(1999), the electric potential experienced by every electron is the same.
Another difference between the method below and the aforementioned Hartree
methods Slater (1930); Kohn (1999) is that the integral of the Hartree density
in our method is not necessarily an integer due to the possibility that the
density correction term may have a nonzero integral.

The saddle-point method below is quite different from what was described in
the previous sections. First, although the diagrams in the saddle-point method
are all connected diagrams, they are not one-particle irreducible (1PI).
Second, unlike the method presented in previous sections, the computation of
the effective action now requires no further functional derivatives of $\beta
W[J]$ with respect to $J$ evaluated at $J_{c}$, while in the formalism
mentioned in previous sections, one needs to compute higher order derivatives
of $\beta W[J_{0}]$ with respect to $J_{0}$ (see Eqs. (67-73) ).

###  VI.1 Evaluation of $e^{-\beta W_{\phi}[J]}$ via expansion around the
saddle-point

The path integral (23)

$e^{-\beta W_{\phi}[J]}\equiv\int
D\phi\;\exp\left\\{-I\left[\phi\right]-iJ{\scriptstyle\circ}\phi\right\\}\;.$

may be evaluated by the saddle point method. The extrema condition gives

$\left.\frac{\delta
I\left[\phi\right]}{\delta\phi(x)}\right|_{\varphi_{c}}=-iJ(x)\;.$

Since

$\displaystyle\left.\frac{\delta
I\left[\phi\right]}{\delta\phi(x)}\right|_{\varphi_{c}}$ $\displaystyle=$
$\displaystyle\left(U{\scriptstyle\circ}\varphi_{c}\right)_{x}-\text{Tr}\left(G_{\phi}\frac{\delta
G_{\phi}^{-1}}{\delta\phi(x)}\right)_{\phi\to\varphi_{c}}$ (191)
$\displaystyle=$
$\displaystyle\left(U{\scriptstyle\circ}\varphi_{c}\right)_{x}-\int
dy\;{\mathcal{G}}_{c}(y,y)\left(iU(y,x)\right)\;,$

where ${\mathcal{G}}_{c}(y,y)\equiv G_{\phi\to\varphi_{c}}(y,y)$, we obtain
with $U(x,y)=U(y,x)$

$\displaystyle J(x)=\left(U{\scriptstyle\circ}(i\varphi_{c})\right)_{x}+\int
dy\;U(x,y){\mathcal{G}}_{c}(y,y)=\int
dy\;U(x,y)\left[i\varphi_{c}(y)+{\mathcal{G}}_{c}(y,y)\right]\;.$ (192)

As $J\to 0$ (the physical condition), when $U$ is invertible such as Coulomb
interaction, one must have $i\varphi_{c}(x)=-{\mathcal{G}}_{c}(x,x)$. Note
that the negative of the diagonal element of the Green’s function
$-{\mathcal{G}}_{c}(x,x)$ is the particle density corresponding to the
following Hamiltonian

$\int
dx\;{\hat{\psi}}^{{\dagger}}(x)\left[-\frac{\nabla^{2}}{2m}+\upsilon_{\rm
ion}({{\bf
x}})-\mu+i(U{\scriptstyle\circ}\varphi_{c})_{x}\right]\hat{\psi}(x)\;.$

The saddle-point equation therefore produces a Hartree-like equation: the
Green’s function depends on the input particle density $i\varphi_{c}(x)$ and
is required to produce the same particle density $i\varphi_{c}(x)$ in the end.

When $J\neq 0$, one can still view (192) as a generalized Hartree equation in
the following sense. Remember that the inverse Green’s function
${\mathcal{G}}^{-1}_{c}(x,y)$ is given by (with
$\delta(x-x^{\prime})=\delta(\tau-\tau^{\prime})\delta({\bf x}-{\bf
x}^{\prime})$)

${\mathcal{G}}^{-1}_{c}(x,x^{\prime})=\left[\partial_{\tau}-\frac{\nabla^{2}}{2m}+\upsilon_{\rm
ion}({{\bf
x}})-\mu+U{\scriptstyle\circ}(i\varphi_{c})\right]\delta(x-x^{\prime})\;,$
(193)

and may be rewritten as

${\mathcal{G}}^{-1}_{c}(x,x^{\prime})=\left[\partial_{\tau}-\frac{\nabla^{2}}{2m}+\upsilon_{\rm
ion}({{\bf
x}})+J(x)-\mu+U{\scriptstyle\circ}n_{H}\right]\delta(x-x^{\prime})\;.$ (194)

That is, we now view the potential as given by $\upsilon_{\rm ion}({\bf
x})+J(x)$ and the Hartree particle density
$n_{H}(x)=i\varphi_{c}(x)-(U^{-1}{\scriptstyle\circ}J)_{x}$. This
interpretation indeed agrees with our equation (192) which can be written as

$0=\int dy\;U(x,y)\left[\left(i\varphi_{c}(y)-\int
dz\;U^{-1}(y,z)J(z)\right)+{\mathcal{G}}_{c}(y,y)\right]\;.$ (195)

That is, for a given $J(x)\neq 0$, one will solve as before the Hartree
equation but with $\upsilon_{\rm ion}({{\bf x}})\to\upsilon_{\rm ion}({{\bf
x}})+J(x)$. Once the Hartree particle density $n_{H}$ is obtained, one obtains
$i\varphi_{c}=n_{H}+U^{-1}{\scriptstyle\circ}J$. Since
$U{\scriptstyle\circ}(i\varphi_{c})$ always appears as a unit inside the
Green’s function ${\mathcal{G}}_{c}(x,x^{\prime})$, for convenience, we define

$J_{c}\equiv\,U{\scriptstyle\circ}(i\varphi_{c})=J+U{\scriptstyle\circ}n_{H}\;.$

Once $\varphi_{c}(x)$ is obtained, we shift $\phi$ by $\varphi_{c}$ and re-
express the exponent in the integrand as

$-I[\phi]-iJ{\scriptstyle\circ}\phi\Rightarrow-I[\varphi_{c}+\phi]-iJ{\scriptstyle\circ}(\varphi_{c}+\phi)\;,$

and then expand around $\varphi_{c}$. To do the expansion, we first rewrite
(20) as

$G_{\phi}^{-1}(x,x^{\prime})={\mathcal{G}}_{c}^{-1}(x,x^{\prime})+i\,b(x)\,\delta(x-x^{\prime})\equiv{\mathcal{G}}^{-1}_{c}(x,x^{\prime})+V(x,x^{\prime})\;,$
(196)

where $b=U{\scriptstyle\circ}\phi$. We may then write down

$G_{\phi}^{-1}={\mathcal{G}}^{-1}_{c}\left[{\mathbf{I}}+{\mathcal{G}}_{c}{\scriptstyle\circ}{\mathbf{V}}\right]\;,$

and

$\ln\left(G_{\phi}^{-1}\right)=\ln\left({\mathcal{G}}^{-1}_{c}\right)+\sum_{k=1}^{\infty}\frac{(-1)^{k-1}}{k}\left[{\mathcal{G}}_{c}{\scriptstyle\circ}{\mathbf{V}}\right]^{k}\;.$
(197)

Note that

$\left[{\mathcal{G}}_{c}{\scriptstyle\circ}{\mathbf{V}}\right]_{x,z}=\int\\!\\!dy\,{\mathcal{G}}_{c}(x,y)V(y,z)=\int
dy{\mathcal{G}}_{c}(x,y)\delta(y-z)\left(i\,b(y)\,\right)={\mathcal{G}}_{c}(x,z)\left(ib(z)\right)\;.$

Consequently,

$\displaystyle\text{Tr}\ln\left(G_{\phi}^{-1}\right)$ $\displaystyle=$
$\displaystyle\text{Tr}\ln\left({\mathcal{G}}^{-1}_{c}\right)+\int\\!dx_{1}\,{\mathcal{G}}_{c}(x_{1},x_{1})(ib(x_{1}))$
(198)
$\displaystyle-\frac{1}{2}\int\\!dx_{1}dx_{2}\,{\mathcal{G}}_{c}(x_{1},x_{2}){\mathcal{G}}_{c}(x_{2},x_{1})(ib(x_{1}))(ib(x_{2}))$
$\displaystyle+\sum_{k=3}^{\infty}\frac{(-1)^{k-1}}{k}\int\\!dx_{1}\ldots
dx_{k}\;{\mathcal{G}}_{c}(x_{k},x_{1})\ldots{\mathcal{G}}_{c}(x_{k-1},x_{k})(ib(x_{1}))\ldots(ib(x_{k}))\;.$

Note that in the final expression of the exponent of the integrand in (23) the
terms linear in $\phi$ (or $b$) cancel out as one may verify and we arrive at

$\displaystyle-I[\varphi_{c}+\phi]-iJ{\scriptstyle\circ}(\varphi_{c}+\phi)$
$\displaystyle=$
$\displaystyle-\frac{1}{2}\varphi_{c}{\scriptstyle\circ}\,U{\scriptstyle\circ}\varphi_{c}-\varphi_{c}{\scriptstyle\circ}b-\frac{1}{2}b{\scriptstyle\circ}\,U^{-1}{\scriptstyle\circ}b+\text{Tr}\ln(G_{\phi}^{-1})$
(199)
$\displaystyle+\frac{1}{2}\text{Tr}\ln(U)-iJ{\scriptstyle\circ}\varphi_{c}-iJ{\scriptstyle\circ}\phi$
$\displaystyle=$
$\displaystyle\frac{1}{2}\text{Tr}\ln(U)-\frac{1}{2}\varphi_{c}{\scriptstyle\circ}\,U{\scriptstyle\circ}\varphi_{c}+\text{Tr}\ln({\mathcal{G}}_{c}^{-1})-iJ{\scriptstyle\circ}\varphi_{c}$
$\displaystyle-\frac{1}{2}b{\scriptstyle\circ}\tilde{\cal
D}_{c}^{-1}{\scriptstyle\circ}b+\sum_{k=3}^{\infty}I^{(k)}[\varphi_{c}]{\scriptstyle\circ}\,b_{1}\ldots{\scriptstyle\circ}\,b_{k}\;,$

where

$\tilde{\cal D}_{c}^{-1}=U^{-1}-D_{c}\;,$
$D_{c}(x,y)={\mathcal{G}}_{c}(x,y)\,{\mathcal{G}}_{c}(y,x)\;,$

and

$I^{(k)}[\varphi_{c}]{\scriptstyle\circ}\,b_{1}\ldots{\scriptstyle\circ}\,b_{k}\equiv\frac{(-1)^{k-1}}{k}\int\\!dx_{1}\ldots
dx_{k}{\mathcal{G}}_{c}(x_{k},x_{1})\ldots{\mathcal{G}}_{c}(x_{k-1},x_{k})\left[(ib(x_{1}))\ldots(ib(x_{k}))\right]\;.$
(200)

We therefore have (based on the linked-cluster theorem)

$\displaystyle\beta W_{\phi}[J]$ $\displaystyle=$
$\displaystyle\frac{1}{2}\varphi_{c}{\scriptstyle\circ}\,U{\scriptstyle\circ}\varphi_{c}-\text{Tr}\ln\left({\mathcal{G}}_{c}^{-1}\right)+iJ{\scriptstyle\circ}\varphi_{c}$
(201) $\displaystyle+\frac{1}{2}\text{Tr}\ln\left(\tilde{\cal
D}_{c}^{-1}{\scriptstyle\circ}\,U\right)-\sum_{n=1}^{\infty}\frac{1}{n!}\langle\left[\sum_{k=3}^{\infty}I^{(k)}[\varphi_{c}]{\scriptstyle\circ}\,b_{1}\ldots{\scriptstyle\circ}\,b_{k}\right]^{n}\rangle_{\rm
conn.}\;,$

where the subscript “${\rm conn.}$” is used to indicate connected Feynman
diagrams.

Consequently, the effective action, defined by

$\Gamma[n]=\beta
W_{\phi}[J]-\frac{1}{2}J{\scriptstyle\circ}\,U^{-1}{\scriptstyle\circ}J-J{\scriptstyle\circ}n=\beta
W_{\phi}[J]+\frac{1}{2}J{\scriptstyle\circ}\,U^{-1}{\scriptstyle\circ}J-iJ{\scriptstyle\circ}\varphi\;$

and with
$\varphi=\varphi_{c}+\langle\phi\rangle\equiv\varphi_{c}+\tilde{\varphi}$ as
well as with $J=J_{c}-U{\scriptstyle\circ}n_{H}$, can now be written as

$\displaystyle\Gamma[n]$ $\displaystyle=$
$\displaystyle-\text{Tr}\ln\left({\mathcal{G}}_{c}^{-1}\right)-J_{c}{\scriptstyle\circ}n_{H}+\frac{1}{2}n_{H}{\scriptstyle\circ}\,U{\scriptstyle\circ}n_{H}+\frac{1}{2}\text{Tr}\ln\left(\tilde{\cal
D}_{c}^{-1}{\scriptstyle\circ}\,U\right)$ (202) $\displaystyle-
iJ{\scriptstyle\circ}{\tilde{\varphi}}-\sum_{n=1}^{\infty}\frac{1}{n!}\langle\left[\sum_{k=3}^{\infty}I^{(k)}[\varphi_{c}]{\scriptstyle\circ}\,b_{1}\ldots{\scriptstyle\circ}\,b_{k}\right]^{n}\rangle_{\rm
conn.}\;.$

This expression should be contrasted with Eq. (45) where the true particle
density is introduced as the natural variable.

From the definition of $\varphi$ (see Eq. (26) we have

$i\varphi(x)=\frac{\int
D\phi\left(i\phi(x)\right)e^{-I[\phi]-iJ{\scriptstyle\circ}\phi}}{\int
D\phi\,e^{-I[\phi]-iJ{\scriptstyle\circ}\phi}}=i\varphi_{c}(x)+\int
dy\,U^{-1}(x,y)\frac{\int
Db\left(ib(y)\right)e^{-I[\varphi_{c}+\phi]-iJ{\scriptstyle\circ}(\varphi_{c}+\phi)}}{\int
Db\,e^{-I[\varphi_{c}+\phi]-iJ{\scriptstyle\circ}(\varphi_{c}+\phi)}}\;,$
(203)

which means that

${i\tilde{\varphi}}=\frac{\int
D\,b\left(iU^{-1}{\scriptstyle\circ}b\right)e^{-\frac{1}{2}b{\scriptstyle\circ}\tilde{\cal
D}_{c}^{-1}{\scriptstyle\circ}b+\sum_{k=3}^{\infty}\frac{1}{k!}I^{(k)}[\varphi_{c}]{\scriptstyle\circ}\,b_{1}\ldots{\scriptstyle\circ}\,b_{k}}}{\int
Db\,e^{-\frac{1}{2}b{\scriptstyle\circ}\tilde{\cal
D}_{c}^{-1}{\scriptstyle\circ}b+\sum_{k=3}^{\infty}\frac{1}{k!}I^{(k)}[\varphi_{c}]{\scriptstyle\circ}\,b_{1}\ldots{\scriptstyle\circ}\,b_{k}}}\;.$
(204)

On the basis of a simple replica argument Negele and Orland (1988), one sees
that the above expression may be written as

$\displaystyle{i\tilde{\varphi}}(x)$ $\displaystyle=$
$\displaystyle\sum_{n=0}^{\infty}\frac{1}{n!}\int
dy\,U^{-1}(x,y)\langle\left[\sum_{k=3}^{\infty}I^{(k)}[\varphi_{c}]{\scriptstyle\circ}\,b_{1}\ldots{\scriptstyle\circ}\,b_{k}\right]^{n}ib(y)\rangle_{\rm
conn.}$ (205) $\displaystyle=$
$\displaystyle\sum_{n=1}^{\infty}\frac{1}{n!}\int
dy\,U^{-1}(x,y)\langle\left[\sum_{k=3}^{\infty}I^{(k)}[\varphi_{c}]{\scriptstyle\circ}\,b_{1}\ldots{\scriptstyle\circ}\,b_{k}\right]^{n}ib(y)\rangle_{\rm
conn.}\;,$

where the $n=0$ term vanishes because $\int
Db\,e^{-\frac{1}{2}b{\scriptstyle\circ}\tilde{\cal
D}_{c}^{-1}{\scriptstyle\circ}b}ib(x)=0$. Eq. (205) implies that
$i\tilde{\varphi}(x)$ at each point $x$ is a functional of $\varphi_{c}$.
Returning to (202), this implies that one may express the effective action
$\Gamma$ in terms of $i\varphi_{c}$ although its canonical argument is
supposed to be $n$. The relation between $n$ and $i\varphi_{c}$ is obtained
through

$n=i\varphi_{c}+i\tilde{\varphi}[i\varphi_{c}]-U^{-1}{\scriptstyle\circ}J[i\varphi_{c}]=n_{H}+i\tilde{\varphi}[i\varphi_{c}]$
(206)

with $i\tilde{\varphi}(x)$ expressed as a functional of $i\varphi_{c}$. This
implies that given an $i\varphi_{c}$, we may obtain $n$ through (i) the
diagrammatic expansion of (205) to produce the corresponding $i\varphi$, and
(ii) equation (192) to find $-U^{-1}{\scriptstyle\circ}J$. Since in general,
we have

$\int d{\bf x}\,\left(i\tilde{\varphi}({\bf x})\right)\neq 0\;,$

one cannot expect the integral of the Hartree electron density to be $N_{e}$,
the total number of electrons. Instead, we have

$N_{e}=\int d{\bf x}\,n({\bf x})=\int d{\bf x}\,\left[n_{H}({\bf
x})+i\tilde{\varphi}({\bf x})\right]$ (207)

and $\int d{\bf x}\,n_{H}({\bf x})$ is not necessarily an integer. In
particular, when $\int d{\bf x}\,n_{H}({\bf x})$ deviates significantly from
$N_{e}$ or $\int d{\bf x}\,|i\tilde{\varphi}({\bf x})|\gg 1$, forcing $\int
d{\bf x}\,n_{H}({\bf x})=N_{e}$ may lead to the occurrence of a self-
consistent solution that differs significantly from the true solution.

Note that $n_{H}({\bf x})=-{\mathcal{G}}_{c}(x,x)$ and in the absence of the
source term $n_{H}({\bf x})=i\varphi_{c}({\bf x})$. Furthermore, since
$(U^{-1}{\scriptstyle\circ}J)_{x}\propto\nabla_{{\bf x}}^{2}J(x)$, $\int
n_{H}({\bf x})\,d{\bf x}=\int(i\varphi_{c}({\bf x}))\,d{\bf x}$. The
constraint (207) for the total number of electrons can also be written as

$N_{e}=\int d{\bf x}\,n({\bf x})=\int d{\bf x}\,\left[n_{H}({\bf
x})+i\tilde{\varphi}({\bf x})\right]=\int d{\bf x}\,\left[i\varphi_{c}({\bf
x})+i\tilde{\varphi}({\bf x})\right]\;.$ (208)

The Hartree-like Green’s function ${\mathcal{G}}_{c}(x,x^{\prime})$ shown in
(193) may be viewed as a functional of $i\varphi_{c}(x)$. When expressing
${\mathcal{G}}_{c}(x,x^{\prime})$ using only single particle orbitals, we
define in Eq. (88) $v({\bf x})\equiv\upsilon_{\rm ion}({\bf
x})-\mu+U{\cdot}(i\varphi_{c})_{{\bf x}}$. For the evaluation of
${\mathcal{G}}_{c}(x,y)$, we solve first the eigensystem (93). The single-
particle wave functions (91) associated with $\hat{h}$ are to be obtained
self-consistently. Basically, one starts with a guess for the electronic
density $i\varphi_{c}({\bf x})$ satisfying $\int d{\bf x}(i\varphi_{c}({\bf
x}))\approx N_{e}$, where $N_{e}$ is the number of electrons. One then obtains
the single-particle wave functions, and then computes the corresponding
Green’s function ${\mathcal{G}}_{c}$, obtains $\int d{\bf
x}(i\tilde{\varphi}({\bf x}))$ and tunes the chemical potential to ensure that
$-\int d{\bf x}\;{\mathcal{G}}_{c}({\bf x},{\bf x})=N_{e}-\int d{\bf
x}(i\tilde{\varphi}({\bf x}))$. One then takes $-{\mathcal{G}}_{c}({\bf
x},{\bf x})$ in place of $i\varphi_{c}({\bf x})$ in the next round of
iteration until convergence is reached.

The procedure of the saddle-point method is now obvious. One starts with an
external potential, and then determines the Hartree electron density via
$n_{H}({\bf x})=-{\mathcal{G}}_{c}(x,x)$,
$n_{H}=i\varphi_{c}-U^{-1}{\scriptstyle\circ}J$ and Eq. (208). This self-
consistent procedure will also provide the physical electron density $n({\bf
x})$. The ground state energy is obtained by using (202) to calculate the
effective action, which is the ground state energy times $\beta$ in the limit
$T\to 0$. When one wishes to obtain the density functional $\Gamma[n]$ at an
electron density other than $n_{T}$, one adds the source term into the
potential $\upsilon_{\rm ion}({\bf x})$ and then solves for the Hartree
density, shown in (194), as outlined above.

###  VI.2 Remark on the single-particle limit

The single electron limit for the Hartree-method is more complicated than for
the method presented in section IV.3. Because $\int n_{H}({\bf x})\,d{\bf x}$
is not necessarily $1$ in the single electron limit, due to Eq. (208), in
general one needs to obtain $n_{1}$ self-consistently. The other issue is that
the diagrammatic expansion contained in (202) does not cover all the diagrams
for a given order of $U$. This means that one can’t use the Hugenholtz diagram
argument to eliminate vertex matrix element even when $n_{1}\leq 1$. To see
this point explicitly, we rewrite (202) as

$\displaystyle\Gamma[n]$ $\displaystyle=$
$\displaystyle-\text{Tr}\ln\left({\mathcal{G}}_{c}^{-1}\right)+(J-J_{c}){\scriptstyle\circ}n_{H}+\frac{1}{2}n_{H}{\scriptstyle\circ}\,U{\scriptstyle\circ}n_{H}+\frac{1}{2}\text{Tr}\ln\left(\tilde{\cal
D}_{c}^{-1}{\scriptstyle\circ}\,U\right)$ (209)
$\displaystyle-\sum_{n=1}^{\infty}\frac{1}{n!}\langle\left[\sum_{k=3}^{\infty}I^{(k)}[\varphi_{c}]{\scriptstyle\circ}\,b_{1}\ldots{\scriptstyle\circ}\,b_{k}\right]^{n}\rangle_{\rm
conn.}-J{\scriptstyle\circ}n\;.$

Evidently, except for the last term above, the rest of the terms must
constitute $\beta W[J]$. Since $J-J_{c}=-U{\scriptstyle\circ}n_{H}$, one can
also write $\beta W[J]$ as

$\displaystyle\beta W[J]$ $\displaystyle=$
$\displaystyle-\text{Tr}\ln\left({\mathcal{G}}_{c}^{-1}\right)-n_{H}{\scriptstyle\circ}\,U{\scriptstyle\circ}n_{H}+\frac{1}{2}n_{H}{\scriptstyle\circ}u{\scriptstyle\circ}n_{H}+\frac{1}{2}\text{Tr}\ln\left(\tilde{\cal
D}_{c}^{-1}{\scriptstyle\circ}\,U\right)$
$\displaystyle-\sum_{n=1}^{\infty}\frac{1}{n!}\langle\left[\sum_{k=3}^{\infty}I^{(k)}[\varphi_{c}]{\scriptstyle\circ}\,b_{1}\ldots{\scriptstyle\circ}\,b_{k}\right]^{n}\rangle_{\rm
conn.}\;\;.$

Diagrammatic expansion of the last two terms shows that the following diagrams

\begin{picture}(96.0,20.0)(0.0,-3.0) \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \end{picture}
\begin{picture}(60.0,20.0)(0.0,-3.0) \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){}
\put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \end{picture}

of order $U^{2}$ are absent, when compared to the regular field theoretic
perturbation calculation. This is not a disadvantage of the method. Instead,
what our derivation shows is that the missing diagrams eventually will be
compensated by the $-n_{H}{\scriptstyle\circ}\,U{\scriptstyle\circ}n_{H}$
term. However, it is obvious that the saddle-point formalism makes the single-
electron limit hard to analyze.

When $n_{1}>0$ but $n_{1}\ll 1$ at the low temperature limit, we know that
$D_{c}(x,y)={\mathcal{G}}_{c}(x,y){\mathcal{G}}_{c}(y,x)$ will be of order
$n_{1}$. This is because

${\mathcal{G}}_{c}(x,y)=\sum_{\alpha}\phi_{\alpha}({\bf
x})\phi_{\alpha}^{*}({\bf
y})e^{-(\varepsilon_{\alpha}-\mu)(\tau_{x}-\tau_{y})}\times\left\\{\begin{array}[]{l
r}(-n_{\alpha})&{\rm if\leavevmode\nobreak\ }\tau_{x}\leq\tau_{y}\\\
(1-n_{\alpha})&{\rm if\leavevmode\nobreak\
}\tau_{x}>\tau_{y}\end{array}\right.\;,$

and whenever $\tau_{x}\leq\tau_{y}$, the propagator is of order $n_{1}$. Since
$D_{c}(x,y)={\mathcal{G}}_{c}(x,y){\mathcal{G}}_{c}(y,x)$, one of the
propagators in the product must be of order $n_{1}$. In principle, one needs
to solve for the occupation number $n_{1}$ of the lowest energy state using
$i\varphi_{c}({\bf
x})-\left(U^{-1}{\scriptstyle\circ}J\right)_{x}=-{\mathcal{G}}_{c}(x,x)$ and
Eq. (208). Nevertheless, the correct one particle limit can be seen if one
starts with a chemical potential $\mu$ such that $n_{1}\approx 0$. In this
case, we have at the $J=0$ limit $i\varphi_{c}(x)=n_{H}({\bf x})\approx 0$ as
well as $D_{c}(x,y)\propto n_{1}\approx 0$. This way, the higher order
exchange-correlation terms may be viewed as the having smaller contributions
and one may control the accuracy by controlling the number of higher order
terms included. Of course, one then has $n({\bf x})\approx
i\tilde{\varphi}({\bf x})$ and the condition $\int d{\bf
x}\left(i\tilde{\varphi}({\bf x})\right)=1-n_{1}$ must be satisfied.

###  VI.3 Obtaining excitations using the Hartree method

In general the excitations are determined by Eqs. (176) and (177). Under the
Hartree formalism described in this section, the natural variable used is the
Hartree density $n_{H}$ rather than the true particle density $n_{T}$. One
thus must transform the variable used in Eqs. (176) and (177) from $n_{T}$ to
$n_{H}$. We describe below how this can be achieved.

At the physical condition ($J=0$), one has

$0=\left.\frac{\delta\Gamma}{\delta n(x)}\right|_{n=n_{T}}=\left.\int
dx_{1}dx_{2}\frac{\delta\Gamma}{\delta
n_{H}(x_{2})}\right|_{n=n_{T}}\left.\frac{\delta n_{H}(x_{2})}{\delta
J_{c}(x_{1})}\right|_{n=n_{T}}\left.\frac{\delta J_{c}(x_{1})}{\delta
n(x)}\right|_{n=n_{T}}\;.$

Note that $\delta n_{H}/\delta J_{c}$ contains no zero mode because
$n_{H}=\delta W_{H}/\delta J_{c}$, where
$W_{H}[J_{c}]\equiv-\text{Tr}\ln\left({{\mathcal{G}}_{c}}^{-1}\right)$, and
$\delta^{2}W_{H}/\delta J_{c}\delta J_{c}$ is known to be strictly negative
from a theorem proved by Valiev and Fernando. Valiev and Fernando (1997) The
strict negative-definiteness of $\delta^{2}W/\delta J\delta J$ can also be
interpreted as the stability condition

$\int dx\,\delta n(x)\,\delta J(x)<0\;,$

which means raising the local one-particle potential leads to an average
decrease of the local particle concentration and vice versa. The strict
negative-definiteness means that $\delta n_{H}/\delta J_{c}$ contains no zero
modes and is invertible. Also, because $n=n_{H}+i\tilde{\varphi}_{c}$, $\delta
n/\delta J_{c}=\delta n_{H}/\delta J_{c}+\delta(i\tilde{\varphi})/\delta
J_{c}$ exists via diagrammatic expansion of $i\tilde{\varphi}$ in terms of
$J_{c}$. The existence of $\delta n/\delta J_{c}$ implies that $\delta
J_{c}/\delta n$ is invertible (i.e., has no zero eigenvalues). Therefore,
after multiplying the inverse of $\delta J_{c}/\delta n$ and $\delta
n_{H}/\delta J_{c}$ on both sides of the above equation, one has

$\left.\frac{\delta\Gamma}{\delta n_{H}(x)}\right|_{n=n_{T}}=0\;,$

which then leads to

$\displaystyle\Gamma^{(2)}(x,y)$ $\displaystyle=$
$\displaystyle\left.\frac{\delta^{2}\Gamma}{\delta n(x)\delta
n(y)}\right|_{n=n_{T}}$ $\displaystyle=$ $\displaystyle\left[\int
dx_{1}dx_{2}dy_{1}dy_{2}\frac{\delta J_{c}(x_{1})}{\delta n(x)}\frac{\delta
n_{H}(x_{2})}{\delta J_{c}(x_{1})}\frac{\delta^{2}\Gamma}{\delta
n_{H}(x_{2})\delta n_{H}(y_{2})}\frac{\delta n_{H}(y_{2})}{\delta
J_{c}(y_{1})}\frac{\delta J_{c}(y_{1})}{\delta n(y)}\right]_{n=n_{T}}\;.$

Using Eq. (172), one may write $\Gamma^{(2)}({\bf x},{\bf y},\omega)$ as the
analytic continuation of the integral of $\Gamma^{(2)}(x,y)$ over time. To
achieve this goal, one first writes down

$\displaystyle\Gamma^{(2)}({\bf x},{\bf y},i\nu_{n})$ $\displaystyle=$
$\displaystyle\int_{0}^{\beta}d(\tau_{y}-\tau_{x})\,e^{i\nu_{n}(\tau_{y}-\tau_{x})}\,\Gamma^{(2)}(x,y)$
(210) $\displaystyle=$
$\displaystyle\int_{0}^{\beta}d(\tau_{y}-\tau_{x})\,e^{i\nu_{n}(\tau_{y}-\tau_{y_{1}}+\tau_{y_{1}}-\tau_{y_{2}}+\tau_{y_{2}}-\tau_{x_{2}}+\tau_{x_{2}}-\tau_{x_{1}}+\tau_{x_{1}}-\tau_{x})}\,\Gamma^{(2)}(x,y)$
$\displaystyle\equiv$ $\displaystyle\int d{\bf x}_{1}d{\bf x}_{2}d{\bf
y}_{1}d{\bf y}_{2}f^{{\bf x}_{1}}_{{\bf x}}\\!(-i\nu_{n})g^{{\bf x}_{2}}_{{\bf
x}_{1}}\\!(-i\nu_{n})\tilde{\Gamma}^{(2)}({\bf x}_{2},{\bf
y}_{2},i\nu_{n})g^{{\bf y}_{2}}_{{\bf y}_{1}}\\!(i\nu_{n})f^{{\bf
y}_{1}}_{{\bf y}}\\!(i\nu_{n})$ $\displaystyle=$ $\displaystyle\int d{\bf
x}_{2}d{\bf y}_{2}h^{{\bf x}_{2}}_{{\bf
x}}(-i\nu_{n})\tilde{\Gamma}^{(2)}({\bf x}_{2},{\bf y}_{2},i\nu_{n})h^{{\bf
y}_{2}}_{{\bf y}}(i\nu_{n})\;,$

where

$\displaystyle\tilde{\Gamma}^{(2)}({\bf x},{\bf y},i\nu_{n})$ $\displaystyle=$
$\displaystyle\int_{0}^{\beta}d(\tau_{y}-\tau_{x})\frac{\delta^{2}\Gamma}{\delta
n_{H}(x)\delta n_{H}(y)}$ $\displaystyle g^{{\bf y}_{2}}_{{\bf
y}_{1}}\\!(i\nu_{n})$ $\displaystyle\equiv$
$\displaystyle\left.\int_{0}^{\beta}d(\tau_{y_{1}}-\tau_{y_{2}})\,e^{i\nu_{n}(\tau_{y_{1}}-\tau_{y_{2}})}\,\frac{\delta
n_{H}(y_{2})}{\delta J_{c}(y_{1})}\right|_{n=n_{T}}$ $\displaystyle f^{{\bf
y}_{1}}_{{\bf y}}\\!(i\nu_{n})$ $\displaystyle\equiv$
$\displaystyle\left.\int_{0}^{\beta}d(\tau_{y}-\tau_{y_{1}})\,e^{i\nu_{n}(\tau_{y}-\tau_{y_{1}})}\,\frac{\delta
J_{c}(y_{1})}{\delta n(y)}\right|_{n=n_{T}}$ $\displaystyle h^{{\bf
y}_{2}}_{{\bf y}}(i\nu_{n})$ $\displaystyle\equiv$ $\displaystyle\int d{\bf
y}_{1}\;g^{{\bf y}_{2}}_{{\bf y}_{1}}\\!(i\nu_{n})f^{{\bf y}_{1}}_{{\bf
y}}\\!(i\nu_{n})\;.$

Therefore, Eq. (179) for finding the excitations becomes

$\int\\!\\!d{\bf z}\,\tilde{\Gamma}^{(2)}({\bf x},{\bf
z},-\omega)\,\tilde{\Delta}({\bf z})=0\;,$ (211)

with $\tilde{\Delta}({\bf z})$ given by

$\tilde{\Delta}({\bf z})=\int d{\bf y}\,h^{{\bf z}}_{{\bf
y}}(-\omega)\Delta({\bf y})\;.$

One therefore obtains the excitation energy via solving Eq. (211). Since
$h^{{\bf z}}_{{\bf y}}(-\omega)$ is invertible, one may also obtain
$\Delta({\bf z})$ via $\tilde{\Delta}({\bf z})$ if desired. Within this
framework, the first two terms of (202) correspond to the non-interacting part
$\Gamma_{0}$ of the effective action in section V. The protocol for obtaining
the excitation energy is then the same as described in section V with the one
exception that the ground state of the Hartree system is not the same as the
ground state of the real system.

##  VII Discussion and Future Directions

In this paper we focus on the auxiliary field method applied to the
development of the density functional. It is natural to inquire into the
physical meaning of the auxiliary, bosonic field $\phi$ introduced in eqs
(15-20). In the path integral treatment of relativistic quantum
electrodynamics, if one were to integrate out the photon field, one generates
the current-current interaction which is quartic in fermionic field. When
viewing this process backwards, one sees that the quartic fermionic
interaction is disentangled by introducing the photon field. The $\phi$ field
here thus plays a similar role to the photon field as it disentangles the
quartic fermionic interaction term. As we will argue below, the $\phi$
variable is closely related to the “time” component of the photon field, i.e.,
the electric potential.

To make the connection between the photon field and $\phi$, let us first seek
the nonrelativistic limit of the Lagrangian density,
${\mathcal{L}}_{EM}=-F^{\mu\nu}F_{\mu\nu}/(16\pi)$ , associated with the
relativistic photon field. Note that in the limit when the speed of light
approaches infinity (removal of terms involving time derivative of the three-
vector potential), the Lagrangian density turns into $(\nabla
A_{0})^{2}/8\pi$, which contains no quadratic derivative with respect to time.
Therefore, at finite temperature (with $it\to\tau$) this implies that the
exponent associated with the photon field path integral will behave as

$i\int dt\,d{\bf x}\;{\mathcal{L}}_{EM}\to\int\\!d\tau\\!\int\\!d{\bf
x}\frac{(\nabla A_{0})^{2}}{8\pi}=-\int\\!d\tau\\!\int\\!d{\bf
x}\frac{1}{8\pi}\,A_{0}\nabla^{2}A_{0}\;.$ (212)

Setting $U({\bf x}-{\bf y})=e^{2}/|{\bf x}-{\bf y}|$ [thus
$U^{-1}(x,y)=-\frac{1}{4\pi e^{2}}\nabla^{2}_{{\bf x}}\delta({\bf x}-{\bf
y})$] and comparing (212) with Eqs. (15-16), we make the identification
$(iU{\scriptstyle\circ}\phi)/e=A_{0}$, because now

$-\frac{1}{2}\phi{\scriptstyle\circ}\,U{\scriptstyle\circ}\phi\to\frac{1}{2}A_{0}\,{\scriptstyle\circ}(e^{2}u^{-1})\,{\scriptstyle\circ}A_{0}=-\int
d\tau\int d{\bf x}\frac{1}{8\pi}\,A_{0}\nabla^{2}A_{0}\;.$