Editors' Suggestion

Emergence of an excitonic collective mode in the dilute electron gas
Yasutami Takada
Phys. Rev. B 94, 245106 (2016) – Published 2 December 2016

In the early 1950s, Bohm and Pines (BP) developed the random-phase approximation(RPA) theory and revealed the plasmon as an emergent collective mode in metals. In
this paper, for the first time ever since the original BP theory, a novel collective modeis predicted to emerge in the dilute electron liquid. The electronic compressibility ofthis mode is negative and thus dielectric catastrophe appears, predicted by a theory
that goes far beyond RPA and fully satisfies various sum rules, constraints, andasymptotic behavior. While the plasmon occurs by the coherent excitations ofmutually uncorrelated electrons and holes, the new mode is associated with theexcitations of a large collection of excitons (or tightly bound electron-hole pairs).
Mathematically, this excitonic collective mode is intimately connected with a newly
discovered singularity in the dielectric function in the electron liquid. This collectivemode should add a new dimension to plasmonics in nanotechnology in the future.
（http://journals.aps.org/prb/highlights）

PHYSICAL REVIEW B 94, 245106 (2016)

Emergence of an excitonic collective mode in the dilute electron gas
Yasutami Takada
Institute for Solid State Physics, University of Tokyo, Kashiwa, Chiba 277-8581, Japan(Received 7 September 2016; revised manuscript received 15 November 2016; published 2 December 2016)
By comparing two expressions for the polarization function (q,iω) given in terms of two different local-field
factors, G+ (q,iω) and Gs (q,iω), we have derived the kinetic-energy-fluctuation (or sixth-power) sum rule forthe momentum distribution function n( p) in the three-dimensional electron gas. With use of this sum rule,
together with the total-number (or second-power) and the kinetic-energy (or fourth-power) sum rules, we haveobtained n( p) in the low-density electron gas at negative compressibility (namely, rs > 5.25 with rs being theconventional density parameter) up to rs ≈ 22 by improving on the interpolation scheme due to Gori-Giorge andZiesche proposed in 2002. The obtained results for n( p) combined with the improved form for Gs (q,ω + i0+ )are employed to calculate the dynamical structure factor S(q,ω) to reveal that a giant peak, even bigger than theplasmon peak, originating from an excitonic collective mode made of electron-hole pair excitations, emerges in
the low-ω region at |q| near 2pF (pF : the Fermi wave number). Connected with this mode, we have discovered asingular point in the retarded dielectric function at ω = 0 and |q| ≈ 2pF .
DOI: 10.1103/PhysRevB.94.245106
I. INTRODUCTION

The electron gas, an assembly of N electrons embedded in a
uniform positive rigid background, has been investigated for avery long time, not only because it is relevant in understanding
the exchange-correlation effects in simple metals and semiconductors, but also because it provides indispensable information
for the actual implementation of the density functional theory(DFT), especially in its local density approximation (LDA) [1].
From a physical point of view, we may claim that thissystem is even more important in the low-density regime(namely, pF a0  1 with pF the Fermi wave number anda0 the lattice constant), because in this regime the firstprinciples Hamiltonian to describe conduction electrons ina metal is universally reduced to the one for the electrongas, leading us to recognize that understanding the physicalproperties of the dilute electron gas is nothing but revealing
the universal electron-correlation behavior of actual dilutemetals.
The dilute electron gas, however, is not well understood; thethermodynamic quantities, such as the correlation energy εcand the compressibility κ, are accurately obtained as a functionof the density parameter rs by the Green’s function Monte
Carlo (GFMC) method [2] and the interpolation formulas toreproduce the GFMC data [3,4], where rs is defined by rs =(4π n/3)−1/3 /aB = (αpF aB )−1 with n = N/ t the electron
density, t the total volume of the system, α = (4/9π )1/3 ≈0.5211, and aB the Bohr radius. The diffusion Monte Carlo(DMC) calculations are done to analyze the ground-statephases, including spin polarization, in the wide range of rs [5]and recently for 0.5  rs  20 in more detail [6], but we do
not know the precise behavior of other physical quantities,
such as the momentum distribution function n( p) for rs > 5.
It is true that some quantum Monte Carlo calculations aredone to obtain n( p) for rs  10 [7,8], but the results are notvery accurate due probably to large size effects and improper
starting trial functions, judging from the assessment of theiraccuracy by the sum rules for n( p) [9,10]. As for the dynamical
structure factor S(q,ω), no reliable data are available in the
dilute metallic regime.
2469-9950/2016/94(24)/245106(14)One of the characteristic features of the dilute electron gas is
the appearance of “dielectric catastrophe” associated with thedivergence of κ at rs = rsc ≡ 5.25 followed by negative κ forrs > rsc [11,12], implying that the static polarization function
(q,0) in the long-wavelength limit is also negative (andthus dielectrically anomalous) owing to the compressibility
sum rule, lim|q|→0 (q,0) = n2 κt [13]. This anomaly for
rs > rsc brings about the curious behavior of the ion-ion paircorrelation function in expanded liquid alkali metals [14–16],
but it does not induce any instabilities in the electron gas itself,
as long as 1/κ changes continuously from positive to negative
through the point of 1/κ = 0 [17], because the electronelectron effective interaction is not determined only by thedielectric function ε(q,iω) = 1 + V (q)(q,iω) with V (q) =4π e2 /(t q 2 ) the bare Coulomb interaction [18], excluding theoccurrence of CDW-type instabilities at rs = rsc [19].
By analyzing εR (q,ω) [= ε(q,ω + i0+ )], the retarded
dielectric function as rs approaches rsc from the positive sideof κ (or rs < rsc ), we identify the physical origin of thisdivergence of κ as the “enhanced excitonic effect” or the effectof strong electron-hole attraction in an electron-hole singlepair excitation [20]. For rs > rsc , by considering causality,
Takayanagi and Lipparini [21] conclude that (q,iω) contains
an anomalous part a (q,iω) made of a single pair of poles atiω = ±i ω̃ex (q) in the form ofa (q,iω) =fex (q)
1 ne2,
V (q) m ω2 − ω̃ex (q)2

(1)

with m the mass of a free electron. These poles at which(q,iω) diverges, however, do not give any prominent
structures in the observable physical quantity S(q,ω), which
can be calculated by S(q,ω) = −ImQRc (q,ω)/π for ω > 0 atzero temperature (T = 0) with the retarded charge response
function QRc (q,ω), given by
QRc (q,ω) = −(q,ω + i0+ ).
1 + V (q)(q,ω + i0+ )(2)

Then a couple of questions arise: what is the physical roleof these poles and how the strong excitonic effect manifests

245106-1©2016 American Physical Society

PHYSICAL REVIEW B 94, 245106 (2016),(7)with Gs (q,iω) the local-field factor as introduced by Richardson and Ashcroft [28]. This factor is different from theconventional local-field factor G+ (q,iω) [18], with which(q,iω) is expressed in another way as0 (q,iω),1 − V (q)G+ (q,iω)0 (q,iω)

(8)

-1

1

-1
S(q,ω)/N [units of (27.21eV) ]uiω − ε p+q + ε p

15
10
5
]
E
of
ts
F
ni
ω [u

1

2

q[

pσ

(q,iω) =p n( p + q) − n( p)

2

0

20

3

F ]

where WI (q,iω) is defined as

5

of

(6)

plasmonnitsWI (q,iω),
1 − V (q)Gs (q,iω)WI (q,iω)0[u

We shall derive this sum rule in Sec. II with giving a concretevalue for the right-hand side of Eq. (5).
Second, once n( p) is known, we can calculate (q,iω) by
the formula given in GW̃ WI as

WI (q,iω) =

10

0

5

15
10
fE ]
ω [units o F

20

FIG. 1. Bird’s-eye view of S(q,ω) in the electron gas at zerotemperature at (a) a usual metallic density rs = 4, in which theplasmon peak dominates, and (b) a very dilute density rs = 22, in
which an excitonic collective peak, much bigger than the plasmon
peak, emerges.

with 0 (q,iω) the Lindhard function, given by0 (q,iω) =

(4)

pσ

plasmon

20
3

15



with σ the spin variable, ε p = p2 /2m, and EF = pF2 /2m theFermi energy. This GZ scheme can give rather accurate n( p)for rs  12. Here we shall improve on it by increasing thenumber of input data for n( p) that we obtain up to rs = 10 in
GW̃ WI . We also use, in addition to the above two, the third
sum rule for the total kinetic-energy fluctuation, written as
ε2p n( p) = N [(ε p − ε p )2  + ε p 2 ].
(5)

(q,iω) =

20

40

q


∂εc
3
,
EF − εc − rs
ε p n( p) = N ε p  = N
5
∂rs
pσ



25

p

pσ

(3)

60

30

0

excitonic collective

80(a) rs=4, T=0nitn( p) = N,(b) rs =22, T=0

100

F ]



120

so
f

itself in S(q,ω) for rs > rsc ? We shall address these questions
by accurately determining ω̃ex (q), “the oscillator strength” ofthe pole fex (q), and S(q,ω) for rs from rsc up to about 20 inwhich the ground state has already been confirmed to be aparamagnetic metal [6].
A straightforward way of obtaining accurate results for
S(q,ω) is to perform the highly self-consistent calculation in
the GW scheme [22] which includes the vertex function
satisfying the Ward identity (WI) [23,24], just as done forS(q,ω) for rs  5 [25], but it turns out that this GW does
not work well in the dielectric-catastrophe regime. Thus it isrevised into the GW̃ WI scheme [10,26] to obtain the selfconsistent results for S(q,ω) as well as n( p) up to rs = 10,
but there still remains the problem of reaching a fully selfconsistent solution for rs > 10.
Confronted with this difficulty, we will take the following
strategy: First, we reconsider the parametrization scheme fordetermining n( p) due to Gori-Giorgi and Ziesche (GZ) [27],who interpolate the accurate data for n( p) in the effective
potential expansion (EPX) method [9] for 1  rs  5 withthat in the Wigner-crystal limit (rs  10), together with thetwo sum rules, one for the total number and the other for the
total kinetic energy [9], expressed asV) ]S(q,ω)/N [units of (27.21e

YASUTAMI TAKADA n0 ( p + q) − n0 ( p)pσiω − ε p+q + ε p,(9)with n0 ( p) = θ (pF − | p|) the step function. The obtained(q,iω) can be used to determine both ω̃ex (q) and fex (q)by the analysis of its divergent property.
Finally, we need to make an analytic continuation of
(q,iω) into (q,ω + i0+ ) to calculate S(q,ω). There are
two ingredients, WI (q,ω + i0+ ) and Gs (q,ω + i0+ ), in thisanalytic continuation. The former can be obtained by the direct
calculation through Eq. (7) with ω + i0+ in place of iω. As forthe latter, we find it better to slightly modify Gs (q,iω) from theoriginal form [28], so that we can not only obtain much moreaccurate εc in a wide range of rs by the adiabatic connection
formula [29] but also eliminate unphysical divergences fromGs (q,ω + i0+ ).
In accordance with the above strategy, we have successfully
performed the calculation and obtained the results of both n( p)and S(q,ω) for rs up to over 20 to reveal that a giant peakother than the conventional plasmon one emerges in S(q,ω)
for rs  10, evolving into a predominant peak for rs > 20, as
shown in Fig. 1. Concomitantly a singular behavior is foundin n( p) for rs  22, which might suggest the occurrence ofan electronic phase transition, though at present we cannot
draw a definite conclusion on this point. Because the newpeak in S(q,ω) is evolved from the shoulder structure dueto the excitonic effect found at usual metallic densities [25],this peak may be called “an excitonic collective mode,” a newconcept in the dilute electron gas. In the dielectric-catastrophe
regime, we also find that εR (q,ω) is controlled by the valueof q ex (= 0) at which ω̃ex (q ex ) = 0, bringing about a singularpoint in εR (q,ω) at ω = 0 and q = q ex for rs > rsc , a featurenever seen in the conventional metals like the electron gas forrs < rsc .

245106-2

EMERGENCE OF AN EXCITONIC COLLECTIVE MODE IN . . .

PHYSICAL REVIEW B 94, 245106 (2016)

This paper is organized as follows: In Sec. II we derive thethird sum rule for n( p), with which we improve on the GZ
scheme to calculate n( p) in the electron gas. Our improved
scheme works very well to obtain n( p) up to rs ≈ 22. In
Sec. III we calculate S(q,ω) with use of the obtained n( p)and the local-field factor Gs (q,ω + i0+ ), the latter of which is
improved on the original form of Richardson and Ashcroft(RA). For rs larger than about 10, a novel peak structureappears in S(q,ω) and it dominates the plasmon peak for
rs > 20. The origin of this peak is investigated by a detailed
study of εR (q,ω). The results obtained in this paper are
summarized in Sec. IV, together with discussions on related
issues. In this paper, we employ units in which  = kB = 1.

A. Hamiltonian

Apart from a constant term, the first-principles Hamiltonian
H for a system of N electrons interacting with each otherthrough the Coulomb interaction in a periodic one-bodypotential is written in the plane-wave basis as [30]
(10)

with
H0 =



ε p c+pσ c pσ ,

(11)

pσ

V =

1 
V (q)c+p+qσ c+p −qσ c p σ c pσ ,2 q =0 pσ p σU=U p, p+G c+pσ c p+Gσ ,

(12)
(13)

pσ G =0

where c pσ is the annihilation operator of an electron with
momentum p and spin σ and G is a reciprocal-lattice vector.
If the electron density n (= N/ t ) is so low as to satisfy thecondition of pF a0  1, we may virtually regard |G|, whichis 2π/a0 or larger, as an infinite number, allowing us toneglect all terms in U . Then the system can be described by
the Hamiltonian reduced in the form of H = H0 + V , which
is nothing but the electron gas and with which we shall beconcerned in this paper.
B. Momentum distribution function

At zero temperature (T = 0), the momentum distributionfunction n( p) is defined as an expectation valuen( p) = c+pσ c pσ ,

10

g (0)

10
10
10

-2

g (0) :Overhauser (1995)

-4

O

-6

gY (0):Yasuhara (1972)
gSND(0):-8Spink, Needs,-10 & Drummond
10(2013)

020g(0):Gori-Giorgi& Perdew (2001)

GP

40 r 60s

80

100

FIG. 2. Comparison of various formulas for the on-top pair
distribution function g(0) plotted as a function of rs over a very
wide range.

II. MOMENTUM DISTRIBUTION FUNCTION

H = H0 + V + U,

0.5

where g(0) is the on-top value of the pair distribution function.
We can rewrite Eq. (15) as
8  αrs 2 pF8
lim n( p) =
g(0),
(16)
| p|→∞
9 π| p|8
dictating the asymptotic behavior of n( p). By using V (q) in
the electron-electron ladder approximation, Yasuhara derived
a formula for g(0) as [32,33] √24αrs /π1,(17)g(0) = gY (0) ≡√2 I1 (2 4αrs /π )where I1 (x) is the modified Bessel function of the firstkind. There are three other parametrization formulas for
g(0); by the phase-shift analysis of the two-electron problem,
Overhauser [34] derived one formula gO (0) and Gori-Giorgi
and Perdew [35] improved on it to give another one gGP (0),both of which are intended to provide accurate results of g(0)for rs  10. More recently, Spink, Needs, and Drummond(SND) [6] provided yet another formula gSND (0), which
reproduces their DMC data for rs  20. The results of thoseformulas are given as a function of rs in Fig. 2. We see
immediately that they are essentially the same for rs  10. For
rs  20, Yasuhara’s and SND’s give about the same results, but
the latter deteriorates very much soon after rs goes beyond 20.
[Actually, SND gives unphysical negative g(0) for rs > 24.26.]
Thus, apparently, the best choice for g(0) is Yasuhara’s formula
gY (0) in Eq. (17) which will be mainly employed over a wide
range of rs hereafter.
C. Sum rules(14)

evaluated with respect to the ground state 0 . This quantity is
independent of spin for a paramagnetic 0 .
The short-range electron correlation determines the asymptotic behavior of n( p) as well as the structure factor S(q) andtheir exact forms are given as [31]
 a 2a BBlim [|q|4 S(q)]lim [| p|8 n( p)] =t
4π | p|→∞8π |q|→∞= g(0),(15)

Let us derive the three sum rules, Eqs. (3)–(5), for n( p). The
static charge response function QRc (q,0) in the electron gas iscalculated by Moroni, Ceperley, and Senatore (MCS) [36] in
DMC for rs  10, giving the numerical data for (q,0). Withuse of Eq. (8), we obtain1
1=− V (q)G+ (q,0).
(q,0)
0 (q,0)(18)This may be regarded as the definition of the static local-field
factor G+ (q,0). The DMC data of (q,0), together with the245106-3

YASUTAMI TAKADAPHYSICAL REVIEW B 94, 245106 (2016)

well-known analytical result of 0 (q,0), determine G+ (q,0),
which is parametrized as a function of rs and q (= |q|) [36].
In particular, in the limit of q → ∞, G+ (q,0) is known as 2q
G+ (q,0) −−−→ C(rs )+ B(rs ),
(19)
q→∞
pF
where the constants B(rs ) and C(rs ) are given by√
1 + (a1 + a2 rs ) rs
B(rs ) = BMCS (rs ) ≡√ ,3 + (b1 + b2 rs ) rs∂εc (rs )π,C(rs ) = αrs −εc (rs ) − rs
4
∂rs

KE ≡

q→∞

2
[1 − g(0)].
3

(20)
(21)

(23)

(25)

(27)

A similar asymptotic expansion can be made for 0 (q,0) withIn replaced by In(0) = 1/(n + 1). Then, by comparing bothsides of Eq. (24) order by order and taking the large-q limit,
we obtain the following set of equations:

I6 =

1
,
3

(28)

∂εc (rs )1 α 2 rs2,+
−εc (rs ) − rs
5
3
∂rs

(29)

5
8
5αrs
2 2
+ I42 +
B(rs ) − + g(0) .
105 3
9π
3 3

mpF
I0 .
π2

(33)

1−x
αrs G(0)√2π 2 n0 − n− 4αrs /π
π 2 π 1 − ln 2 n0 − n− (1 − x)2, (35)+ b(rs )αrs 3 F (0) G(0)xx− = a− (rs )with F (0)/2 = 8.984373 andx+ = a+ (rs )

0

I4 =

(32)

where x− and x+ are, respectively, introduced as(22)with x = | p|/pF , z = |q|/2pF , and n(x) = n( p). This equation can be asymptotically expanded as
mpF I2
1 I4
1 I6
+
+
+ ··· ,
(26)
WI (q,0) = t 2
π
z2
3 z4
5 z6

I2 =

3I6 − 9I42 EF ,Following Gori-Giorgi and Ziesche (GZ) [27], we considern(x) in the parametrization form of⎧n0 − n−⎪⎨n0 −G(x− ) for x < 1,G(0)n(x) =(34)n+⎪⎩
G(x+ )
for x > 1,
G(0)

1
1
=
+ [Gs (q,0) − G+ (q,0)]V (q). (24)
WI (q,0)0 (q,0)where the nth-power integral In is defined as ∞In =dx n(x) x n .

(ε p − ε p )2  =D. Parametrization form of n( p)By combining Eq. (22) with Eq. (18), we obtain

Now WI (q,0) in Eq. (7) is reduced to



mpF ∞
x  x + z 
,
ln
dx
n(x)
WI (q,0) = t
2π 2 0z x − z

(31)

q→0

where the other local-field factor Gs (q,0) was investigated by
Niklasson [37] through the equation of motion method and is
known to behave asymptotically as
Gs (q,0) −−−→

KE ≡ ε p  = 3I4 EF ,

WI (q,0) −−→ WI (0,0) = twith a1 = 2.15, a2 = 0.435, b1 = 1.57, b2 = 0.409, and εc (rs )(in Ry) in the parametrization scheme due to Perdew andWang [4].
On the other hand, we obtain another expression for
(q,0)−1 with the use of Eq. (6) as1
1=− V (q)Gs (q,0),(q,0)WI (q,0)KE, the fluctuation of kinetic energy KE, and WI (0,0)
as

(30)

These equations correspond to Eqs. (3)–(5), respectively, andthus the sum rules are now proved. Note that with use of In and
EF = 1/(α 2 rs2 ) (in Ry), we can write the average kinetic energy+αrs G(0) x − 1√
2π 2 n+4αrs /π

3π (1 − ln 2) n+ π
(x − 1)4 F∞ (x).
g(0)G(0) 4αrs(36)Here G(x) is the Kulik function [38], defined as ∞R (u) u R(u) − R(y) G(x) =
du, (37)
√R(u) u + yu−y0
y=x/ R(u)with R(u) = 1 − u arctan(1/u). This function appears in thecalculation of n(x) in the random-phase approximation (RPA)
and its asymptotic behavior may be summarized as follows; inthe small-x limit,π √ 
+ 3 x ln x + · · · ,(38)G(x  1) = G(0) + π
4
with
 ∞
1
−R (u)
arctan = 3.353337 · · · , (39)duG(0) =R(u)u0while in the large-x limit, we obtainπ
1(1 − ln 2) 2 + · · · .
(40)6x
In accordance with this asymptotic behavior of G(x), n(x)in Eq. (34) behaves in the following way at various limits: In
the limit of x → 0,
F (0)  π 4 αrs 2 2
n(x) −−→ n0 −
x .
(41)
x→0
2
α
b(rs )G(x  1) =By a simple interpolation between high- and low-density limitsof the coefficient of the x 2 term in Eq. (41), GZ deduced b(rs )245106-4

EMERGENCE OF AN EXCITONIC COLLECTIVE MODE IN . . .

PHYSICAL REVIEW B 94, 245106 (2016)

to be the form of

1.0

1/2

.
b(rs ) = 1 + 0.0009376925rs13/4

n 0(rs ), n−(rs ), n (rs )We adopt this form for b(rs ). In the limit of x → 1,n(x) −−−−−→
n± ± A± (rs )|1 − x| ln |1 − x|,
±
x→1+0

0.8

+

(43)

(x − 1)2 + exp[a∞ (rs )],(46)(x − 1)2 + exp[−a∞ (rs )]where a∞ (rs ) is an adjustable parameter related to thecoefficient of the next-leading term. Note that F∞ (x) ≡ 1 fora∞ (rs ) ≡ 0. In actual calculations, we find that a∞ (rs ) plays
an important role in fulfilling the third sum rule, Eq. (30).
F∞ (x) =E. Input data for the parametrization scheme

There are six free parameters, namely, n0 , n± , a± , and a∞ , inthe definition of Eq. (34). The basic idea of GZ is to determine
the first three by the interpolation between the accurate data ofn( p) for 1  rs  5 obtained by EPX [9] and the result in thelimit of rs → ∞ in which GZ assume that n( p) is reduced tothat in the Wigner lattice, nW ( p), expressed as 2 3/2


pF
4π
p2
,
(47)
exp −
nW ( p) =3 π mω
mω
−3/2

with ω = 0.8833rs
in a.u. In improving on the GZ scheme,
we also adopt nW ( p) as basic information for rs  10, but
the EPX data are upgraded by the data for 1  rs  10
obtained by GW (or actually GW̃ WI ) [10]. Then the actual
parametrization formulas for n0 (rs ) and n± (rs ) are slightly
different from those in GZ. Our revised results for them aregiven, respectively, by
5/2

1 + t1 rs2 + t2 rs

13/4

1 + t3 rs2 + t4 rs

(48)

,

0.4
n − (rs)

0.2

n (rs )

020+

If the electron-hole symmetry strictly holds in the electronic
excitations, we obtain A− (rs ) = A+ (rs ), as was assumed by
GZ, but at low densities in which electron-hole excitations are
not restricted only near the Fermi surface and/or the dispersion
of a quasielectron is different from that of a quasihole, thesymmetry will be broken in general, requiring us to determine
A− (rs ) and A+ (rs ) separately. Finally in the limit of x  1,
8  αrs 2 1 g(0)n(x) −−−→
.
(45)
x→∞ 9
π
x 8 F∞ (x)2
If F∞ (x) approaches unity at x → ∞, Eq. (45) is reduced toEq. (16), as it should be. In fact, GZ took F∞ (x) ≡ 1, but we
can include the effect of the next-leading term of O(x −10 ) atx  1 by giving a form for F∞ (x) as

n 0(rs)0.6

with A± (rs ) “the Fermi-edge coefficient” or the coefficient of
the logarithmic singularity of the derivative of n(x) at x =1 ± 0+ , given by1  π √  αrs 1/2
+ 3
A± (rs ) =a± (rs )4 4π√≈ 0.256 rs a± (rs ).
(44)

n0 (rs ) =: Present parametrization: GZ parametrization: Data by GWΓ; : Data by EPX(42)

40 r 60s

80

100

FIG. 3. Parametrized n0 , n− , and n+ as a function of rs in bothpresent (solid curves) and Gori-Giorgi and Ziesche (dashed curves)
schemes. The input data obtained by EPX (triangles) and GW (solid
circles) are also shown.

with
t1 = 0.013 813 294 1,
t2 = 0.006 504 281 94,
t3 = 0.025 275 492 1,
t4 = 0.001 015 523 77,v1 =0.198 200 080, v2 = 0.001 839 431 25, v3 = 0.286 719 080,
v4 = 0.013 899 539 9,
v5 = 0.000 473 501 718,
q1 =0.088 519 000, q2 = 0.002 142 766 78, q3 = 0.387 350 521,
q4 = 0.079 971 845 3, and q5 = 0.000 551 585 578.
The obtained n0 (rs ) and n± (rs ) are plotted in Fig. 3. We seethat our results are about the same as those in GZ for rs < 10,
but they are considerably different, especially for the case ofn+ (rs ); ours change much more smoothly with the increaseof rs than those in GZ and eliminate the strange behavior ofn− (rs ) at rs ≈ 16 in GZ.
Incidentally, if we assume that n( p) converges to nW ( p)in the low-density limit, we should assess this assumption bychecking the three sum rules with this nW ( p) at rs → ∞. In
fact, the first two sum rules are easily found to be satisfied,
as already mentioned in GZ, but a problem exists on the2third one; this nW ( p) gives I4 =√ 1/(2εW ) and I6 = 5/(4εW )2with εW = pF /(mω) = 4.1692/ rs . By substituting those I4
and I6 into Eq. (30), we obtain limrs →∞ B(rs ) = 1.1868,
while limrs →∞ BMCS (rs ) = 1.064, indicating the difference inabout 10%. Because BMCS (rs ) was determined with the input
data only for rs  10, we consider it more appropriate toamend B(rs ) slightly from BMCS (rs ) into a form fulfilling theasymptotic value of 1.1868 as√
1 + [a1 + (a2 + a3 rs )rs ] rs
B(rs ) =(51)√ ,3 + [b1 + (b2 + b3 rs )rs ] rs
with a1 = 2.161, a2 = 0.4599, a3 = 0.006679, b1 = 1.594,
b2 = 0.4388, and b3 = 0.005627. In Fig. 4, this modified B(rs )is plotted in comparison with BMCS (rs ) to find that they are
essentially the same for rs  10. The function C(rs ) in Eq. (21)is also shown to see that C(rs ) correctly behaves asymptoticallyin the Wigner-lattice limit.

5/2

n− (rs ) =1 + v1 rs + v2 rs

13/4

1 + v3 rs + v4 rs2 + v5 rs

,

(49)

F. Calculated results(50)

The remaining three parameters, a± (or equivalently A± )and a∞ , are determined so as to accurately satisfy the three sum
rules at each rs . The results for n( p) calculated at rs = 8 with

5/2

n+ (rs ) =q1 rs + q2 rs
1/2

1 + q3 rs

3/2+ q4 rs

13/4+ q5 rs.

245106-5

YASUTAMI TAKADA

1.0

B: Wigner-lattice limit

B(rs ): Present proposal

0.8

BMCS (rs)

0.6

-1

10

rs =
1
2
3
4
5
6
8
10
12
14
16

n(p)

B(rs ), C(rs )100

PHYSICAL REVIEW B 94, 245106 (2016)

C(rs): Wigner-lattice limit

0.4

C(rs )

(a)

0.2
-2

With g Y (0)

10

-1

10

1

0

10 r 10
s

10

2

10

0.0

3

FIG. 4. Parametrized B(rs ) and C(rs ) as a function of rs . The
former is compared with the one proposed by Moroni, Ceperley, andSenatore, together with the value in the Wigner-lattice limit.

1.0

rs = 8

n(p)

0.8
0.6

0.2

: with g Y (0): with g GP (0): with g (0)
O
: GZ
: GWΓ

0.0

0.5

0.4

1.0

1.5

p (units of pF )2.0

FIG. 5. An example of n( p) calculated with use of g(0) in
Yasuhara (solid), GP (dotted-dashed), and Overhauser (dashed) forms
at rs = 8. Our results are compared with that in the GZ scheme aswell as that in GW .

18
20

0.6

1.0

1.5

2.0

p (units of pF )
rs=16

(b)
21
22

0.5

n(p)

use of three different forms, namely, Yasuhara, Overhauser,
and GP, of g(0) are shown in Fig. 5, together with those in the
GZ scheme and GW [10]. All these results are essentiallythe same, assuring that the present scheme works very wellto obtain accurate enough n( p) satisfying three sum rules
simultaneously.
In Fig. 6, the results of n( p) calculated with the choice ofg(0) = gY (0) are shown for a wide range of rs . In the metallicdensity regime (1  rs  5), the obtained n( p) is about thesame as that given in EPX [9] or GW [22]. Even in the
dielectric-catastrophe regime, n( p) behaves in a normal andwell-known manner, as long as rs is smaller than about 12. With
the further increase of rs , however, it begins to behave rather
differently from the normal one for | p| < pF and eventuallyfor rs > 20, it exhibits a novel feature in the sense that n( p)becomes virtually flat for | p| < pF , implying that the quasihole

0.7

0.5

0.4
0.3
0.2
0.1

With g Y (0)

0.0

0.5

1.0

1.5

p (units of pF )2.0

FIG. 6. Obtained n( p) for (a) 1  rs  16 and (b) 16  rs  22.

excitation energy at the center of the Fermi sphere is about thesame as that near the Fermi surface. For | p| > pF , on theother hand, this anomalous flat behavior is never seen in n( p),indicating strong electron-hole asymmetric excitation spectra.
Incidentally, a sign of this asymmetry has already been seen in
the one-electron spectral function obtained in GW at rs = 8(see Fig. 3(b) in Ref. [10]) and here we find that this tendencyof asymmetry becomes so enhanced for rs > 20 to provide thevery characteristic flat behavior of n( p) in Fig. 6(b).
In Fig. 7, we plot the results for A± and a∞ determined so
as to satisfy the three sum rules, Eqs. (28)–(30), at each rs .
Three different forms for g(0) are employed, but the obtained
A± is virtually independent of its choice; the difference in g(0)is mostly compensated by the difference in a∞ . The difference
between A+ and A− indicates the degree of electron-holeasymmetry near the Fermi surface; for rs < 12, they are
essentially the same, but for rs  12, A− rapidly increases,
while A+ does not change much. This asymmetry is probably
the physical reason why the GZ scheme, in which A− = A+
was assumed, could not give a convergent result of n( p) for
rs > 12.
In the parametrization scheme in Eq. (34), the behaviorof n( p) for | p| < pF is controlled only by the parameter A−
and thus its anomalous flat feature is directly connected with

245106-6

EMERGENCE OF AN EXCITONIC COLLECTIVE MODE IN . . .

PHYSICAL REVIEW B 94, 245106 (2016)

8

1.4

0.5
0.4

1.2

A-−1

<KE>/EF , ΔKE/EF

0.3
0.2

A+, A- , a∞

4

1.0

A-

0.1
018

19

2

20 r 21
s

22

23

A+
: with gY(0): with gGP(0): with gO(0): A + =A − in GZ

-4 0

5

10 r
s

15

0.2
0

20

FIG. 7. Parameters A± and a∞ determined so as to satisfy thethree sum rules for n( p) at each rs . Three different forms for g(0)are employed. For comparison, A± in GZ are also shown. (In GZ,
a∞ = 0.) The inset displays A−1− for the region of rs in which A− isvery large.

the rapid increase of A− for rs > 12. As shown in the inset
of Fig. 7, A− actually diverges at rs = r̄sc = 22.4 ± 0.4, theerror of which is deduced from the difference in r̄sc causedby the different choice of g(0). For rs beyond r̄sc , our scheme
does not work and fails to determine the three parameters,
A± and a∞ , simultaneously satisfying three sum rules. We
might be able to interpret this divergence of the Fermi-edge
coefficient or the nonanalytic behavior at the Fermi surface
as an indication of the occurrence of some non-Fermi-liquid
phase, but at present, we do not know exactly whether this is thecase or not, mainly because we cannot rule out the possibilityof inadequacy of the postulated form in Eq. (34) for the caseof A−  1 (or rs > 20). In any case, for rs  12, it is certainthat some unusual situation due to electron correlation begins
to appear, as we can see in Fig. 8 in which the fluctuation of
the kinetic energy KE becomes larger than the average one
KE.

with the retarded dielectric function εR (q,ω), given by
ε (q,ω) = 1 + V (q) (q,ω)
RR= 1 + V (q)0.6
0.4 I2

RWI (q,ω).
1 − GRs (q,ω)V (q)RWI (q,ω)

5

(53)

Here Eq. (6) is employed to write R (q,ω) in terms of
RWI (q,ω) and GRs (q,ω), the latter of which is given by the
analytic continuation of Gs (q,iω) to the real axis in the

0.2

I4

0.0
0

5

10 r
s

10 rs 15

15

20

20

complex upper ω plane (iω → ω + i0+ ). Thus our task to
calculate S(q,ω) begins with obtaining both RWI (q,ω) andGRs (q,ω).
R
B. Polarization function WI(q,ω)

The function WI (q,iω) in Eq. (7) is rewritten as
ε p+q − ε p
WI (q,iω) = 2
n( p) 2ω+
(ε p+q − ε p )2
pσ
mpF
PWI (z,iu),
2π 2with the function PWI (z,iu), introduced as ∞1
u2 + (x + z)2
PWI (z,iu) =,xdx n(x) ln 2
2z 0
u + (x − z)2
= t

(54)

(55)

where x = | p|/pF , z = |q|/2pF , and u = (ω/4EF )/z. By
analytic continuation, RWI (q,ω) is obtained as
RWI (q,ω) = t

At T = 0, with use of Eq. (2) for QRc (q,ω), S(q,ω) can becast into the form of
1 1
1
,
(52)
S(q,ω) = −Im R
π V (q)ε (q,ω)I6

FIG. 8. The average kinetic energy and the fluctuation of the
kinetic energy as a function of rs . No difference is seen by the different
choice of g(0) because of the fulfillment of the three sum rules. Inset:
Corresponding values of In for n = 0,2,4, and 6.

III. DYNAMICAL STRUCTURE FACTOR
A. Basic equation to determine S(q,ω)

I0

0.8

0.4

a∞

1.2
1.0

0.6 ΔKE

0
-2

<KE>

0.8

In

6

: with gY (0): with gGP(0): with g O(0)mpF
PWI (z,u + i0+ ),
2π 2

(56)

with PWI (z,u + i0+ ), given by ∞1
PWI (z,u + i0+ ) =xdx n(x){ln[u2 − (x + z)2 + i0+ ]2z 0− ln[u2 − (x − z)2 + i0+ ]},(57)

where ln(ω) is so defined as ln(ω) = ln(|ω|) + iarg(ω) witharg(ω) the argument of ω in the upper ω plane being between
0 and π . In performing the numerical evaluation of this integral
over the variable x, care must be exerted at both x = |u − z|
and x = u + z in order to accurately obtain the principal valuesof the integral.
With n( p)[= n(x)] obtained in Sec. II F, we can calculate
RWI (q,ω) through Eqs. (56) and (57). The overall structureof the calculated RWI (q,ω) does not depend much on rs , in

245106-7

YASUTAMI TAKADAPHYSICAL REVIEW B 94, 245106 (2016)
C. Dynamical local-field factor G sR (q,ω)

R
Re ΠWI(q,ω) /(Ωt mpF /2π 2 )

2.0
1.5

(a)

1.0

rs =22

0.5
0.0

-0.5
-1.0
-1.5
0

q=pF
q=0.5pF

2

4

2.5
R
Im ΠWI(q,ω) /(Ωt mpF /2π 2)There is a long history of research on the dynamical localfield factor and its parametrization [29,39–46]. Those studies
are mainly motivated by the construction of the exchangecorrelation kernel fxc (r,r ; ω) in the time-dependent density
functional theory (TDDFT) and thus they are concerned witheither G+ (q,iω) or GR+ (q,ω), but not with Gs (q,iω). So far aparametrization form for Gs (q,iω) is given only by RA [28],proposing: Π R (q,ω)
WIR: Π 0 (q,ω)

q=2pF q=2.5pF

q=1.5pF

2.0
1.5

q=pF

: Π R (q,ω)
WIR: Π 0 (q,ω)

rs =22as (i) = λ(∞)
+
s

q=2pF

0.5

q=2.5pF

=
λ(∞)
s

0.00

2

4

(∞)
λ(0)
s − λs
,1 + β1 γs  + (β2 γs )2(∞)
where λ(0)are, respectively, defined ass and λs


1
π
κF
(0)
λs =
−
,
αrs I0κ

q=1.5pF

1.0

as (i)z2 + 2[1 − g(0)]bs (i)z8 /3,
1 + cs (i)z2 + bs (i)z8(58)

with  ≡ zu = ω/4EF . This form of Gs (q,iω) satisfies
the rigorous limit due to Niklasson [37] at |q| → ∞ asGs (q,iω) → 2[1 − g(0)]/3, irrespective of ω.
The function as (i) is so determined as to satisfy theconstraints in the small-|q| regime, namely, the compressibility
sum rule at  = 0 and the third-moment sum rule [47–49] at
 → ∞, leading to the following form:

6 8 10 12 14ω (units of E F )(b)
q=0.5pF

Gs (q,iω) =6 8 10 12 14ω (units of E F )

(59)

(60)

∂εc (rs )3 2+ 2εc (rs ) ,− π αrs rs
5 5
∂rs

(61)

where the ratio of the compressibilities with and without theinteractions, κ and κF , can be calculated as

FIG. 9. The retarded polarization function RWI (q,ω) in comparison with R0 (q,ω) at rs = 22 for (a) real and (b) imaginary parts.

spite of the appearance of an anomalous feature in n( p) for
rs > 20; as shown in Fig. 9, even at rs = 22, the results of
RWI (q,ω) given as a function of ω at various values of q (thesolid curves) retain the main features of those of R0 (q,ω)
the Lindhard function in the retarded form (the dotted-dashedcurves). We can understand this insensitivity of RWI (q,ω) to
the anomalous feature in n( p) by paying attention to the factthat RWI (q,ω) is mostly controlled by the sum rules for n( p);the details of n( p) are smeared out by the integral over thevariable x in Eq. (57).
Note, however, that there is a very important difference
between R0 (q,ω) and RWI (q,ω); ImR0 (q,ω) does not vanish
only in the electron-hole single-pair excitation region in (q,ω)
space, while basically ImRWI (q,ω) is not zero but positive
definite for ω > 0 at any q due to the inclusion of the multiple
electron-hole pair excitations by the use of n( p) instead ofn0 ( p). For this reason, as long as Eq. (53) is employed, we
do not need to introduce an artificial finite broadening width
γ in the numerical calculation of the dielectric function, in
sharp contrast with the usual treatment of ω → ω + iγ in thecalculation of the plasmon-peak structure in S(q,ω) with using
R0 (q,ω) or its simple extension in which no account is takenof the multiple electron-hole pair excitations.

κFαrs
α 2 rs3∂ 2 εc (rs )∂εc (rs )
.
=1−
+
rs
−2
2
κ
π
6
∂rs
∂rs

(62)

As for the two functions, cs (i) and bs (i), in Eq. (58), by theconsideration of the condition to obtain the maximum valueof 4[1 − g(0)]/3 for Gs (q,iω) at z2 ≈  → ∞, RA gave thefollowing forms:
(∞)

λs
4
− α1RA + 34 1−g(0)3 λ(∞)3scs (i) =−,4 1 − g(0)1 + γs 

bs (i) =(63)

as (i)3as (i) − 2[1 − g(0)][3cs (i) + 4/(1 + )]/3
1
×
,
(64)
1 + γ1  + γ2 2 + γ3 3 + 4

with the parameter γs determined as
γs =

3 1
9 λ(∞)s+1−.
16 1 − g(0)
4 αRA

(65)

In Eqs. (59), (63), and (64), six parameters, αRA , β1 , β2 ,γ1 , γ2 , and γ3 , are at our disposal, but RA made the following choice: αRA = 0.9, β1 = 0, β2 = 1, γ1 = γ3 = 4, and
γ2 = 6.
A similar parametrization scheme was proposed by RA for
Gn (q,iω) [≡ G+ (q,iω) − Gs (q,iω)]. As a matter of fact, with

245106-8

PHYSICAL REVIEW B 94, 245106 (2016)

2.0 Ghost Exciton Mode

Original Richardson-Ashcroft

20

1.8

10

~ (q) (units of E )ωexFεc − εcexact (units of meV)

EMERGENCE OF AN EXCITONIC COLLECTIVE MODE IN . . .

0Modified Richardson-Ashcroft:
Simple Version
Modified Richardson-Ashcroft:
Improved Version

-10
-20
0

5

10 rs 1520

25

rs =201.6

rs =16

1.4
1.2
1.0

rs =12
rs =8

0.8
0.6
0.4

FIG. 10. Difference between the “exact” value of εc due to
Perdew and Wang and the one calculated through the adiabatic
connection formula with using the local-field factor in the originalRA (the dotted-dashed curve), a simple version of the modified RA(the dashed curve), and an improved version of the modified RA (thesolid curve) schemes.

rs =5.3

rs=5.5

0
0.20

0.5

1.0
1.5
2.0
q (units of pF )(b) c /v = 2_ |κ /κ|
0.15

cex/vFuse of Eq. (24), Gn (q,iω) is given exactly as(q,iω)Gn (q,iω) = Gexact
n

F

π

F

rs=4.5rs =5.24

0.10

(66)

0.05

We have calculated Gexact(q,iω) through Eq. (66) withn
WI (q,iω) given by Eqs. (54) and (55) and compared with
GRA
n (q,iω), the parametrized form for Gn (q,iω) due to RA,
to find that GRAn (q,iω) is in fact very reliable; in particular,
no appreciable difference can be seen for the quantities given
by integrals over q and ω, such as εc (rs ) obtained through the
adiabatic connection formula, Eq. (27) in Ref. [29], in which
we have calculated with use of either Gs (q,iω) + GRA
n (q,iω)(q,iω).
With
thisassessmentofthe RAor Gs (q,iω) + Gexactnparametrization scheme, we can expect that the parametrizedGs (q,iω) will also be reliable. Incidentally, this RA scheme
comes to be known to possess a very good feature [29]; it
provides a quite accurate εc (rs ) in reference to its “exact” value
εcexact (rs ) given by Perdew and Wang [4] (see the dotted-dashedcurve in Fig. 10).
By analytic continuation (iω → ω + i0+ ), we obtainR
Gs (q,ω) from Gs (q,iω) in Eq. (58), but with βi s chosen
in RA, GRs (q,ω) at small q suffers from the existenceof unphysical divergence at ω = ±4EF /γs . A simple way
to eliminate this unphysical divergence while keeping thegood feature of providing accurate εc (rs ) is found to choose
αRA = 0.958/(1 + 0.006rs ), β1 = 0.8, and β2 = 0.4 with γi s
unchanged. This simple version of the modified RA scheme
gives better εc (rs ) than the original RA scheme (see Fig. 10). A
closer analysis of this scheme, however, reveals that unphysicaldivergence appears for 2pF  |q|  3pF with γi s chosen as
those in RA.
Under these circumstances, we have made an extensive search of suitable values for all free parameters in
Eqs. (59), (63), and (64) with considering the following four
criteria: (i) GRs (q,ω) contains no unphysical divergence forreal ω, (ii) GRs (q,ω) should be analytic in the upper complex
ω plane, in accord with the causality condition [39,40], (iii)
Gs (q,iω) accurately provides εc in a wide range of rs , and (iv)GRs (q,ω) combined with RWI (q,ω) gives the results of S(q,ω)

rs =6

0.2

ex

1
1
1
≡
−
.
V (q) 0 (q,iω) WI (q,iω)

rs =22(a)0

rs=5.2
rs=5.1
rs=5.0rs =4.8

2.5

rs =4.0
rs=6.5rs=6.2

rs=6.0rs=5.85
rs=5.7
rs =5.5 rs =5.4~
ω (q)/qcex=lim
rs =5.3
q_>0 ex
rs =5.26

0.1 |κ /κ| 0.2
F0.3

FIG. 11. (a) Dispersion relation of the ghost exciton mode atvarious values of rs . (b) Sound velocity of the ghost exciton mode inunits of the Fermi velocity vF plotted as a function of |κF /κ|.

in good agreement with those in GW [25] for 1  rs  5.
As a result of this search, we have chosen the following set ofparameters, which we call an improved version of the modified
RA scheme: β1 = β2 = 0.8, γ1 = γ2 = γ3 = 1, and
α0


 ,
αRA (rs ) =(67)
1 + α1 + α2 / 1 + α3 rs + α4 rs2 rs
with α0 = 1.2110, α1 = 0.01064, α2 = 0.10296, α3 =−0.2910, and α4 = 0.2540.
D. Ghost exciton mode ω̃ex (q)With WI (q,iω) and Gs (q,iω) thus determined, we cancalculate (q,iω) through Eq. (6) to find that it actually
diverges at iω = ±i ω̃ex (q) for rs > rsc in accord with thediscussion in Ref. [21]. Numerically, ω̃ex (q) can be determinedby the search of zero of (q,iω)−1 for ω > 0.
In Fig. 11(a), the obtained results of ω̃ex (q) are shown
as a function of q (= |q|) at several values of rs in the
dielectric-catastrophe regime of rs > rsc . For small q, ω̃ex (q)increases in proportion to q, just like a sound mode, and “thesound velocity” cex changes in proportion to |κF /κ| withthe coefficient 2vF /π (vF : the Fermi velocity) for rs nearrsc , as shown in Fig. 11(b). Note that this behavior of cex isquite similar to that for “the exciton mode” identified by thepeak position of Im εR (q,ω) for rs < rsc [20]. Thus it is more
appropriate to call the mode of ω̃ex (q) on the imaginary ω axis

245106-9

YASUTAMI TAKADAPHYSICAL REVIEW B 94, 245106 (2016)

rs =2210
rs =20rs=16

fex(q)

8
6

1.0
c /v
0.8 ex F
0.6
0.4
Fraction ofstrength in0.2
the f-sum rule
0.0
5 10 rs 15 20

cex/vF, Fraction12

rs =12
rs=8

4

rs=6

2

rs =5.5
rs =5.3

0

0.5

1.0
1.5
2.0
q (units of pF )2.5

FIG. 12. Oscillator strength of the ghost exciton mode. Inset:
Fraction of this oscillator strength to the full strengths in the limit ofq → 0, together with the sound velocity in units of vF , plotted as afunction of rs .

“the ghost exciton mode” than “the ghost plasmon” which was
suggested in Ref. [21].
The oscillator strength fex (q) in the definition of a (q,iω)in Eq. (1) is determined by the evaluation of [ω −ω̃ex (q)](q,iω) in the limit of ω → ω̃ex (q) and the obtained
results are shown in Fig. 12. Incidentally, the asymptotic
behavior of (q,iω) in Eq. (6) is found to be(q,iω) −−−→
ω→∞

4π ne2 1,
V (q) m ω2

(68)from which we can prove the f -sum rule for S(q,ω) as ∞N q2dω ωS(q,ω) =
.
(69)
2m
−∞
On the other hand, a (q,iω) behaves asa (q,iω) −−−→
ω→∞1 ne2 fex (q).
V (q) m ω2

(70)

Comparison of Eq. (70) with Eq. (68) indicates that the ghost
exciton mode contributes to the full oscillator strengths withthe fraction of fex (q)/4π , which becomes more than 80% at
q → 0 for rs > 14. Thus, in such a dilute electron gas, wecome to notice that not the plasmon but the exciton mode willplay a main role in the whole excitation spectra.
E. Calculated results of S(q,ω)

RWI (q,ω)

With
and GRs (q,ω) already given, we cancalculate S(q,ω) through Eqs. (52) and (53). As mentioned
before, in the usual metallic region (1  rs  5), the form ofGRs (q,ω) is finely tuned to reproduce the results of S(q,ω) thatare already obtained in Ref. [25]. As illustrated in Fig. 1(a)for rs = 4, the structure of S(q,ω) is featured by the single
plasmon peak and the exciton contribution appears not as apeak but only as a shoulder structure in the low-ω region. (To

FIG. 13. The dynamical structure factor at (a) rs = 8 and (b)rs = 16. Inset: Bird’s-eye view.

be more definite, see the structure “a” specified in Fig. 1 inRef. [25].)
Even in the dielectric-catastrophe regime (rs > rsc ), themain feature of S(q,ω) hardly changes at rs near rsc , but as
rs increases further, the shoulder structure due to the excitonic
effect gradually evolves into a broad peak, as illustrated inFig. 13(a) at rs = 8. The peak becomes sharper with thefurther increase of rs and appears as a clear peak for rs ≈ 10.
Eventually for rs  14, combined with the plasmon peak,

245106-10

2.0
1.8
1.6

10

Excitonic Collective Mode(a) Re εR(q,ω) @ rs=4

rs =5.3
rs =5.5
rs =6

1.4
1.2
1.0

8

rs =8

0.8
0.6

rs =12
rs=16
rs =20
rs=22

0.4
0.2
0

PHYSICAL REVIEW B 94, 245106 (2016)ω (units of EF)

ωex (q) (units of EF )EMERGENCE OF AN EXCITONIC COLLECTIVE MODE IN . . .

0.5

1.0
1.5
2.0
q (units of p F)

64

plasmon2
2.5

0FIG. 14. Dispersion relation of the excitonic collective mode as
determined from the peak position in S(q,ω) in (q,ω) space. The
shaded area represents the electron-hole single-pair excitation region
and the dotted curve indicates the boundary of ω = (2pF q − q 2 )/2m.

F. Retarded dielectric function ε R (q,ω)

Let us make a detailed analysis of the retarded dielectric
function εR (q,ω) in order to investigate the relation between
the excitonic-peak position (or the excitonic collective mode)ωex (q) and the ghost exciton mode ω̃ex (q).
In the conventional metals without the dielectric catastrophe
as represented by the electron gas at rs = 4, the overallbehavior of εR (q,ω) is well known and rather simple; as shownin Fig. 15(a) for its real part in (q,ω) space, there is only a single
singular point at q = ω = 0 in this function, characterized by2
εR (q,ω) ≈ 1 − ωpl/ω2 in the limit of ω → 0 with keeping2/q 2
vF q/ω  1 (the ω limit) and εR (q,ω) ≈ 1 + (κ/κF )qTF
in the limit of q → 0 with vF q/ω  1 (the q limit), where

1.0

1.5 2.0 2.5
q (units of pF)

3.0

10

(b) Re εR(q,ω) @ rs =88
ω (units of EF)

S(q,ω) is characterized by a twin-peak structure, as shownat rs = 16 in Fig. 13(b). At rs = 20 or larger, “the excitonicpeak” grows even bigger than the plasmon one [see Fig. 1(b)].
In Fig. 14, we have plotted ωex (q) the excitonic peak
position (including the shoulder position, if a peak structure is
not well identified) at each q in (q,ω) space. This plot showsthat ωex (q) is always inside the single-pair excitation region, as
it should be for the excitonic effect working on an electron-holesingle-pair excitation.
An even more interesting fact is that if ωex (q) is in theregion of ω  (2pF q − q 2 )/2m, the corresponding peak inS(q,ω) is very sharp. In particular, ωex (q) becomes lowest if
it lies just on this boundary and the highest peak in S(q,ω)
appears in its very vicinity. Furthermore, with the increase ofrs , this lowest-energy peak position (qmin ,ωex (qmin )) becomes
lower approximately in proportion to rs along the boundary
curve, so that with the further increase of rs up to about
30, we may imagine a critical situation of qmin → 2pF andωex (qmin ) → 0 by extrapolation. If this critical situation wereactually realized, we might expect the occurrence of a very
exotic phase transition brought about by the spontaneous
excitation of a macroscopic number of excitons or “the exciton
condensation.”

0.5

6excitonic collective
plasmon4
2

ghost exciton

00.5

1.0 1.5 2.0 2.5
q (units of pF)

3.0

FIG. 15. Two-dimensional contour map of the real part ofεR (q,ω) at (a) rs = 4 and (b) rs = 8.

ωpl (= 4π e2 n/m) is the plasmon energy at q → 0 andqTF (= 4e2 mpF /π ) is the Thomas-Fermi screening constant.
Associated with this singular point at the origin of (q,ω)
space, the plasmon emerges as a collective mode satisfying
εR (q,ω) ≈ 0 at ω ≈ ωpl for small q.
In the dielectric-catastrophe regime (rs > rsc ), there appears
the ghost exciton mode, but as long as ω̃ex (q) is positive, themode is not a pole on the real axis in the complex ω plane,
indicating its irrelevance from a physical point of view. At
both q = 0 and qex (= 0) at which ω̃ex (q) becomes zero [seeFig. 11(a)], however, the pole is situated on the real axis andthus it must have direct physical relevance.
In fact, we find that εR (q,ω) contains a couple of singularities, one at the origin of (q,ω) space and the other at
q = qex and ω = 0, as illustrated in Fig. 15(b) for the case of
rs = 8. The singularity occurring at the origin behaves in muchthe same way as that in the conventional metals, excluding a
possibility of the observation of any exotic effects associatedwith this singularity. The only difference in the case of negative

245106-11

YASUTAMI TAKADAPHYSICAL REVIEW B 94, 245106 (2016)10(a) Re εR(q,ω) @ rs=22
ω (units of EF)

8
6 excitonic collective
plasmon4

ghost exciton

2
0

0.5

1.0 1.5 2.0 2.5
q (units of pF)

3.0

10

(b) Im εR(q,ω) @ rs=22
FIG. 16. (a) Calculated qex in units of pF , Aex , and Bex as afunction of rs . (b) Comparison of εR (q,ω) between the full calculationin Eq. (53) and the approximate one in Eq. (71) at rs = 22 forω = 0.24EF at which energy S(q,ω) becomes maximum due to theexcitonic collective mode.

where the coefficients, Aex and Bex , are given by Aex =pF ∂f (qex ,0)/∂qex and Bex = − limω→0 EF Imf (qex ,ω)/ω withf (q,ω), defined in reference to Eq. (53) by
 2
pF
1
π 1
−
f (q,ω) =GRs (q,ω).
(72)
2 αrs PWI (q,ω)
q
Note that f (qex ,0) = 0 by the definition of qex .
In Fig. 16(a), the obtained qex , Aex , and Bex are plotted
as a function of rs . At rs = rsc , Aex vanishes, removing
this singularity from εR (q,ω), but Aex is positive definite
for rs > rsc and thus this singularity always exists in the
dielectric-catastrophe regime. With the increase of rs , qex
increases monotonically and at rs ≈ 9.51 it becomes equal
to 2pF , which is the boundary at the electron-hole single-pair
excitation. Due to this transition across the boundary, both Aex
and Bex exhibit anomalous behavior at rs ≈ 9.51. Because Bex
is much decreased at rs ≈ 9.51 and becomes smaller further for
rs > 9.51, the effect of this singularity on εR (q,ω) is graduallyenhanced with the increase of rs .

ω (units of EF)κ from that of the conventional metals is seen in the sign ofRe εR (q,ω) in the q limit and concomitantly the reduction of
the oscillator strength for the plasmon pole, but this reduction
never becomes perfect (see Fig. 12), allowing the existence ofthe plasmon at any value of rs .
The singularity at q = qex and ω = 0, on the other hand, has
not been discussed in the literature and thus a detailed analysis
is needed here. In the vicinity of this singularity, εR (q,ω) isexpressed as
 2
1
pF
, (71)εR (q,ω) ≈ 1 +q
Aex (q − qex )/pF − iBex ω/EF

8
6

plasmonexcitonic
collective
ghost exciton

4
2

0

0.5

1.0 1.5 2.0 2.5
q (units of pF)

3.0

FIG. 17. Two-dimensional contour map of εR (q,ω) at rs = 22 for(a) real and (b) imaginary parts.

In Fig. 16(b), we have given an example of the comparison
between the approximate form for εR (q,ω) in Eq. (71) and thefull result of εR (q,ω) in Eq. (53) calculated at rs = 22, fromwhich we can confirm the accuracy of this approximate form in
the neighborhood of the singular point at (q,ω) = (qex ,0). This
figure also shows that the excitonic-peak position providing
the maximum height in S(q,ω), (qmax ,ωex (qmax )), is actuallycoincident with the point in (q,ω) space at which εR (q,ω) ≈ 0,
validating to call ωex (q) “the excitonic collective mode.” In
this way, the excitonic collective mode is brought about by the
singular behavior of εR (q,ω) associated with the singularity at(qex ,0) which is originally induced by the ghost exciton mode.
Finally in Fig. 17, εR (q,ω) given in Eq. (53) is plotted at
rs = 22 in the whole (q,ω) space in order to fully represent
the typical behavior of εR (q,ω) in the dielectric-catastrophe
regime. As a general feature, εR (q,ω) contains a coupleof singular points. Associated with each singularity, there
appears a collective mode at which εR (q,ω) ≈ 0, leading to the245106-12

EMERGENCE OF AN EXCITONIC COLLECTIVE MODE IN . . .

PHYSICAL REVIEW B 94, 245106 (2016)

double-peak structure in S(q,ω), one for the plasmon and theother for the excitonic collective mode, as shown in Fig. 1(b).

et al. [51] found that HF described their experimental results
well in the metallic state at negative κ.
(iii) As mentioned in Secs. II F and III E, we find indicationsof some exotic electronic phase transition of the 3D electron
gas for rs > 20. We need to make a further study of thepossibilities of such a phase transition, but at the same time,
we need to take account of the possibility of a partiallyspin-polarized state for 20 < rs < 40 [5], suggesting that weinvestigate the fate of the excitonic collective mode in the
spin-polarized state as well as the coupling between charge
and spin channels, a difficult but very intriguing problem left
in the future.
(iv) Finally, we remark on the similarity and difference
between the plasma and the excitonic collective modes. Due
to εR (q,ω) ≈ 0 at each mode in the relation of D(q,ω) =εR (q,ω)E(q,ω), there exists a wavelike charge-density oscillation as an eigenmode, accompanied by a finite electricfield E(q,ω) even in the absence of an external force [i.e.,D(q,ω) = E(q,ω) + 4π P(q,ω) ≈ 0]. The associated electric
polarization P(q,ω) is made of collective electron-hole excitations and works as a restoring force. Those charge-density
oscillations can be quantized as bosons, just as the plasma
oscillations quantized into the plasmons. In spite of thosesimilarities, there is an important difference in the correlation
among excited electrons and holes: for the plasmons, because
its energy is larger than EF , all excited electrons and holes areuncorrelated, as can be well treated in RPA. For the excitonic
collective modes, on the other hand, the excitations occur asa collection of excitons, indicating no polarization field in thelong-wavelength limit and thus the absence of the mode at|q| → 0; the polarization appears, if |q| is of the order of 2pF
or the inverse of the exciton binding radius. Incidentally, if
we consider the exotic phase mentioned in (iii) in terms ofa spontaneous excitation of this charge-density oscillation,
it might be regarded as a CDW state in the sense of notnesting-driven but Overhauser [52,53].

IV. CONCLUSION AND DISCUSSION

In this paper, we have calculated the momentum distributionfunction n( p) and the dynamical structure factor S(q,ω) in thethree-dimensional (3D) electron gas at low densities, i.e., rs >5.25, in which the electronic compressibility κ is negative.
The obtained n( p) is considered to be sufficiently accurate,
because it rigorously satisfies not only the exact asymptotic
behavior but also the three sum rules, Eqs. (3)–(5), the lastof which is successfully derived in this work. The obtainedS(q,ω) is also considered to be sufficiently accurate, because
the two ingredients to construct S(q,ω), namely, RWI (q,ω)
and GRs (q,ω), are both reliably determined in accordance withall the known constraints that must be fulfilled for correct
evaluation of those quantities. For rs > 10, the calculatedS(q,ω) exhibits a peak structure due to the excitonic collectivemode in addition to the well-known plasmon peak. With
the increase of rs , this new peak structure grows steadily
and eventually for rs > 20 it dominates the plasmon peak.
Associated with this excitonic collective mode, a singularity
is found for the first time in εR (q,ω) at q ≈ 2pF and ω = 0,
which is the direct physical consequence of the appearance ofthe ghost exciton mode for rs > 5.25.
Four comments are in order:(i) Although we have treated the 3D electron gas in thispaper, exactly the same physics will be found in the 2D electron
gas in which the dielectric catastrophe regime appears for
rs > 2. Actually, negative κ has already been observed in the2D system [50,51]. Thus, in order to confirm the emergence ofthe excitonic collective mode by the measurement of S(q,ω)
in some suitably designed experiments, the 2D system may be
more recommended, although detailed calculation for the 2D
system is left for the future.
(ii) Physically, if the excitonic collective mode is a dominantpolarization process in the charge response to an external pointcharge perturbation, we may imagine that the screening effectwill be much more reduced than that in the usual metals due tothe fact that the excitation of tightly bound electron-hole pairs
does not contribute to screening, indicating that the HartreeFock (HF) approximation may work well in the evaluation ofsome physical quantities. This might be the reason why Ilani[1] See, for example, G. F. Giuliani and G. Vignale, Quantum
Theory of the Electron Liquids (Cambridge University Press,
Cambridge, 2005).
[2] D. M. Ceperley and B. J. Alder, Phys. Rev. Lett. 45, 566 (1980).
[3] S. H. Vosko, L. Wilk, and M. Nusair, Can. J. Phys. C 58, 1200
(1980).
[4] J. P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992).
[5] G. Ortiz, M. Harris, and P. Ballone, Phys. Rev. Lett. 82, 5317
(1999).
[6] G. G. Spink, R. J. Needs, and N. D. Drummond, Phys. Rev. B
88, 085121 (2013).
[7] G. Ortiz and P. Ballone, Phys. Rev. B 50, 1391 (1994).

ACKNOWLEDGMENTS

This work is supported by a Grant-in-Aid on Innovative
Area “Materials Design through Computics: Complex Correlation and Non-Equilibrium Dynamics” (Grant No. 22104011)
from Ministry of Education, Culture, Sports, Science, andTechnology, Japan.

[8] M. Holzmann, B. Bernu, C. Pierleoni, J. McMinis, D. M.
Ceperley, V. Olevano, and L. Delle Site, Phys. Rev. Lett. 107,110402 (2011).
[9] Y. Takada and H. Yasuhara, Phys. Rev. B 44, 7879 (1991).
[10] H. Maebashi and Y. Takada, Phys. Rev. B 84, 245134 (2011).
[11] See, for example, G. D. Mahan, Many-Particle Physics, 3rd ed.
(Kluwer Academic/Plenum Press, New York, 2000), p. 355.
[12] Y. Takada, Phys. Rev. B 43, 5979 (1991).
[13] D. Pines and P. Noziéres, The Theory of Quantum Liquids(Benjamin, New York, 1966), Vol. 1.
[14] K. Matsuda, K. Tamura, and M. Inui, Phys. Rev. Lett. 98, 096401(2007).

245106-13

YASUTAMI TAKADAPHYSICAL REVIEW B 94, 245106 (2016)[15] H. Maebashi and Y. Takada, J. Phys. Soc. Jpn. 78, 053706(2009).
[16] H. Maebashi and Y. Takada, J. Phys.: Condens. Matter 21,
064205 (2009).
[17] O. V. Dolgov, D. A. Kirzhnits, and E. G. Maksimov, Rev. Mod.
Phys. 53, 81 (1981).
[18] C. A. Kukkonen and A. W. Overhauser, Phys. Rev. B 20, 550(1979).
[19] Lack of knowledge on this structure of the effective electronelectron interaction erroneously leads to proposing various typesof CDW instability; for example: A. M. J. Schakel, Phys.
Rev. B 64, 245101 (2001).
[20] Y. Takada, J. Superconductivity 18, 785 (2005).
[21] K. Takayanagi and E. Lipparini, Phys. Rev. B 56, 4872 (1997).
[22] Y. Takada, Phys. Rev. Lett. 87, 226402 (2001).
[23] J. C. Ward, Phys. Rev. 78, 182 (1950).
[24] Y. Takahashi, Nuovo Cimento 6, 371 (1957).
[25] Y. Takada and H. Yasuhara, Phys. Rev. Lett. 89, 216402 (2002).
[26] Y. Takada, Molecular Phys. 114, 1041 (2016).
[27] P. Gori-Giorgi and P. Ziesche, Phys. Rev. B 66, 235116 (2002).
[28] C. F. Richardson and N. W. Ashcroft, Phys. Rev. B 50, 8170
(1994).
[29] M. Lein, E. K. U. Gross, and J. P. Perdew, Phys. Rev. B 61,
13431 (2000).
[30] See, for example, Y. Takada, Phys. Rev. B 47, 3482 (1993).
[31] J. C. Kimball, J. Phys. A: Math. Gen. 8, 1513 (1975).
[32] H. Yasuhara, Solid State Commun. 11, 1481 (1972).
[33] H. Yasuhara, J. Phys. Soc. Jpn. 36, 361 (1974).
[34] A. W. Overhauser, Can. J. Phys. 73, 683 (1995).

[35] P. Gori-Giorgi and J. P. Perdew, Phys. Rev. B 64, 155102(2001).
[36] S. Moroni, D. M. Ceperley, and G. Senatore, Phys. Rev. Lett.
75, 689 (1995).
[37] G. Niklasson, Phys. Rev. B 10, 3052 (1974).
[38] I. O. Kulik, Zh. Eksp. Teor. Fiz. 40, 1343 (1961) [Sov. Phys.
JETP 13, 946 (1961)].
[39] E. K. U. Gross and W. Kohn, Phys. Rev. Lett. 55, 2850 (1985);
57, 923 (1986).
[40] B. Dabrowski, Phys. Rev. B 34, 4989 (1986).
[41] N. Iwamoto and E. K. U. Gross, Phys. Rev. B 35, 3003 (1987).
[42] A. Holas and K. S. Singwi, Phys. Rev. B 40, 158 (1989).
[43] R. Nifosı̀, S. Conti, and M. P. Tosi, Phys. Rev. B 58, 12758(1998).
[44] Z. Qian and G. Vignale, Phys. Rev. B 65, 235121 (2002).
[45] K. Morawetz, Phys. Rev. B 66, 075125 (2002).
[46] L. A. Constantin and J. M. Pitarke, Phys. Rev. B 75, 245127(2007).
[47] R. D. Puff, Phys. Rev. 137, A406 (1965).
[48] K. N. Pathak and P. Vashishta, Phys. Rev. B 7, 3649 (1973).
[49] N. Iwamoto, Phys. Rev. A 30, 3289 (1984).
[50] J. P. Eisenstein, L. N. Pfeiffer, and K. W. West, Phys. Rev. Lett.
68, 674 (1992); Phys. Rev. B 50, 1760 (1994).
[51] S. Ilani, A. Yacoby, D. Mahalu, and H. Shtrikman, Phys. Rev.
Lett. 84, 3133 (2000); Science 292, 1354 (2001).
[52] A. W. Overhauser, Phys. Rev. 128, 1437 (1962); 167, 691 (1968);
Adv. Phys. 27, 343 (1978).
[53] J. P. Perdew and T. Datta, Phys. Status Solidi (b) 102, 283 (1980);
H. T. Tran and J. P. Perdew, Am. J. Phys. 71, 1048 (2003).

245106-14