# Energy Band Structure of Multistream Quantum Electron System

M. Akbari-Moghanjoughi  Faculty of Sciences, Department of Physics, Azarbaijan
Shahid Madani University, 51745-406 Tabriz, Iran<EMAIL_ADDRESS>
###### Abstract

In this paper, using the quantum multistream model, we develop a method to
study the electronic band structure of plasmonic excitations in streaming
electron gas with arbitrary degree of degeneracy. The multifluid quantum
hydrodynamic model is used to obtain $N$-coupled pseudoforce differential
equation system from which the energy band structure of plasmonic excitations
is calculated. It is shown that inevitable appearance of energy bands
separated by gaps can be due to discrete velocity filaments and their
electrostatic mode coupling in the electron gas. Current model also provides
an alternative description of collisionless damping and phase mixing, i.e.,
collective scattering phenomenon within the energy band gaps due to mode
coupling between wave-like and particle-like oscillations. The quantum
multistream model is further generalized to include virtual streams which is
used to calculate the electronic band structure of one-dimensional plasmonic
crystals. It is remarked that, unlike the empty lattice approximation in free
electron model, energy band gaps exist in plasmon excitations due to the
collective electrostatic interactions between electrons. It is also shown that
the plasmonic band gap size at first Brillouin zone boundary maximizes at the
reciprocal lattice vector, $G$, close to metallic densities. Furthermore, the
electron-lattice binding and electron-phonon coupling strength effects on the
electronic band structure are discussed. It is remarked that inevitable
formation of energy band structure is a general characteristics of various
electromagnetically and gravitationally coupled quantum multistream systems.

###### pacs:

52.30.-q,71.10.Ca, 05.30.-d

##  I Introduction

Plasmons are high frequency elementary quantized excitations of electron
plasma oscillations [1, 2]. They play inevitable role in many fundamental
properties of plasmas semiconductors and metallic nanoparticles from electric
and heat transport phenomena to optical and dielectric response, etc. [3, 4].
Dynamics of these quantized electromagnetic quasiparticles make an ideal
platform for miniaturization of ultrafast terahertz device communications [5],
where conventional integrated circuits fail to operate. They also have
numerous other interesting applications in nanotechnology [6], plasmonics [7,
8, 9], optoelectronics [10], etc. for engineering low-dimensional nano-
fabricated semiconductor industry [11, 12, 13]. Energy conversion by plasmons
is an new way of solar power extraction due to its high efficiency in
photovoltaic and catalytic devices. Use of the collective oscillations of
electrons instead of single particles makes huge amount of energy extraction
in an operation step in plasmonic solar devices [14, 15].In local surface
plasmon resonance (LSPR) [16] process, the surface electrons, the so called
hot electrons, are collectively excited by electromagnetic radiations in UV-
VIS range generating a huge amount of energy transfer. The hot electron
current are collected in an appropriate contacts of nanoparticle surfaces by
an efficient electron collecting material like $TiO_{2}$ in Schottky
configuration [17].

Collective charge screening effect which manifests itself as the
characteristic optical edge in metallic surfaces already have may applications
in metallic alloys making then optically unique among other solids. Collective
electron excitations rule almost every aspect of solid from optical to
dielectric response in plasmas [18, 19, 20] and condensed matter. Recent
infrared spectroscopic techniques shows that Low dimensional semiconductors
[21] such as gapped graphene also demonstrate interesting surface plasmon
effects. The collective electron transport property of graphene makes it an
ideal element for multilayer composite devices such as compact ultrafast
switches, optical modulators, optical lattices, photodetectors, tandem solar
cells and biosensors [22, 23, 24]. The first theoretical development of the
idea of collective electron excitations by Bohm and Pines dates back to mid-
nineties, when they coined plasmon name for such excitations due to the long-
range electromagnetic nature of interactions [25, 26, 27, 28, 29]. The
theoretical as well as experimental aspects of collective electron dynamics in
quantum level has been the subject of intense investigations over the past few
decades [30, 31, 32, 33, 34, 35], due to its fundamental importance in many
field of physics and chemistry.

Pioneering developments of quantum statistical and kinetic theories [36, 37,
38, 39] had a long tradition furnishing a pavement for modern theories of
quantum plasmas [40, 41, 42, 43]. Many interesting new aspects of collective
quantum effects in astrophysical and laboratory plasmas has been recently
investigated using quantum plasma theories [44, 45, 46, 47, 48, 49, 50, 51,
52, 53, 54, 55, 56, 57]. The quantum kinetic theories like time-dependent
density functional theories (TDFT) are, however, less analytic as compared to
the quantum hydrodynamic analogues, due mostly to mathematical complexity
which require large scale computational programming. Recent investigation
reveals [58] that quantum hydrodynamic approaches based on the density
functional formalism [43] can reach beyond the previously thought kinetic
limitations, such as the collisionless damping if accurately formulated. One
of the most effective hydrodynamic formalism for studying the quantum aspects
of plasmas is the Schrödinger-Poisson model [59, 60], based on the Madelung
quantum fluid theory which originally attempted for the single-electron
quantum fluid modeling [61]. It has been recently shown that the analytic
investigation of linearized Schrödinger-Poisson system for arbitrary
degenerate electron gas provides routes to some novel quantum feature of
collective plasmon excitations [62, 63]. In current study we use the
multistream model in order to investigate the band structure plasmon
excitations in streaming plasmas and plasmonic lattices.

##  II Mathematical Model

Starting with a one dimensional collision-less multi-fluid quantum
hydrodynamic model for electron gas with an arbitrary degree of degeneracy,
the set of equations read

$\displaystyle\frac{{\partial{n_{s}}}}{{\partial
t}}+\frac{{\partial{n_{s}}{v_{s}}}}{{\partial x}}=0,$ (1a)
$\displaystyle\frac{{\partial{v_{s}}}}{{\partial
t}}+{v_{s}}\frac{{\partial{v_{s}}}}{{\partial
x}}=\frac{{{e}}}{{{m_{e}}}}\frac{{\partial\phi}}{{\partial
x}}-\frac{1}{{{m_{e}}}}\frac{{\partial{\mu}}}{{\partial
x}}+\frac{{{\hbar^{2}}}}{{2m_{e}^{2}}}\frac{\partial}{{\partial
x}}\left({\frac{1}{{\sqrt{{n_{s}}}}}\frac{{{\partial^{2}}\sqrt{{n_{s}}}}}{{\partial{x^{2}}}}}\right),$
(1b) $\displaystyle\frac{{{\partial^{2}}\phi}}{{\partial{x^{2}}}}=4\pi
e\sum\limits_{s}{{n_{s}}},$ (1c)

in which the dependent variables, $n_{s}$ and $v_{s}$ refer to the number
density and fluid velocity of given electron stream, indexed by $s$, and
$\phi$ is the electrostatic potential. The last term in momentum equation
arises due to quantum Bohm potential, by elimination of which the system (1)
reduces to the classical Dawson’s multistream model [64]. Moreover, $\mu$ is
the chemical potential of the electron gas which is related to the electron
number density using an appropriate equation of state (EoS) and is used to
close the hydrodynamic system (1). For isothermal electron gas of arbitrary
degeneracy the EoS is

$\displaystyle
n(\nu,T)=\frac{{{2^{7/2}}\pi{m_{e}^{3/2}}}}{{{h^{3}}}}{F_{1/2}}(\nu)=-\frac{{{2^{5/2}}{{(\pi
m_{e}{k_{B}}T)}^{3/2}}}}{{{h^{3}}}}{\rm{L}}{{\rm{i}}_{3/2}}[-\exp(\nu)],$ (2a)
$\displaystyle
P(\nu,T)=\frac{{{2^{9/2}}\pi{m_{e}^{3/2}}}}{{3{h^{3}}}}{F_{3/2}}(\nu)=-\frac{{{2^{5/2}}{{(\pi
m_{e}{k_{B}}T)}^{3/2}}({k_{B}}T)}}{{{h^{3}}}}{\rm{L}}{{\rm{i}}_{5/2}}[-\exp(\nu)],$
(2b)

in which $P$ is the statistical pressure of the gas and $\nu=\beta\mu$ with
$\beta=1/k_{B}T$, where, $F_{k}$ is the Fermi integral of order $k$ given as

${F_{k}}(\nu)=\int_{0}^{\infty}{\frac{{{x^{k}}}}{{\exp(x-\nu)+1}}}dx.$ (3)

In terms of polylog function, ${\rm{Li}}_{k}$, the Fermi integrals are defined
as

${F_{k}}(\nu)=-\Gamma(k+1){\rm{L}}{{\rm{i}}_{k+1}}[-\exp(\nu)],$ (4)

in which $\Gamma$ is the conventional gamma function. It is seen that, the
thermodynamic identity, $\partial P/\partial\mu=n$, holds for the electron gas
in the thermodynamic equilibrium. Note also that we ignore the chemical
potential dependence on stream index in the gas and it is assumed that the
index $s$ characterizes only the velocity spectrum in the system. For a
multispecies plasmas, however, this index may apply to the chemical potential
of species. The hydrodynamic model (1) may be cast into a more simple form as
the effective Schrödinger-Poisson model [59] for the system, using the
Madelung transformations ${\cal
N}_{s}(x,t)=\sqrt{n_{s}(x,t)}\exp[iS_{s}(x,t)/\hbar]$ and
$v_{s}(x,t)=(1/m_{e})\partial S_{s}(x,t)/\partial x$. By using the later
definition, the continuity and momentum balance, after separation of
real/imaginary parts become

$\displaystyle{m_{e}}\frac{{\partial{n_{s}}}}{{\partial
t}}+\frac{{\partial{n_{s}}}}{{\partial x}}\frac{{\partial{S_{s}}}}{{\partial
x}}+{n_{s}}\frac{{{\partial^{2}}{S_{s}}}}{{\partial{x^{2}}}}=0,$ (5a)
$\displaystyle\frac{{{\partial^{2}}{S_{s}}}}{{\partial t\partial
x}}+\frac{1}{m_{s}}\frac{{\partial{S_{s}}}}{{\partial
x}}\frac{{{\partial^{2}}{S_{s}}}}{{\partial{x^{2}}}}=\frac{{{e}\partial\phi}}{{\partial
x}}-\frac{{\partial{\mu}}}{{\partial x}}+\frac{{\partial{B_{s}}}}{{\partial
x}},$ (5b)
$\displaystyle{B_{s}}=\frac{{{\hbar^{2}}}}{{8{m_{e}}n_{s}^{2}}}\left[{2{n_{s}}\frac{{{\partial^{2}}{n_{s}}}}{{\partial{x^{2}}}}-{{\left({\frac{{\partial{n_{s}}}}{{\partial
x}}}\right)}^{2}}}\right],$ (5c)

which combining with ${\cal
N}_{s}(x,t)=\sqrt{n_{s}(x,t)}\exp[iS_{s}(x,t)/\hbar]$ together with the
Poisson’s equation leads to the following Schrödinger-Poisson system [65]

$\displaystyle i\hbar{\frac{{\partial{\cal N}_{s}}}{{\partial
t}}=-\frac{{{\hbar^{2}}}}{{2{m_{e}}}}\frac{{{\partial^{2}}{\cal
N}_{s}}}{{\partial{x^{2}}}}-e\phi{\cal N}_{s}+{\mu}{\cal N}_{s}},$ (6a)
$\displaystyle\frac{{{\partial^{2}}\phi}}{{\partial{x^{2}}}}=4\pi
e\sum\limits_{s}{|{\cal N}_{s}{|^{2}}}.$ (6b)

For our purpose let us consider a particular case of multistream model in the
linear perturbation limit in which every stream of electrons is a
monoenergetic beam interacting with others through the electrostatic
potential. One may linearize the system (7) through the assumptions
$n_{s}(x,t)=n_{0}+n_{1s}$, $\phi_{s}(x,t)=0+\phi_{1}$, $S_{s}(x)=S_{0}+p_{s}x$
and $\mu=\mu_{0}$ with $\mu_{0}$ being the equilibrium chemical potential of
the gas and $p_{s}$ being the constant momentum of given stream. Also, $S_{0}$
results in a constant phase which is ignored in this analysis. In the steady
state limit we are also able to decompose the state function into
spatiotemporal product of variables as, ${\cal
N}_{s}(x,t)=\psi_{s}(t)\psi_{s}(x)\exp(ik_{s}x)$, where,
$\psi_{s}(x)=\sqrt{n_{s}(x)}$ with $k_{s}=p_{s}/\hbar$ being the de Broglie
wavenumber of given stream. Hence, eliminating the first-order index, the
normalized linear system of coupled equations read

$\displaystyle\frac{{{d^{2}}{\Psi_{s}(x)}}}{{d{x^{2}}}}+2i{k_{s}}\frac{{d{\Psi_{s}(x)}}}{{dx}}+\Phi(x)+(E-k_{s}^{2}){\Psi_{s}(x)}=0,$
(7a)
$\displaystyle\frac{{{d^{2}}\Phi(x)}}{{d{x^{2}}}}-\sum\limits_{s}f_{s}\Psi_{s}(x)=0,$
(7b)

where $E=(\epsilon-\mu_{0})/E_{p}$ is the normalized multistream system energy
with $E_{p}=\hbar\omega_{p}$ being the plasmon energy and
$\omega_{p}=\sqrt{4\pi e^{2}n_{0}/m_{e}}$ the electron plasma frequency. Also,
$f_{s}$ represents the momentum ($\hbar k_{s}$) distribution function for
given stream $s$ with the property $\sum\limits_{s}f_{s}=1$. However in
current multistream model we have chosen $f_{s}=1$, for simplicity. Moreover,
$\epsilon=\hbar\omega$ defines the common energy eigenvalue of the multistream
system. Moreover, we used the normalization scheme,
$\Psi_{s}(x)\to\Psi_{s}(x)/\sqrt{n_{0}}$, with $n_{0}$ being the equilibrium
electron number density and $\Phi(x)=e\phi/E_{p}$. The space coordinate $x$ is
normalized to the plasmon length, $l_{p}=2\pi/k_{p}$, with
$k_{p}=\sqrt{2m_{e}E_{p}}/\hbar$ being the characteristic plasmon wavenumber.
Therefore, the de Broglie wavenumber is normalized to the plasmon length and
temperature to the plasmon temperature $T_{p}=E_{p}/k_{B}$. The system (7),
plus the temporal term proportional to $\exp(-i\omega t)$, describes the
steady state evolution of an electron gas in the linear limit. To obtain the
state functions ${\cal N}(x,t)=\sum\limits_{s}\Psi_{s}(x)\exp(-i\Omega t)$
($\Omega=\omega/\omega_{p}$) and $\Phi(x)$ one has to evaluate $N$-coupled
differential equations through the electrostatic potential each given for an
electron stream. Note that $\Psi_{s}(x)$ characterize the pure states of the
multistream system from which the mixed states are calculated. The fluid
velocity of each stream satisfy the relation
$v_{s}(x,t)=j_{s}(x,t)/n_{s}(x,t)$, where,
$j_{s}(x,t)=i\hbar/(2m_{e})[\partial{\cal N}_{s}(x,t)/\partial x\times{\cal
N}^{*}_{s}(x,t)-\partial{\cal N}^{*}_{s}(x,t)/\partial x\times{\cal
N}_{s}(x,t)]$ is the current density of given stream. This velocity is also
given by the relation $v_{s}(x,t)=(\hbar/m_{e}){\rm{Im}}[\partial{\cal
N}_{s}(x,t)/\partial x/{\cal N}_{s}(x,t)]$, which is identical with the
pseudoparticle velocity in the pilot-wave theory for guiding equation. In this
linear limit we have $v_{s}=\hbar k_{s}/m_{e}$. The multistream velocity is
obtained through the state function as
$v(x,t)=(\hbar/m_{e}){\rm{Im}}[\partial{\cal N}(x,t)/\partial x/{\cal
N}(x,t)]$ by solving the $N$-coupled equations (7).

Figure 1: (a) Variation in the plasmon energy $E_{p}=\hbar\omega_{p}$ in terms
of electron number density in logarithmic scale. (b) Variations in plasmon
length $1/k_{p}$ with electron number density in logarithmic scale. (c) The
electron concentration in terms of equilibrium chemical potential of arbitrary
degenerate electron gas. (d) Variation of plasmon temperature
$T_{p}=E_{p}/k_{B}$ in terms of the chemical potential of the electron gas for
different values of the electron temperature. The increase in thickness of
curves in plot (d) indicates the increase in varied parameter above the panel.

Figure 1 shows the variations in characteristic parameters of the plasmon
system with electron number density and an arbitrary degenerate electron gas.
In Fig. 1(a) it is shown that the plasmon energy varies up to few
electronvolts from the classical to fully degenerate electron gas. Note that
only in highly doped semiconductors and metals this energy becomes
significant. The plasmon length variations in nanometer unit in terms of
electron number density is shown in Fig. 1(b). This length sharply decreases
with increase in number density to a fraction of a nanometer in a typical
metal. The chemical potential variation is shown in Fig. 1(c). In the fully
degenerate (zero temperature) limit in typical metals the chemical potential
at $E=0$ or $\epsilon=\mu$ characterizes the fundamental Fermi energy level of
the system which is assumed to be constant. Hence, current model is most
appropriate for metals and nano-metallic density regime and beyond. Figure
1(d), on the other hand, shows the variation plasmon temperature, $T_{p}$ with
the variations in chemical potential for different values of the electron
temperature. It is remarked that, the plasmon temperature increases with
increase in the chemical potential but decreases with increase in electron
temperature.

##  III One-Stream Model and the Doppler Shift

Despite the simplicity of the model (7), it will be shown that it is useful in
describing some fundamental physical phenomenon corresponding to the plasmon
system. Consider a single stream ($s=1$) described by the following system

$\displaystyle\frac{{{d^{2}}{\Psi_{1}(x)}}}{{d{x^{2}}}}+2i{k_{1}}\frac{{d{\Psi_{1}(x)}}}{{dx}}+\Phi(x)+(E-k_{1}^{2}){\Psi_{1}(x)}=0,$
(8a) $\displaystyle\frac{{{d^{2}}\Phi(x)}}{{d{x^{2}}}}-\Psi_{1}(x)=0,$ (8b)

satisfying the energy dispersion relation
${E}={1}/{{{k^{2}}}}+{\left({k+{k_{1}}}\right)^{2}}$ which reduces to the
plasmon dispersion relation in the limit $k_{1}=0$. It is remarked that the
particle branch of the energy dispersion is Doppler shifted due to the
streaming electrons. However, the wave-like branch is not affected by the
electron drift. It has been shown that in an inertial frame moving along with
the electron beam the traveling wave solution to the system (8) becomes
identical with that of the electron gas in rest frame [68] with the beam speed
replacing the energy eigenvalues.

Figure 2: Dispersion curves of (a) static and (b) streaming free electron
(dashed curve) and plasmon (solid curve) excitations. (c) Variation of state-
functions $\Psi(x)$ (thin curve) and $\phi(x)$ (thick curve) in streaming
electron gas with the de Broglie wavenumber $k_{1}=0.5$ at stable orbital
$E=3.5$. (d) Quasiparticle velocity corresponding to the one-stream electron
gas in (c) at stable $E=3.5$ (thick curve) and unstable $E=2$ (thin curve)
orbital.

Figure 2 shows the energy dispersion and state functions of plasmon
excitations for given parameters. The energy dispersion plasmon excitations
for the case of $k_{1}=0$ (solid curves) along with the free electron
dispersion (dashed curve) are shown in Fig. 2(a). There are stable plasmon
excitations above the critical value $E=2$ (as shown by horizontal line) which
are double-tone due to both particle-like ($k>1$) and wave-like ($k<1$)
phenomena. However, below this critical line the excitation wavenumbers become
complex and energy exchange occurs between the partcile-like and wave-like
branches, as discussed in Ref. [68]. Figure 2(b) depicts the energy dispersion
of excitations for $k_{1}=0.5$. It is remarked that the free electron
dispersion undergoes a Doppler shift and two critical minimum values for
energy appear, namely, $E_{m1}\simeq 1.1786$ and $E_{m2}\simeq 3.1944$. For
$E>E_{m2}$ the plasmon excitation with four real wavenumbers are stable.
However, for $E_{m1}<E<E_{m2}$ only two of the excitations wavenumbers are
real, hence, excitations are unstable. For $E<E_{m1}$ all four wavenumbers
become complex and plasmon excitations become unstable again. However, there
is a fundamental difference between the two unstable regimes $E_{m1}<E<E_{m2}$
and $E<E_{m1}$, as will be discussed later. Figure 2(c) shows the profiles of
state functions, namely $\Psi(x)$ (thin curve) and $\Phi(x)$ (thick curve),
for given stable oscillation parameter values. These state function have been
obtained by numerical solution of (8) with initial conditions,
$\Phi(0)=\Phi^{\prime}(0)=\Psi^{\prime}(0)=0$ and $\Psi(0)=1$. The variation
of quasiparticle velocity for stable orbital ($E=3.5$ as thick curve) and
unstable orbital ($E=2$ as thin curve) is depicted in Fig. 2(d). Evidently,
there are oscillations in the velocity profiles which are damped for unstable
energy orbital $E=2$.

##  IV Two-Stream Model and Energy Band Formation

Let us consider the following symmetric two stream system

$\displaystyle\frac{{{d^{2}}{\Psi_{1}}}}{{d{x^{2}}}}+2i{k_{1}}\frac{{d{\Psi_{1}}}}{{dx}}+\Phi+(E-k_{1}^{2}){\Psi_{1}}=0$
(9a)
$\displaystyle\frac{{{d^{2}}{\Psi_{2}}}}{{d{x^{2}}}}+2i{k_{2}}\frac{{d{\Psi_{2}}}}{{dx}}+\Phi+(E-k_{2}^{2}){\Psi_{2}}=0$
(9b) $\displaystyle\frac{{{d^{2}}\Phi}}{{d{x^{2}}}}-\Psi_{1}-\Psi_{2}=0,$ (9c)

where $k_{1}$ and $k_{2}$ are de Broglie wavenumbers of the streams. The
energy dispersion relation can be obtained by assuming plane-wave expansions,
$\Psi_{1}(x)=\Psi_{11}\exp(ikx)$, $\Psi_{2}(x)=\Psi_{21}\exp(ikx)$ and
$\Phi_{1}(x)=\Phi_{1}\exp(ikx)$, which leads to the following eigenvalue
system

$\left({\begin{array}[]{*{20}{c}}{E-{{(k+{k_{1}})}^{2}}}&0&1\\\
0&{E-{{(k+{k_{2}})}^{2}}}&1\\\
1&1&{{k^{2}}}\end{array}}\right)\left({\begin{array}[]{*{20}{c}}{{\Psi_{11}}}\\\
{{\Psi_{21}}}\\\
{{\Phi_{1}}}\end{array}}\right)=\left({\begin{array}[]{*{20}{c}}0\\\ 0\\\
0\end{array}}\right).$ (10)

The energy dispersion relation of two-stream system reads

$E_{\pm}={k^{2}}+\frac{1}{{{k^{2}}}}+\frac{1}{2}\left({k_{1}^{2}+k_{2}^{2}}\right)+k\left({{k_{1}}+{k_{2}}}\right)\pm\frac{{\sqrt{4+{k^{4}}{{\left({{k_{1}}+{k_{2}}}\right)}^{2}}{{\left({{k_{1}}+{k_{2}}+2k}\right)}^{2}}}}}{{2{k^{2}}}}.$
(11)

Figure 3: (a) Energy band structure of symmetric two-stream plasmon
excitations. (b) Energy band structure of asymmetric two-stream plasmon
excitations. Quasiparticle velocity of two-stream excitations in symmetric
band gap. (d) Quasiparticle velocity of two-stream excitations in asymmetric
band gap showing the phase mixing effect due to mode coupling of the energy
bands.

It is remarked that, in the two-stream model extra energy band appears due to
mode coupling between the streams. Figure 3 shows the structure of energy
bands in two-stream model. In Fig. 3(a) the two streams have same but opposite
velocity and the upper and lower energy bands are symmetric. The imbalanced
two stream is shown in Fig. 3(b), the asymmetry of which is caused by the
Doppler shift in particle-like branch. The quasiparticle orbital velocities in
each stream are shown in Figs. 3(c) and 3(d) for symmetric and asymmetric
cases. It is seen that in the asymmetric two stream model the orbital velocity
of quasiparticle in electron beam with $k_{1}=1$ reverses at $E=4$ and merges
with the other stream. The later phenomenon, which we may call the phase
mixing effect, is a novel feature of the two-stream model caused by the
collective wave-particle interactions in the energy band gaps. Note that in
the preceding analysis (and the following) we consider equal density
distribution for streams, for the sake of simplicity. In fact a generalized
momentum distribution function, such as the Maxwell-Boltzmann, may be used in
the Poisson’s equation as weight function of stream probability functions
$\Psi_{N}$. However, in the limit of full degeneracy the momentum distribution
becomes unity.

##  V Energy Band Structure in Multistream Model

To this end, it is straightforward to generalize the model to include a large
number of streams each characterized by their de Broglie’s wavenumber,
$k_{N}$. Therefore, $N$-coupled differential equation system read

$\displaystyle\frac{{{d^{2}}{\Psi_{1}}}}{{d{x^{2}}}}+2i{k_{1}}\frac{{d{\Psi_{1}}}}{{dx}}+\Phi+(E-k_{1}^{2}){\Psi_{1}}=0,$
(12a)
$\displaystyle\frac{{{d^{2}}{\Psi_{2}}}}{{d{x^{2}}}}+2i{k_{2}}\frac{{d{\Psi_{2}}}}{{dx}}+\Phi+(E-k_{2}^{2}){\Psi_{2}}=0,$
(12b) $\displaystyle\hskip 28.45274pt\vdots\hskip 56.9055pt\vdots\hskip
56.9055pt\vdots\hskip 56.9055pt$ (12c)
$\displaystyle\frac{{{d^{2}}{\Psi_{N}}}}{{d{x^{2}}}}+2i{k_{N}}\frac{{d{\Psi_{N}}}}{{dx}}+\Phi+(E-k_{N}^{2}){\Psi_{N}}=0,$
(12d)
$\displaystyle\frac{{{d^{2}}\Phi}}{{d{x^{2}}}}-\Psi_{1}-\Psi_{2}-\cdots-\Psi_{N}=0,$
(12e)

Fourier analysis of (12) leads to the following eigenvalue system

$\left({\begin{array}[]{*{20}{c}}{E-{{(k+{k_{1}})}^{2}}}&0&\ldots&{\begin{array}[]{*{20}{c}}{\begin{array}[]{*{20}{c}}{}\hfil&{}\hfil&0&{\begin{array}[]{*{20}{c}}{}\hfil&{}\hfil\end{array}}\end{array}}&1\end{array}}\\\
0&{E-{{(k+{k_{2}})}^{2}}}&\ldots&{\begin{array}[]{*{20}{c}}{\begin{array}[]{*{20}{c}}{}\hfil&{}\hfil&0&{\begin{array}[]{*{20}{c}}{}\hfil&{}\hfil\end{array}}\end{array}}&1\end{array}}\\\
\vdots&\vdots&\ddots&{\begin{array}[]{*{20}{c}}{\begin{array}[]{*{20}{c}}{}\hfil&{}\hfil&{\begin{array}[]{*{20}{c}}\vdots&{}\hfil\end{array}}&{}\hfil\end{array}}&\vdots\end{array}}\\\
{\begin{array}[]{*{20}{c}}0\\\
1\end{array}}&{\begin{array}[]{*{20}{c}}\ldots\\\
1\end{array}}&{\begin{array}[]{*{20}{c}}0\\\
\cdots\end{array}}&{\begin{array}[]{*{20}{c}}{\begin{array}[]{*{20}{c}}{E-{{(k+{k_{N}})}^{2}}}&1\end{array}}\\\
{\begin{array}[]{*{20}{c}}{}\hfil&{\begin{array}[]{*{20}{c}}{}\hfil&1&{\begin{array}[]{*{20}{c}}{}\hfil&{\begin{array}[]{*{20}{c}}{}\hfil&{{k^{2}}}\end{array}}\end{array}}\end{array}}\end{array}}\end{array}}\end{array}}\right)\left({\begin{array}[]{*{20}{c}}{{\Psi_{11}}}\\\
{\begin{array}[]{*{20}{c}}{{\Psi_{21}}}\\\ \vdots\\\
{{\Psi_{N1}}}\end{array}}\\\
{{\Phi_{1}}}\end{array}}\right)=\left({\begin{array}[]{*{20}{c}}0\\\
{\begin{array}[]{*{20}{c}}0\\\ \vdots\\\ 0\end{array}}\\\
0\end{array}}\right).$ (13)

The system (13) can be evaluated numerically to any number of streams in order
to calculate the structure of energy bands in multistream electron system.

Figure 4: (a) The periodic $N$-coupled ($N=10$), $k_{N}=Nk_{1}$, energy band
structure of multistream electron gas with $k_{1}=1$. (b) The periodic
10-coupled, $k_{N}=Nk_{1}$, energy band structure of multistream electron gas
with $k_{1}=2$. (c) The random $N$-coupled ($N=10$), energy band structure of
multistream electron gas with $k_{1}=1$. (d) The random 10-coupled energy band
structure of multistream electron gas with $k_{1}=2$.

Figure 4 shows the electronic energy band structure of a 10-coupled ($N=10$)
multistream system. The periodic stream profile in Fig. 4(a) corresponds to
the especial case of $k_{N}=Nk_{1}$ with $N=10$ and $k_{1}=1$ illustrating
different bands which are separated via band gaps. As mentioned previously the
origin of band gaps in this case is the mode coupling between different
electron streams, very similar to the electronic band structure in crystalline
materials, as will be shown in the next section. Note that band gaps increase
with increase in energy, $E$, but decrease with increase in wavenumber, $k$.
The periodic 10-coupled band structure for $k_{1}=2$ is depicted in Fig. 4(b).
It is remarked that increase in the stream speed leads to increase in the
energy band widths. Figures 4(c) and 4(d) depict the band structure of
10-coupled randomly distributed multistream system in the arbitrary degenerate
electron gas. The increase in the number of streams and arbitrary velocity
distribution can lead to formation of a very complex energy band structure
which describes the plasmon excitations in the multistream system.

We have already seen that presence of discrete electron streams leads to
formation of energy band structure and energy gaps leading to complex
wavenumbers with the imaginary part representing the growing/damping features.
It has been shown [67] that in collisionless electron systems the complex
wavenumbers does not lead to dissipation of energy but the exchange between
the particle-like and wave-like excitations very similar to the Landau damping
phenomenon. In extreme limit of our multistream theory when every electron
(velocity) specie constitute an individual stream with the weigh function,
$f_{s}$, being the Maxwell-Boltzmann velocity distribution, large number of
energy minibands form which are separated by tiny band gaps through which
phase-mixing, i.e. quasiparticle scattering, can occur. In this case energy
exchanges can take place between the collective electrostatic excitations and
single electron oscillations by irreversible energy transfer from wave-like
oscillations to the particle-like ones. the multistream model has been
originally used by Dawson in order to give a physical interpretation of the
Landau Damping effect [64]. Indeed the effect can occur by resonant electrons
which have speeds close to the plasmon excitation phase speed, i.e., $v\simeq
E/\hbar k$, residing in a band gap.

To better understand the collisionless damping phenomenon, we consider the
traveling wave solutions of a electron beam (stream) with normalized drift
speed, $\gamma$, $\Phi(x-2\gamma t)$ and ${\cal N}(x,t)=\Psi(x-2\gamma
t)\exp[i\gamma(x-2\gamma t)]$ obtained in Ref. [68], in which

$\left\\{{\begin{array}[]{*{20}{c}}{{\Phi}(x-2\gamma t)}\\\ {{\Psi}(x-2\gamma
t)}\end{array}}\right\\}=\frac{1}{{2\alpha}}\left\\{{\begin{array}[]{*{20}{c}}{{\Psi_{0}}+k_{p}^{2}{\Phi_{0}}}&{-\left({{\Psi_{0}}+k_{w}^{2}{\Phi_{0}}}\right)}\\\
{-\left({{\Phi_{0}}+k_{w}^{2}{\Psi_{0}}}\right)}&{{\Phi_{0}}+k_{p}^{2}{\Psi_{0}}}\end{array}}\right\\}\left\\{{\begin{array}[]{*{20}{c}}{\cos[{k_{w}}(x-2\gamma
t)]}\\\ {\cos[{k_{p}}(x-2\gamma t)]}\end{array}}\right\\},$ (14)

in which $\Phi_{0}$ and $\Psi_{0}$ define the initial values at $x=2\gamma t$
and the characteristic wave-like and particle-like wavenumbers $k_{w}$ and
$k_{p}$ are given,respectively, as

${k_{w}}=\sqrt{(E_{d}-\alpha)/2},\hskip
8.53581pt{k_{p}}=\sqrt{(E_{d}+\alpha)/2},\hskip
8.53581pt\alpha=\sqrt{{E_{d}^{2}}-4},$ (15)

where $E_{d}=\gamma^{2}-\mu$ and $k_{w}k_{p}=1$. Inside the gap the
wavenumbers become complex, i.e., $k=k_{r}+ik_{i}$ with $k_{r}$ and $k_{i}$
being the real and imaginary parts of the wavenumbers. It has been shown that
[68], while the real parts are equal ($k_{wr}=k_{pr}$), the imaginary part of
wave-like excitation is always negative, ($k_{wi}<0$) and that of the
particle-like is always positive, ($k_{pi}=-k_{wi}>0$) for space-time range
$x>2\gamma t$. Therefore, the wave-like/particle-like oscillations grow in
space for energy values with imaginary de Broglie wavenumbers (inside the
energy gaps) where the electron streams experience the so-called quantum drift
instability [68]. On the other hand, a close inspection of the solution (14)
reveals that the wave-like/particle-like oscillations experience damping in
time, simultaneously. Generally speaking, particle-like/wave-like excitations
of arbitrary degenerate electron beam always undergo spacial/temporal
growing/damping inside the energy gaps. The above description of multistream
electron behavior may be regarded as an elegant quantum description of the
colissionless Landau damping effect due to the resonant wave-particle
interactions and energy exchange between wave-like and particle-like
oscillations within the energy band gaps. It is remarkable however that the
stream velocity defined through, $v_{s}={\mathop{\rm Im}\nolimits}[{\cal
N}_{x}(x,t)/{{\cal N}}(x,t)]=\gamma$, is invariant under the wave-particle
processes, indicating total energy conservation. It can be shown that
spacial/temporal growth/damping of wave-like/particle-like behavior is
intrinsic behavior of an electron beam in quantum tunneling process where
quantum drift instability takes place (similar to instability of electron
stream excitations within the energy band gaps). The detailed study of the
relationship between the wave-particle phenomenon and the collective tunneling
through a potential barrier may be the subject of a future study in the
framework of complex energy band structure and is beyond the scope of current
research.

##  VI Band Structure of 1D Plasmonic Crystals

In this section we would like to generalize the theory of multistream model to
band structure of plasmon excitations in periodic system like such as
plasmonic crystals [66]. Considering a lattice of constant $a$ the crystal is
characterized by reciprocal lattice vectors $G_{N}=NG_{1}$ with $G_{1}=2\pi/a$
being the first reciprocal lattice vector and $N$ is an integer number. Now,
we model the electronic excitations through the $N$-coupled virtual streams
(16) in which the reciprocal lattice wavevectors, $G_{N}$, play the role of de
Broglie’s wavenumber of virtual streams. Therefore, we have the following
$N$-coupled virtual stream system with a solution of type ${\cal
N}_{N}(x,t)=\Psi_{N}(x)\exp(iG_{N}x-i\Omega t+i\Theta_{N})$ in which $G_{N}/2$
characterize the N-th Brillouin zone boundary and $\Theta_{N}$ is the
arbitrary phase angle of the given stream.

$\displaystyle\frac{{{d^{2}}{\Psi_{1}}}}{{d{x^{2}}}}+2iG_{1}\frac{{d{\Psi_{1}}}}{{dx}}+\Phi+(E-G_{1}^{2}){\Psi_{1}}=0,$
(16a)
$\displaystyle\frac{{{d^{2}}{\Psi_{2}}}}{{d{x^{2}}}}+2iG_{2}\frac{{d{\Psi_{2}}}}{{dx}}+\Phi+(E-G_{2}^{2}){\Psi_{2}}=0,$
(16b) $\displaystyle\hskip 28.45274pt\vdots\hskip 56.9055pt\vdots\hskip
56.9055pt\vdots\hskip 56.9055pt$ (16c)
$\displaystyle\frac{{{d^{2}}{\Psi_{N}}}}{{d{x^{2}}}}+2i{G_{N}}\frac{{d{\Psi_{N}}}}{{dx}}+\Phi+(E-G_{N}^{2}){\Psi_{N}}=0,$
(16d)
$\displaystyle\frac{{{d^{2}}\Phi}}{{d{x^{2}}}}-\Psi_{1}-\Psi_{2}-\cdots-\Psi_{N}=0,$
(16e)

the Fourier analysis of (16) leads to the following eigenvalue system

$\left({\begin{array}[]{*{20}{c}}{E-{{(k+{G_{1}})}^{2}}}&0&\ldots&{\begin{array}[]{*{20}{c}}{\begin{array}[]{*{20}{c}}{}\hfil&{}\hfil&0&{\begin{array}[]{*{20}{c}}{}\hfil&{}\hfil\end{array}}\end{array}}&1\end{array}}\\\
0&{E-{{(k+{G_{2}})}^{2}}}&\ldots&{\begin{array}[]{*{20}{c}}{\begin{array}[]{*{20}{c}}{}\hfil&{}\hfil&0&{\begin{array}[]{*{20}{c}}{}\hfil&{}\hfil\end{array}}\end{array}}&1\end{array}}\\\
\vdots&\vdots&\ddots&{\begin{array}[]{*{20}{c}}{\begin{array}[]{*{20}{c}}{}\hfil&{}\hfil&{\begin{array}[]{*{20}{c}}\vdots&{}\hfil\end{array}}&{}\hfil\end{array}}&\vdots\end{array}}\\\
{\begin{array}[]{*{20}{c}}0\\\
1\end{array}}&{\begin{array}[]{*{20}{c}}\ldots\\\
1\end{array}}&{\begin{array}[]{*{20}{c}}0\\\
\cdots\end{array}}&{\begin{array}[]{*{20}{c}}{\begin{array}[]{*{20}{c}}{E-{{(k+{G_{N}})}^{2}}}&1\end{array}}\\\
{\begin{array}[]{*{20}{c}}{}\hfil&{\begin{array}[]{*{20}{c}}{}\hfil&1&{\begin{array}[]{*{20}{c}}{}\hfil&{\begin{array}[]{*{20}{c}}{}\hfil&{{k^{2}}}\end{array}}\end{array}}\end{array}}\end{array}}\end{array}}\end{array}}\right)\left({\begin{array}[]{*{20}{c}}{{\Psi_{11}}}\\\
{\begin{array}[]{*{20}{c}}{{\Psi_{21}}}\\\ \vdots\\\
{{\Psi_{N1}}}\end{array}}\\\
{{\Phi_{1}}}\end{array}}\right)=\left({\begin{array}[]{*{20}{c}}0\\\
{\begin{array}[]{*{20}{c}}0\\\ \vdots\\\ 0\end{array}}\\\
0\end{array}}\right).$ (17)

Note that we assumed an empty lattice approximation in which the atomic
lattice potential is negligible compared to the plasmon energy. Such
assumption can be valid in fully degenerate regime due to effective charge
screening or in the case on weak-potential plasmonic lattices. The system (17)
may be evaluated numerically for a finite number of Brillouin zone
approximation. Analytical solution to energy band dispersion exists for
$3$-coupled system with lattice momentum $-G$, $0$ and $+G$. For $N=3$
approximation-order, one obtains

$\displaystyle{E_{1}}=\frac{1}{{3{k^{2}}}}\left[{3+2{G^{2}}{k^{2}}+3{k^{4}}-\frac{\tau}{\delta}-\delta}\right],$
(18a)
$\displaystyle{E_{2}}=\frac{1}{{6{k^{2}}}}\left[{2\left({3+2{G^{2}}{k^{2}}+3{k^{4}}}\right)+\frac{\tau}{\delta}\left({1+{\rm{i}}\sqrt{3}}\right)+\delta\left({1-{\rm{i}}\sqrt{3}}\right)}\right],$
(18b)
$\displaystyle{E_{3}}=\frac{1}{{6{k^{2}}}}\left[{2\left({3+2{G^{2}}{k^{2}}+3{k^{4}}}\right)+\frac{\tau}{\delta}\left({1-{\rm{i}}\sqrt{3}}\right)+\delta\left({1+{\rm{i}}\sqrt{3}}\right)}\right],$
(18c)
$\displaystyle\delta={\left[{{G^{6}}{k^{6}}-36{G^{4}}{k^{8}}+\frac{1}{2}\sqrt{{{\left({54-2{G^{6}}{k^{6}}+72{G^{4}}{k^{8}}}\right)}^{2}}-4{{\tau}^{3}}}}-27\right]^{1/3}},$
(18d) $\displaystyle\tau={9+{G^{4}}{k^{4}}+12{G^{2}}{k^{6}}}.$ (18e)

Figure 5: (a) The electronic band structure in 1D plasmonic lattice with $G=2$
in three-zone ($N=3$) empty-lattice approximation showing direct and indirect
band gaps. (b) The electronic band structure in 1D plasmonic lattice with
$G=3$ in three-zone ($N=3$) empty-lattice approximation showing the band gap
at first Brillouin zone boundary, $k=G/2$. (c) Variations of $\Delta E_{0}$
and $\Delta E_{d}$ in terms of the reciprocal lattice vector, $G=2\pi/a$,
where, $a$ is the lattice constant. (d) Variations of $\Delta E_{id}$, $\Delta
E_{12}$ and the first conduction band height $\Delta E_{01}$ in terms of the
reciprocal lattice vector, $G=2\pi/a$.

Figure 5 depicts the results of calculation for $N=3$ plasmonic crystal
lattice dispersion in empty lattice approximation. It is remarkable that the
band gap is still present in the absence of lattice potential due to the
collective electron interactions contrary to the free electron model of solids
[3]. Figure 5(a) shows the energy band structure with $G=2$. Direct and
indirect band gaps are evident which are the result of mode coupling between
different virtual streams. There are in fact three band the upper one not
shown in the figure. Figure 5(b) shows the band structure for $G=3$ and the
band gap $\Delta E_{12}$ between the energy bands $E_{1}$ and $E_{2}$ at the
first Brillouin zone, $k=G/2$. In terms of reciprocal lattice vector $G$ one
obtains $\Delta E_{d}=3G^{2}/2$ for the direct gap at $k=0$, and

$\displaystyle\Delta
E_{id}=\frac{1}{{3{G^{2}}}}\left[{12+\frac{{7{G^{4}}}}{4}+\frac{{\left({1-{\rm{i}}\sqrt{3}}\right)\left({36+{G^{8}}}\right)}}{\zeta}+\left({1+{\rm{i}}\sqrt{3}}\right)\zeta}\right],$
(19a) $\displaystyle\Delta
E_{12}=\frac{{108-36{\rm{i}}\sqrt{3}+\left({3-{\rm{i}}\sqrt{3}}\right){G^{8}}+3{\zeta^{2}}+{\rm{i}}\sqrt{3}{\zeta^{2}}}}{{3{G^{2}}\zeta}},$
(19b)
$\displaystyle\zeta={\left({6{\rm{i}}{G^{4}}\sqrt{3}\sqrt{36-4{G^{4}}+{G^{8}}}-{G^{12}}}-216\right)^{1/3}},$
(19c)

for the corresponding indirect gap and direct gap at first Brillouin zone, in
three zone, $N=3$, approximation. Interesting features appear for variations
of these gap with the reciprocal lattice vector (lattice constant). Figure
5(c) shows that $\Delta E_{0}$ increases monotonically with increase in $G$
and consequently decrease in lattice constant $a$. However, with increase of
$G$ the value of direct gap at $k=0$ first increases and reaches a maximum
value at $G\simeq 1.7726$ and then passes through a minimum value at $G\simeq
2.17572$. Figure 5(d) depicts variations in the indirect gap and the direct
gap size at first Brillouin zone boundary. It is remarked that the indirect
gap maximizes at $G\simeq 1.52319$ and closes at $G\simeq 3.11302$. The value
of $\Delta E_{12}$ maximizes at $G\simeq 1.56508$. The first plasmon
conduction band $\Delta E_{01}$ which occurs at $k=G$ is an important gap for
fully degenerate electron gas (at zero temperature), where all the electrons
are packed under the Fermi level ($E=0$ or $\epsilon=\mu_{0}$). The variation
of this gap in terms of the reciprocal lattice vector is shown also in Fig.
5(d). It is noticed that the gap maximized at the value $G\simeq 1.20944$.
Analytical expression for this gap for $N=3$ is

$\displaystyle\Delta{E_{01}}=\frac{1}{{6{G^{2}}}}\left[{6+10{G^{4}}-\frac{{2\left({9+13{G^{8}}}\right)}}{\delta}-2\delta}\right],$
(20a)
$\displaystyle\delta={\left[{\frac{1}{2}\sqrt{{{\left({54+70{G^{12}}}\right)}^{2}}-4{{\left({9+13{G^{8}}}\right)}^{3}}}-35{G^{12}}}-27\right]^{1/3}}.$
(20b)

Figure 6: Variations in imaginary parts of plasmon excitation wavenumbers in
plasmonic crystal in terms of reciprocal lattice vector $G$ at orbital (a)
$E=0.8$ (b) $E=1.8$ (c) $E=2.5$ and (d) $E=3$. The possive branches are small
wavelength particle-like and negative large wavelength wave-like plasmon
excitation wavenumbers corresponding to each energy orbital.

Figure 6 shows the variations in imaginary wavenumber in terms of reciprocal
lattice vector for different energies. The symmetric nature of figure
indicates that excitations are dissipation free but energy exchange occurs
between particle-like and wave-like oscillations. The imaginary part of
wavevector at the band gaps play important role in Zener tunneling phenomenon
in semiconductor diodes [4]. Recently, the propagation of a single stream
electron beam studied in Ref. [68], reveals that the plasmonic excitations in
the intrinsic energy gap in plasmonic excitations, due to wave-particle branch
coupling below the critical value $E<2E_{p}$, leads to spacial growing/damping
of the wave-like/particle-like excitations. It is however, concluded that the
collective wave-particle interactions of dual-nature plasmonic excitations
inside energy band gaps is accompanied by enhancement of wave-like amplitude
and reduction in particle-like one through the space. This is a novel aspect
of collective quantum interaction phenomenon detailed investigation of which
is required in a future research. The positive/negative branches of imaginary
components in Fig. 6 are particle-like/wave-like damp/growth rates. These
imaginary wavenumbers always appear due to coupling of a wave-like excitation
dispersion branch with that of a particle-like which leads to the appearance
of energy gaps between separate bands.

Figure 7: (a) The electronic extended-zone energy band structure in free
electron model. (b) The plasmon excitation band in the limit $G\to\infty$. (c)
The extended-zone electronic band structure of plasmonic crystal with $G=3$.
(d) The extended-zone electronic band structure of plasmonic crystal with
$G=2$.

Figure 7 shows the calculation results for band structure in 10-coupled
approximation ($N=10$). Our result in Fig. 7(c) may be directly compared to
the empty lattice band structure of free electron model with $G=3$ shown in
Fig. 7(a). It is remarked that the band gap in our model occurs through
collective electronic excitations even in the absence of lattice potential.
The free electron lattice band structure may be obtained in our model by
setting the coupling electrostatic field $\phi$ to zero. Figure 7(b) shows
limiting case of $G\gg 1$ which is obviously the energy dispersion curve for
ordinary plasmon excitation in electron gas. Figures 7(c) and 7(d) show the
band structure for different values of $G$. It is remarked that band gaps take
place at the Brillouin zone boundary the size of which decreases with increase
in $k$ but increase with increase in $E$. The exact similarity between Figs.
7(c) and 4(b) shows that the problem of excitations in plasmonic crystals and
quantum multistream have the same root. As a further generalization one may
consider the effect of dynamic ions [69] on the plasmonic lattice or super-
lattice structure.

Figure 8: (a) Variation of the first two gap size at first Brilloun-zone
boundary ($k=G/2$) with the reciprocal lattice vector size, $G=2\pi/a$ with
$a$ being the lattice constant. (b) Variation of the first conduction band and
first two higher gap sizes at $k=G$ with the reciprocal lattice vector size,
$G=2\pi/a$. (c) Variation of the first two gap size at second Brilloun-zone
boundary ($k=3G/2$) with the reciprocal lattice vector size, $G=2\pi/a$. (d)
Variation of the first conduction band and two higher gap size at $k=2G$ with
the reciprocal lattice vector size, $G=2\pi/a$.

In Fig. 8 we have shown the variations in band gaps in various Brillouin zone
boundaries with respect to the reciprocal lattice vector. Figure 8(a) shows
the gap size variation at first zone boundary at $k=G/2$ for $\Delta
E_{12}=E_{2}-E_{1}$ and $\Delta E_{34}=E_{4}-E_{3}$. It is seen that $\Delta
E_{12}$ and $\Delta E_{34}$, respectively, maximize at $G\simeq 1.74903$ and
$G\simeq 1.44203$, that is, the gaps maximize at values of $G$ which tend to
decrease with increase of the energy bands. Comparing Figs. 5(d) and 8(a), it
is also noted that, the size of the band gap $\Delta E_{12}$ strongly depends
on the number of lattice sites in $N$-coupled approximation in the plasmonic
crystal. The band gap variations at $k=G$ are shown in Fig. 8(b). The first
conduction band, $\Delta E_{01}=E_{1}$ ($E=0$ is the Fermi level for fully
degenerate electron gas at zero temperature), at $k=G$ maximizes at the value
$G=1.37936$. It is seen that the gap sizes are larger for upper bands, while,
relatively lower in size with respect to those of the first zone. It is noted
that $\Delta E_{23}$ and $\Delta E_{45}$ at $k=3G/2$, respectively, maximize
at $G\simeq 1.1031$ and $G\simeq 0.9581$. Figures 8(c) and 8(d) reveal that
gaps at higher zone boundaries maximize at lower $G$-values. It is noted that
$\Delta E_{12}$ and $\Delta E_{34}$, respectively, maximize at $G\simeq
1.0061$ and $G\simeq 0.8282$. Also, $\Delta E_{01}$, $\Delta E_{12}$ and
$\Delta E_{34}$ at $k=2G$, respectively, maximize at $G\simeq 0.9709$,
$G\simeq 0.7738$ and $G\simeq 0.6537$. The variation of band structure in
terms of lattice parameter may play a fundamental role in effective plasmonic
crystal for technology. Note that the reciprocal lattice vector is related to
the equilibrium electron number-density through $n_{0}=1/a^{3}=G^{3}/8\pi^{3}$
where $G$ is normalized with respect to the plasmon wavevector
$k_{p}=\sqrt{2m_{e}\omega_{p}/\hbar}$ with $\omega_{p}=\sqrt{4\pi
e^{2}n_{0}/m_{e}}$ being the characteristic plasmon frequency. For instance,
Fig. 8(b) shows that the first conduction band gap maximizes for $G\simeq 1.4$
(in $k_{p}$ unit), i.e., for a equilibrium electron density of
$n_{0}=0.00368568n_{p}\simeq 1.28\times 10^{22}$cm-3 with
$n_{p}=\pi^{3}e^{6}m_{e}^{3}/16\hbar^{6}$ being the plasmon density, which is
slightly larger than the number density of a typical metallic elements [3] and
is very close to the critical screening point $E_{p}=2E_{F}$ [57]. On the
other hand, the maximum value of the valence-conduction gap is $\Delta
E_{01}\simeq 0.3$ (in $E_{p}$ unit). As an example the plasmon energy of
metallic Sodium is $E_{p}\simeq 5.9$eV giving the energy gap size of $\Delta
E_{01}\simeq 1.77$eV at $k=G$. The corresponding gap size at $k=2G$ is $\Delta
E_{01}\simeq 0.7$eV and becomes much smaller at higher electron momentum for
large $N$ at boundaries $k=NG$. Note that one has to take into account the
dynamic effects of lattice ions [69] which leads to sinking the conduction
band into the Fermi electron sea. The calculated amount of ion potential
effect on the energy band gap in the first-order perturbation approximation is
known to be constant [3] and independent of the number $N$ at boundaries
$k=NG$.

##  VII Electron-Lattice Binding Effect

In this section we would like to study the effect of electronic binding to
lattice sites on the energy band structure of plasmon excitations, in the
empty lattice approximation. To this end, we consider the following normalized
and linearized non-Hermitian system, which includes the spacial damping
effect. After the separation of spatiotemporal variables, one obtains

$\displaystyle\frac{{{\partial^{2}}{{\cal
M}_{G}}}}{{\partial{x^{2}}}}+2\xi\frac{{\partial{{\cal M}_{G}}}}{{\partial
x}}+\Phi{{\cal M}_{G}}+E{{\cal M}_{G}}=0,$ (21a)
$\displaystyle\frac{{{\partial^{2}}\Phi}}{{\partial{x^{2}}}}+2\xi\frac{{\partial{\Phi}}}{{\partial
x}}-\sum\limits_{G}{|{\cal M}_{G}|}=0,$ (21b)

where, $\xi$ denotes the strength of plasmon oscillation damping, due to the
electronic binding to the periodic lattice sites. However, in this simplified
model, we do not want to go into details of the dependence of the damping
parameter on other electron gas parameters, such as the equilibrium electron
number-density and temperature. It is evident that, the system (21) should
admit the general solution, ${\cal
M}_{G}(x)=\Psi_{N}(x)\exp(iNGx-\xi|x-Na|+i\Theta_{N})$, with the time
dependent solution as ${\cal N}_{G}(x,t)={\cal M}_{G}(x)\exp(-i\Omega t)$,
where the functions $\Psi_{G}(x)$ and $\Phi(x)$ satisfy the following
$N$-coupled system

$\displaystyle\frac{{{d^{2}}{\Psi_{1}}}}{{d{x^{2}}}}+2iG_{1}\frac{{d{\Psi_{1}}}}{{dx}}+\Phi+(E-G_{1}^{2}-\xi^{2}){\Psi_{1}}=0,$
(22a)
$\displaystyle\frac{{{d^{2}}{\Psi_{2}}}}{{d{x^{2}}}}+2iG_{2}\frac{{d{\Psi_{2}}}}{{dx}}+\Phi+(E-G_{2}^{2}-\xi^{2}){\Psi_{2}}=0,$
(22b) $\displaystyle\hskip 28.45274pt\vdots\hskip 56.9055pt\vdots\hskip
56.9055pt\vdots\hskip 56.9055pt$ (22c)
$\displaystyle\frac{{{d^{2}}{\Psi_{N}}}}{{d{x^{2}}}}+2i{G_{N}}\frac{{d{\Psi_{N}}}}{{dx}}+\Phi+(E-G_{N}^{2}-\xi^{2}){\Psi_{N}}=0,$
(22d)
$\displaystyle\frac{{{d^{2}}\Phi}}{{d{x^{2}}}}-\Psi_{1}-\Psi_{2}-\cdots-\Psi_{N}=0,$
(22e)

Fourier analysis of which results in the following eigenvalue system

$\left({\begin{array}[]{*{20}{c}}{E-{Q_{1}}}&0&\ldots&{\begin{array}[]{*{20}{c}}{\begin{array}[]{*{20}{c}}{}\hfil&{}\hfil&0&{\begin{array}[]{*{20}{c}}{}\hfil&{}\hfil\end{array}}\end{array}}&1\end{array}}\\\
0&{E-{Q_{2}}}&\ldots&{\begin{array}[]{*{20}{c}}{\begin{array}[]{*{20}{c}}{}\hfil&{}\hfil&0&{\begin{array}[]{*{20}{c}}{}\hfil&{}\hfil\end{array}}\end{array}}&1\end{array}}\\\
\vdots&\vdots&\ddots&{\begin{array}[]{*{20}{c}}{\begin{array}[]{*{20}{c}}{}\hfil&{}\hfil&{\begin{array}[]{*{20}{c}}\vdots&{}\hfil\end{array}}&{}\hfil\end{array}}&\vdots\end{array}}\\\
{\begin{array}[]{*{20}{c}}0\\\
1\end{array}}&{\begin{array}[]{*{20}{c}}\ldots\\\
1\end{array}}&{\begin{array}[]{*{20}{c}}0\\\
\cdots\end{array}}&{\begin{array}[]{*{20}{c}}{\begin{array}[]{*{20}{c}}{E-{Q_{N}}}&1\end{array}}\\\
{\begin{array}[]{*{20}{c}}{}\hfil&{\begin{array}[]{*{20}{c}}{}\hfil&1&{\begin{array}[]{*{20}{c}}{}\hfil&{\begin{array}[]{*{20}{c}}{}\hfil&{{k^{2}}}\end{array}}\end{array}}\end{array}}\end{array}}\end{array}}\end{array}}\right)\left({\begin{array}[]{*{20}{c}}{{\Psi_{11}}}\\\
{\begin{array}[]{*{20}{c}}{{\Psi_{21}}}\\\ \vdots\\\
{{\Psi_{N1}}}\end{array}}\\\
{{\Phi_{1}}}\end{array}}\right)=\left({\begin{array}[]{*{20}{c}}0\\\
{\begin{array}[]{*{20}{c}}0\\\ \vdots\\\ 0\end{array}}\\\
0\end{array}}\right),$ (23)

where $Q_{N}=(k+G_{N})^{2}-\xi^{2}$.

Figure 9: Variation in the electronic band structure of plasmon excitations in
a plasmonic crystal with $G=2\pi/a$ for (a) $\xi=0.1$, (b) $\xi=0.5$, (c)
$\xi=1$ and (d) $\xi=1.5$.

Figure 9 shows the effect of electron-lattice binding on the plasmonic band
structure of periodic system with $G=2$ and different values of the binding
strength parameter, $\xi$. It is clearly remarked that the increase in binding
strength leads to the overall shift of plasmonic excitation bands to higher
energies. It is further remarked that the energy gaps at the long wavelength
limit $k=0$ decrease sharply by increase of the binding parameter. This is due
to the significant effect of the electron-lattice binding on the wave-like
branch rather than the particle-like one. It is also clearly remarked that the
ground state electronic valence-conduction gap, $\Delta E_{01}$ at $k=NG$,
through which the valence electrons can tunnel in the nearly free electron
model, becomes smaller for larger values of the electron momentum, $\hbar k$.
For higher binding strength regime ($\xi>1$), shown in Fig. 9(d), which we
call the tight-binding limit where the electrons are tightly bound to the
lattice sites, the first conduction plasmon energy band shift to much higher
energies with the band inaccessible to Fermi electrons at $E=0$
($\epsilon=\mu_{0}$) at zero temperature limit, thus, leading to insulating
solid-state plasmon gap. Therefore, the critical value of binding parameter,
$\xi$, may provide a quantitative measure for the Mott metal-insulator
transition phenomenon in terms of the ground state gap energy, $\Delta
E_{01}$, at zero temperature limit, where electro-hole process can occur. At
finite temperature, on the other hand, electrons can excite to much higher
energy bands and collective phenomenon become more pronounced.

Figure 10: (a) Variation of the first two gap size at $k=0$ with the
reciprocal lattice vector size, $G=2$. (b) Variation of the first two gap size
at first Brilloun-zone boundary ($k=G/2$) with the reciprocal lattice vector
size, $G=2$. (c) Variation of the first two gap size at secong Brilloun-zone
boundary ($k=G$) with the reciprocal lattice vector size, $G=2$. (d) Variation
of the conduction band height, $\Delta E_{10}$, for $k=0,G$ with the
reciprocal lattice vector size, $G=2$.

Figure 10 shows the variations in various band gap sizes for 1D plasmonic
crystal with reciprocal lattice vector size $G=2$ in terms of the electron
binding parameter. Figure 10(a) shows the long wavelength ($k=0$) gap size
variation. It is remarked that, with increase in the binding parameter the
energy gaps $\Delta E_{45}$ and $\Delta E_{23}$ decrease sharply and saturate
to the same value for large $\xi$, which is also clearly evident from Fig. 9.
Figure 10(b) shows that the energy gaps $\Delta E_{34}$ and $\Delta E_{12}$
also decrease with increase in the value of $\xi$, but, with lower rate
compared to gaps at $k=0$. The variations in the gaps $k=G$ of Fig. 10(c)
becomes very small, indicating that the energy gaps at larger $k$ are less
affected by the binding parameter variations. Variation in the first plasmon
conduction bands $\Delta E_{01}$ for values of wavevectors $k=0,G$ are shown
in Fig. 10(d). It is seen that $\Delta E_{01}$ at $k=0$ and $k=G$ increase
with increase in $\xi$ and become identical for large $\xi$. It is concluded
that in the tight-binding limit, $\xi\gg 1$, the band gaps tend to close and
we obtain a free electron-like dispersion, similar to Fig. 7(a), with a very
large ground state gap $\Delta E_{01}\gg 1$.

##  VIII Plasmon-Phonon Coupling Effect

As a final remark, we would like to consider the effect of heavy species like
dynamic ions on the band structure of 1D plasmonic crystals. In such a case we
have a $N+1$-coupled system which may be written as

$\displaystyle\frac{{{d^{2}}{\Psi_{1}}}}{{d{x^{2}}}}+2iG_{1}\frac{{d{\Psi_{1}}}}{{dx}}+\Phi+(E-G_{1}^{2}){\Psi_{1}}=0,$
(24a)
$\displaystyle\frac{{{d^{2}}{\Psi_{2}}}}{{d{x^{2}}}}+2iG_{2}\frac{{d{\Psi_{2}}}}{{dx}}+\Phi+(E-G_{2}^{2}){\Psi_{2}}=0,$
(24b) $\displaystyle\hskip 28.45274pt\vdots\hskip 56.9055pt\vdots\hskip
56.9055pt\vdots\hskip 56.9055pt$ (24c)
$\displaystyle\frac{{{d^{2}}{\Psi_{N}}}}{{d{x^{2}}}}+2i{G_{N}}\frac{{d{\Psi_{N}}}}{{dx}}+\Phi+(E-G_{N}^{2}){\Psi_{N}}=0,$
(24d)
$\displaystyle\sigma\frac{{{d^{2}}{\Psi_{i}}}}{{d{x^{2}}}}-\Phi+(E+\mu){\Psi_{i}}=0,$
(24e)
$\displaystyle\frac{{{d^{2}}\Phi}}{{d{x^{2}}}}-\Psi_{1}-\Psi_{2}-\cdots-\Psi_{N}+\Psi_{i}=0,$
(24f)

where $\Psi_{i}$ denotes the ion wavefunction and $\sigma=m_{e}/m_{i}$ is the
electron to ion mass ratio. As before, the Fourier analysis of the
$N+1$-coupled system (24) leads to the desired energy bands.

Figure 11: Band structure of 1D plasmonic crystal with $G=2$ in the presence
of dynamic ions for (a) classical electron gas with $\mu=0$ and (b) degenerate
electron gas with $\mu=0.8$. The variation in the gap between (c) first band
and (d) second band from top of Fermi level with the change in normalized
chemical potential.

Figure 11 shows the band structure for given values of the chemical potential
and reciprocal lattice vector. The effect of heavy species such as dynamic
ions on the energy dispersion of plasmonic excitations has recently been
studied in [69]. For the electron gas in semiconductor regime with $\mu=0$, as
shown in Fig. 11(a), it is remarked that a nearly flat band (phonon-like low
energy band) appears at $E=0$ (the Fermi level). The flatness of a conduction
band is an indication of decreased mobility of electrons which is caused by
electrostatic coupling of electrons to inertial ions. However, existence of
such band is important for long wavelength phonon-assisted plasmon excitations
in metals in the zero temperature limit. For a fully degenerate electron gas
with $\mu=0.8$, the lowest energy band resides well below the Fermi energy
level embracing a large amount of nearly free electrons in the gas, as seen in
Fig. 11(b). The variation in gap between first band and the top of Fermi level
$\Delta E_{01}$ is shown in Fig. 11(c) in terms of the normalized chemical
potential. It is seen that the first band at $k=G$ touches the Fermi level at
exact value of $\mu=0$ above/under which value the first band is belove/above
the Fermi energy level. Moreover, Fig. 11(d) depicts the gap of second band
from Fermi level at $k=G$ for $G=2$. The second band never touches the Fermi
level for the chemical potential values in semiconductors and metallic density
regime. However, the second energy band gap from $E=0$ decreases as the
chemical potential increases. It should be noted that plasmonic crystals,
unlike ordinary solids, can constitute from different charged species other.
Also, the lattice sites can be interfaces between different plasmonic
(metallic and semiconductor) layers, so called superlattice configuration,
which do not contribute electrostatic potential to energy band structure.
Therefore, current empty lattice model may well apply to a wide range of
plasmonic crystal and superlattice configurations.

We have already considered the multistream electrostatic systems which only
include electrons and ions. However, the simplified current model may be
further generalized to include multispecies complex plasmas with a wide range
of mass and charge-state spectrum or even gravitationally coupled uncharged
quantum fluids. As discussed earlier, the energy band gap structure formation
in quantum multistream system is expected be the origin of collisionless
quantum stream instability and Landau damping effects. Therefore it is
concluded that, these effects not only are characteristics of electrostatic
systems, but also are inevitable in uncharged mass/spin multistream systems
coupled through gravitational/magnetic potentials. The fundamental difference
between electrostatic and gravitational Landau damping is that for
gravitational case the damping occurs for wavenumbers larger than the Jeans
wavenumber below which Jeans instability occurs.

##  IX Conclusion

We used the multifluid model to study the plasmonic excitations in electron
gas with arbitrary degree of degeneracy by reducing the quantum hydrodynamic
model into the $N$-coupled pseudoforce system. The energy band structure of a
multistream system was obtained by linearizing the coupled differential
equations which indicated that the energy bands form due to discrete stream
filaments in the system and mode coupling by collective electrostatic
interactions. Such velocity filaments may also be the root to collisionless
damping and stream instability by gap opening very similar to the crystalline
solids. Current model, generalized to virtual streams, was used to calculate
the electronic band structure in one-dimensional plasmonic crystal. The
dependence of energy band gaps on the lattice spacing is also studied in
detail. The electronic band structure of a electron system can have essential
effect on many characteristics of collective excitations in plasmonic crystals
and metallic superlattices. We further studied the effect electron-lattie
binding effect on the energy band structure of plasmonic crystals which
indicates that with increase in the strength of the electron binding the first
energy conduction band shifts to higher values where unaccessible to electrons
at the Fermi sea. Inclusion of dynamic inertial ions in the plasmonic
crystals, on the other hand, reveals that for degenerate electrons a flat-like
ground state energy band appears inside the Fermi sea of free electrons due to
electrostatic interaction between free electrons and ions decreasing the
electron mobility substantially at this level. Therefore, current model of
plasmonic excitations is capable of incorporating a wide range realistic
features of electron dynamics in one-dimensional periodic structures and show
akin similarities in band structure between multistream electron gas and
plasmonic crystals.

##  X Data Availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * [1] F. F. Chen, Introduction to Plasma Physics and Controlled Fusion, 2nd ed. (Plenum Press, New York, London, 1984).
  * [2] N. A. Krall and A. W. Trivelpeice, ”Principles of Plasma Physics”, (San francisco Press, San francisco 1986).
  * [3] C. Kittel, Introduction to Solid State Physics, (John Wiely and Sons, New York, 1996), 7th ed.
  * [4] N. W. Ashcroft and N. D. Mermin, Solid State Physics (Saunders College Publishing, Orlando, 1976).
  * [5] S. Ummethala, T. Harter, K. Koehnle, et al., Nat. Photonics 13, 519(2019); doi.org/10.1038/s41566-019-0475-6
  * [6] P. A. Markovich, C.A. Ringhofer, and C. Schmeister, Semiconductor Equations (Springer, Berlin, 1990).
  * [7] C. Gardner, SIAM, J. Appl. Math. 54 409(1994).
  * [8] G. Manfredi, Phys. Plasmas 25, 031701(2018); https://doi.org/10.1063/1.5026653
  * [9] S. A. Maier, Plasmonics: Fundamentals and Aplications, Springer Science Business Media LLC (2007).
  * [10] H. Haug and S. W. Koch, ”Quantum theory of the optical and electronic properties of semiconductors”, World Scientific, 2004,
  * [11] C. Hu, Modern Semiconductor Devices for Integrated Circuits (Prentice Hall, Upper Saddle River, New Jersey, 2010) 1st ed.
  * [12] K. Seeger, Semiconductor Physics (Springer, Berlin, 2004) 9th ed.
  * [13] H. A. Atwater, Sci. Am. 296, 56(2007).
  * [14] C. Calvero, Nat. Photonic 8, 95(2014).
  * [15] Jacob B. Khurgin, Nanophotonics 9(2), 453(2020).
  * [16] H. A. Atwater and A. Polman, Nat. Mater. 9, 205(2010).
  * [17] Y. Tian and T. Tatsuma, J. Am. Chem. Soc. 127, 7632(2005).
  * [18] S. Ichimaru, Rev. Mod. Phys. 54, 1017 (1982).
  * [19] S. Ichimaru, H. Iyetomi, and S. Tanaka, Phys. Rep. 149, 91 (1987).
  * [20] S. Ichimaru, Statistical Physics: Condensed Plasmas (Addison Wesely, New York, 1994).
  * [21] A. D. Yofee, Adv. Phys., 42, 173-262(1993), DOI: 10.1080/00018739300101484
  * [22] Jianing Chen, Michela Badioli, Pablo Alonso-González, Sukosin Thongrattanasiri, Florian Huth, Johann Osmond, Marko Spasenović, Alba Centeno, Amaia Pesquera, Philippe Godignon, Amaia Zurutuza Elorza, Nicolas Camara, F. Javier García de Abajo, Rainer Hillenbrand, Frank H. L. Koppens, Nat. 487 77(2012).
  * [23] Z. Fei, A. S. Rodin, G. O. Andreev, W. Bao, A. S. McLeod, M. Wagner, L. M. Zhang, Z. Zhao, M. Thiemens, G. Dominguez, M. M. Fogler, A. H. Castro Neto, C. N. Lau, F. Keilmann, D. N. Basov, Nat. 487 82(2012); doi:10.1038/nature11253
  * [24] Hugen Yan, Tony Low, Wenjuan Zhu, Yanqing Wu, Marcus Freitag, Xuesong Li, Francisco Guinea, Phaedon Avouris, Fengnian Xia, Nat. Photonics 7 394(2013); doi:10.1038/nphoton.2013.57
  * [25] D. Bohm and D. Pines, Phys. Rev. 92 609(1953).
  * [26] Bohm, D. Phys. Rev. 85, 166–179 (1952).
  * [27] Bohm, D. Phys. Rev. 85, 180–193 (1952).
  * [28] D. Pines, Phys. Rev. 92 609(1953).
  * [29] P. Levine and O. V. Roos, Phys. Rev, 125 207(1962).
  * [30] N. H. March, Phys. Rev. 92, 510(1953).
  * [31] Katuro Sawada, Phys. Rev. 106, 372(1957).
  * [32] W. Kohn and L. J. Sham, Phys. Rev. 137, A1697(1965).
  * [33] A. L. Fetter and J. D. Walecka, Quantum Theory of Many-Particle Systems,. McGraw-Hill 1971.
  * [34] G. D. Mahan, Many-particle physics, 2nd edition, chapter 5 (Plenum press, New York, 1990).
  * [35] D. Pines and P. Nozieres, The Theory of Quantum Liquids (Addison-Wesley, 1968).
  * [36] E. Fermi and E. Teller, Phys. Rev. 72, 399 (1947).
  * [37] S. Chandrasekhar, ”An Introduction to the Study of Stellar Structure”, Chicago, Ill. (The University of Chicago press), (1939), p.392.
  * [38] F. Hoyle and W. A. Fowler, Astrophys. J. 132, 565(1960).
  * [39] Y. Klimontovich and V. P. Silin, in Plasma Physics, edited by J. E. Drummond (McGraw-Hill, New York, 1961).
  * [40] F. Haas, Quantum Plasmas: An Hydrodynamic Approach (Springer, New York, 2011).
  * [41] G. Manfredi, “How to model quantum plasmas,” Fields Inst. Commun. 46, 263–287 (2005); in Proceedings of the Workshop on Kinetic Theory (The Fields Institute, Toronto, Canada 2004): http://arxiv.org/abs/quant–ph/0505004.
  * [42] F. Haas, Quantum Plasmas: An Hydrodynamic Approach (Springer, New York, 2011).
  * [43] M. Bonitz, Zh. A. Moldabekov, and T. S. Ramazanov, Phys. Plasmas 26, 090601 (2019); doi.org/10.1063/1.5097885
  * [44] P. K. Shukla and B. Eliasson, Phys. Rev. Lett. 99, 096401(2007).
  * [45] L Stenflo Phys. Scr. T50 15(1994).
  * [46] P. K. Shukla, B. Eliasson, and L. Stenflo Phys. Rev. E 86, 016403(2012).
  * [47] G. Brodin and M. Marklund, New J. Phys. 9, 277(2007).
  * [48] M. Marklund and G. Brodin, Phys. Rev. Lett. 98, 025001(2007).
  * [49] N. Crouseilles, P. A. Hervieux, and G. Manfredi, Phys. Rev. B 78, 155412 (2008).
  * [50] F. Haas, G. Manfredi, P. K. Shukla, and P.-A. Hervieux, Phys. Rev. B, 80, 073301 (2009).
  * [51] B. Eliasson and P. K. Shukla, Phys. Scr. 78, 025503 (2008).
  * [52] L. Stenflo, Phys. Scripta 14, 320(1967).
  * [53] L. Stanton and M. S. Murillo, Phys. Rev. E 91, 033104(2015).
  * [54] D. Michta, F. Graziani, and M. Bonitz, Contrib. Plasma Phys. 55, 437 (2015).
  * [55] S. Ali, H. Terças, J. T. Mendonca, Phys. Rev. B, 83(15), 153401(2011); doi:10.1103/PhysRevB.83.153401
  * [56] S. Ali, Shafiq Ur Rehman, and Z. J. Ding, Phys. Plasmas 25, 082115 (2018); doi: 10.1063/1.5041297
  * [57] M. Akbari-Moghanjoughi, Alireza Abdikian, and Arash Phirouznia1, Phys. Plasmas 27, 042107 (2020); doi.org/10.1063/5.0004857
  * [58] G. Manfredi, 86(2), 825860201(2020); doi:10.1017/S0022377820000240
  * [59] G. Manfredi and F. Haas, Phys. Rev. B 64, 075316 (2001);
  * [60] J. Hurst, K. L. Simon, P. A. Hervieux, G. Manfredi and F. Haas, Phys. Rev. B 93, 205402(2016).
  * [61] E. Madelung, Quantentheorie in hydrodynamischer Form. Zeitschr. f. Phys. 40, 322–326 (1926).
  * [62] M. Akbari-Moghanjoughi, Phys. Plasmas 22, 022103 (2015); ibid. 22, 039904 (E) (2015).
  * [63] M. Akbari-Moghanjoughi, Phys. Plasmas, 26, 012104 (2019); doi.org/10.1063/1.5078740
  * [64] J. Dawson, Phys. Fluids 4, 869 (1961).
  * [65] F. Haas, G. Manfredi, and M. Feix, Phys. Rev. E 62, 2763 (2000); doi.org/10.1103/PhysRevE.62.2763
  * [66] M. Akbari-Moghanjoughi, Phys. Plasmas, 26, 022111 (2019); doi.org/10.1063/1.5083150
  * [67] M. Akbari-Moghanjoughi, Phys. Plasmas, 28, 022109 (2021); doi.org/10.1063/5.0039067
  * [68] M. Akbari-Moghanjoughi, Phys. Plasmas, 28, 082109 (2021); doi: 10.1063/5.0057662
  * [69] M. Akbari-Moghanjoughi, Phys. Plasmas, 27, 112107 (2020); doi.org/10.1063/5.0026298


11institutetext:  Microsoft Research, Redmond, Washington<EMAIL_ADDRESS>22institutetext: Computer Science Department, University of New Mexico,
Albuquerque and the Santa Fe Institute, Santa Fe, New Mexico<EMAIL_ADDRESS>
# On the 2-colorability of random hypergraphs

Dimitris Achlioptas and Cristopher Moore Supported by NSF grant PHY-0071139,
the Sandia University Research Program, and Los Alamos National Laboratory.

###### Abstract

A 2-coloring of a hypergraph is a mapping from its vertices to a set of two
colors such that no edge is monochromatic. Let $H_{k}(n,m)$ be a random
$k$-uniform hypergraph on $n$ vertices formed by picking $m$ edges uniformly,
independently and with replacement. It is easy to show that if $r\geq
r_{c}=2^{k-1}\ln 2-(\ln 2)/2$, then with high probability $H_{k}(n,m=rn)$ is
not 2-colorable. We complement this observation by proving that if $r\leq
r_{c}-1$ then with high probability $H_{k}(n,m=rn)$ is 2-colorable.

##  1 Introduction

For an integer $k\geq 2$, a $k$-uniform hypergraph $H$ is an ordered pair
$H=(V,E)$, where $V$ is a finite non-empty set, called the set of vertices of
$H$, and $E$ is a family of distinct $k$-subsets of $V$, called the edges of
$H$. For general hypergraph terminology and background see [5]. A $2$-coloring
of a hypergraph $H=(V,E)$ is a partition of its vertex set $V$ into two
(color) classes so that no edge in $E$ is monochromatic. A hypergraph is
2-colorable if it admits a 2-coloring.

The property of 2-colorability was introduced and studied by Bernstein [6] in
the early 1900s for infinite hypergraphs. The 2-colorability of finite
hypergraphs, also known as “Property B” (a term coined by Erdős in reference
to Bernstein), has been studied for about eighty years (e.g. [4, 10, 11, 15,
16, 19, 20]). For $k=2$, i.e. for graphs, the problem is well understood since
a graph is 2-colorable if and only if it has no odd cycle. For $k\geq 3$,
though, much less is known and deciding the 2-colorability of $k$-uniform
hypergraphs is NP-complete [17].

In this paper we discuss the 2-colorability of random $k$-uniform hypergraphs
for $k\geq 3$. (For the evolution of odd cycles in random graphs see [12].)
Let $H_{k}(n,m)$ be a random $k$-uniform hypergraph on $n$ vertices, where the
edge set is formed by selecting uniformly, independently and with replacement
$m$ out of all possible ${n\choose k}$ edges. We will study asymptotic
properties of $H_{k}(n,m)$ when $k\geq 3$ is arbitrary but fixed while $n$
tends to infinity. We will say that a hypergraph property $A$ holds with high
probability (w.h.p.) in $H_{k}(n,m)$ if
$\lim_{n\rightarrow\infty}\Pr[H_{k}(n,m)\mbox{ has }A]=1$. The main question
in this setting is:

As $m$ is increased, when does $H_{k}(n,m)$ stop being 2-colorable?

It is popular to conjecture that the transition from 2-colorability to
non-2-colorability is sharp. That is, it is believed that for each $k\geq 3$,
there exists a constant $r_{k}$ such that if $r<r_{k}$ then $H_{k}(n,m=rn)$ is
w.h.p. 2-colorable, but if $r>r_{k}$ then w.h.p. $H_{k}(n,m=rn)$ is not
2-colorable. Determining $r_{k}$ is a challenging open problem, closely
related to the satisfiability threshold conjecture for random $k$-SAT.
Although $r_{k}$ has not been proven to exist, we will take the liberty of
writing $r_{k}\geq r^{*}$ to denote that for $r<r^{*}$, $H_{k}(n,rn)$ is
2-colorable w.h.p. (and analogously for $r_{k}\leq r^{*}$).

A relatively recent result of Friedgut [13] supports this conjecture as it
gives a non-uniform sharp threshold for hypergraph 2-colorability. Namely, for
each $k\geq 3$ there exists a sequence $r_{k}(n)$ such that if
$r<r_{k}(n)-\epsilon$ then w.h.p. $H_{k}(n,rn)$ is 2-colorable, but if
$r>r_{k}(n)+\epsilon$ then w.h.p. $H_{k}(n,rn)$ is not 2-colorable. We will
find useful the following immediate corollary of this sharp threshold.

###### Corollary 1

If

${\liminf_{n\rightarrow\infty}}\Pr[H_{k}(n,r^{*}n)\mbox{ is
2-colorable}]>0\enspace,$

then for $r<r^{*}$, $H_{k}(n,rn)$ is 2-colorable w.h.p.

Alon and Spencer [3] were the first to give bounds on the potential value of
$r_{k}$. In particular, they observed that the expected number of 2-colorings
of $H_{k}(n,m=rn)$ is $o(1)$ if $2(1-2^{1-k})^{r}<1$, implying

$r_{k}<2^{k-1}\ln 2-\frac{\ln 2}{2}\enspace.$ (1)

Their main contribution, though, was providing a lower bound on $r_{k}$.
Specifically, by applying the Lovász Local Lemma, they were able to show that
if $r=c\,2^{k}/k^{2}$ then w.h.p. $H_{k}(n,rn)$ is 2-colorable, for some small
constant $c>0$.

In [1], Achlioptas, Kim, Krivelevich and Tetali reduced the asymptotic gap
between the upper and lower bounds of [3] from order $k^{2}$ to order $k$. In
particular, they proved that there exists a constant $c>0$ such that if $r\leq
c\,2^{k}/k$ then a simple, linear-time algorithm w.h.p. finds a 2-coloring of
$H_{k}(n,rn)$. Their algorithm was motivated by algorithms for random $k$-SAT
due to Chao and Franco [7] and Chvátal and Reed [8]. In fact, those algorithms
give a similar $\mathrm{\Omega}(2^{k}/k)$ lower bound on the random $k$-SAT
threshold which, like $r_{k}$, can also be easily bounded as
$\mathrm{O}(2^{k})$.

Very recently, the authors eliminated the gap for the random $k$-SAT
threshold, determining its value within a factor of two [2]. The proof amounts
to applying the “second moment” method to the set of satisfying truth
assignments whose complement is also satisfying. Alternatively, one can think
of this as applying the second moment method to the number of truth
assignments under which every $k$-clause contains at least one satisfied
literal and at least one unsatisfied literal, i.e. which satisfy the formula
when interpreted as a random instance of Not-All-Equal $k$-SAT (NAE $k$-SAT).

Here we extend the techniques of [2] and apply them to hypergraph
2-colorability. This allows us to determine $r_{k}$ within a small additive
constant.

###### Theorem 1.1

For every $\epsilon>0$ and all $k\geq k_{0}(\epsilon)$,

$r_{k}\geq 2^{k-1}\ln 2-\frac{\ln 2}{2}-\frac{1+\epsilon}{2}\enspace.$

Our method actually yields an explicit lower bound for $r_{k}$ for each value
of $k$ as the solution to a simple equation (yet one without a pretty closed
form, hence Theorem 1.1). Below we compare this lower bound to the upper bound
of (1) for small values of $k$. The gap converges to 1/2 rather rapidly.

$\begin{array}[]{c|ccccccccccc}k&3&4&5&7&9&11&12\\\ \hline\cr\mbox{Lower
bound}&\;\;3/2&\;49/12&\;9.973&\;43.432&\;176.570&\;708.925&\;1418.712\\\
\mbox{Upper
bound}&\;\;2.409&5.191&10.740&44.014&177.099&709.436&1419.219\end{array}$
Table 1: Upper and lower bounds for $r_{k}$

##  2 Second moment and NAE $k$-SAT

We prove Theorem 1.1 by applying the following version of the second moment
method (see Exercise 3.6 in [18]) to the number of 2-colorings of
$H_{k}(n,m=rn)$.

###### Lemma 1

For any non-negative integer-valued random variable $X$,

$\Pr[X>0]\,\geq\,\frac{{\bf E}[X]^{2}}{{\bf E}[X^{2}]}\enspace.$ (2)

In particular, if $X$ is the number of 2-colorings of $H_{k}(n,m=rn)$, we will
prove that for all $\epsilon>0$ and all $k\geq k_{0}(\epsilon)$, if
$r=2^{k-1}\ln 2-\ln 2/2-(1+\epsilon)/2$ then there exists some constant
$C=C(k)$ such that

${\bf E}[X^{2}]<C\times{\bf E}[X]^{2}\enspace.$

By Lemma 1, this implies $\Pr[X>0]=\Pr[H_{k}(n,rn)\mbox{ is
2-colorable}]>1/C$. Theorem 1.1 follows by invoking Corollary 1.

This approach parallels the one taken recently by the authors for random NAE
$k$-SAT [2]. Naturally, what differs is the second-moment calculation which
here is prima facie significantly more involved.

We start our exposition by outlining the NAE $k$-SAT calculation of [2]. This
serves as a warm up for our calculations and allows us to state a couple of
useful lemmata from [2]. We then proceed to outline the proof of our main
result, showing the parallels with NAE $k$-SAT and reducing the proof of
Theorem 1.1 to the proof of three independent lemmata.

The first such lemma is specific to hypergraph 2-colorability and expresses
${\bf E}[X^{2}]$ as a multinomial sum. The second one is a general lemma about
bounding multinomial sums by a function of their largest term and is perhaps
of independent interest. It generalizes Lemma 2 of [2], which we state below.
After applying these two lemmata, we are left to maximize a three-variable
function parameterized by $k$ and $r$. This is analogous to NAE $k$-SAT,
except that there we only have to deal with a one-variable function, similarly
parameterized. That simpler maximization, in fact, amounted to the bulk of the
technical work in [2]. Luckily, here we will be able to get away with much
less work: a convexity argument will allow us to reduce our three-dimensional
optimization precisely to the optimization in [2].

###  2.1 Proof outline for NAE $k$-SAT

Let $Y$ be the number of satisfying assignments of a random NAE $k$-SAT
formula with $n$ variables and $m=rn$ clauses. It is easy to see that ${\bf
E}[Y]=2^{n}(1-2^{1-k})^{rn}$. Then ${\bf E}[Y^{2}]$ is the sum, over all
ordered pairs of truth assignments, of the probability that both assignments
in the pair are satisfying. It is not hard to show that if two assignments
assign the same value to $z=\alpha n$ variables, then the probability that
both are satisfying is

$p(\alpha)=1\,-\,2^{1-k}\left(2-\alpha^{k}-(1-\alpha)^{k}\right)\enspace.$

Since there are $2^{n}{n\choose z}$ pairs of assignments sharing $z$
variables, we have

$\frac{{\bf E}[Y^{2}]}{{\bf E}[Y]^{2}}=\sum_{z=0}^{n}{n\choose
z}\left[\frac{1}{2}\,\left(\frac{p(z/n)}{(1-2^{1-k})^{2}}\right)^{r}\right]^{n}\enspace.$

To bound such sums within a constant factor, we proved the following in [2].

###### Lemma 2

Let $f$ be any real positive analytic function and let

$S=\sum_{z=0}^{n}{n\choose z}\,f(z/n)^{n}\enspace.$

Define $0^{0}\equiv 1$ and let $g$ on $[0,1]$ be

$g(\alpha)=\frac{f(\alpha)}{\alpha^{\alpha}\,(1-\alpha)^{1-\alpha}}\enspace.$

If there exists $\alpha_{\max}\in(0,1)$ such that $g(\alpha_{\max})\equiv
g_{\max}>g(\alpha)$ for all $\alpha\neq\alpha_{\max}$, and
$g^{\prime\prime}(\alpha_{\max})<0$, then there exist constants $B$ and $C$
such that for all sufficiently large $n$

$B\times g_{\max}^{n}\,\leq\,S\leq\,C\times g_{\max}^{n}\enspace.$

Thus, using Lemma 2, bounding ${\bf E}[X^{2}]/{\bf E}[X]^{2}$ reduces to
maximizing

$\ell_{r}(\alpha)=\frac{1}{2\,\alpha^{\alpha}\,(1-\alpha)^{1-\alpha}}\,\left(\frac{p(\alpha)}{(1-2^{1-k})^{2}}\right)^{r}\enspace.$
(3)

Note now that $\ell_{r}(1/2)=1$ for all $r$ and that our goal is to find $r$
such that $\ell_{r}(\alpha)\leq 1$ for all $\alpha\in[0,1]$. Indeed, in [2] we
showed that

###### Lemma 3

[2] For every $\epsilon>0$, and all $k\geq k_{0}(\epsilon)$, if

$r\leq 2^{k-1}\ln 2-\frac{\ln 2}{2}-\frac{1+\epsilon}{2}$

then $\ell_{r}(1/2)=1>\ell_{r}(\alpha)$ for all $\alpha\neq 1/2$ and
$\ell_{r}^{\prime\prime}(1/2)<0$.

Thus, for all $r,k,\epsilon$ as in Lemma 3, we see that Lemma 2 implies ${\bf
E}[Y^{2}]/{\bf E}[Y]^{2}<C\times\ell_{r}(1/2)^{n}=C$, concluding the proof.

##  3 Proof outline for hypergraph 2-colorability

Let $X$ be the number of 2-colorings of $H_{k}(n,rn)$. Let $q=1-2^{1-k}$ and

$p(\alpha,\beta,\gamma)=1-\alpha^{k}-(1-\alpha)^{k}-\beta^{k}-(1-\beta)^{k}+\gamma^{k}+(\alpha-\gamma)^{k}+(\beta-\gamma)^{k}+(1-\alpha-\beta+\gamma)^{k}.$

We will prove that

###### Lemma 4

There exists a constant $A$ such that

$\frac{{\bf E}[X^{2}]}{{\bf
E}[X]^{2}}\leq\frac{1}{A^{2}}\sum_{z_{1}+\cdots+z_{4}=n}{n\choose
z_{1},z_{2},z_{3},z_{4}}\,\left(\frac{1}{4}\left(\frac{p\\!\left(\frac{z_{1}+z_{2}}{n},\frac{z_{1}+z_{3}}{n},\frac{z_{1}}{n}\right)}{q^{2}}\right)^{r}\right)^{n}\enspace.$

Similarly to NAE $k$-SAT we would like to bound this sum by a function of its
maximum term. To do this we will establish a multidimensional generalization
of the upper bound of Lemma 2.

###### Lemma 5

Let $f$ be any real positive analytic function and let

$S\;\;=\sum_{z_{1}+\cdots+z_{d}=n}{{n\choose
z_{1},\cdots,z_{d}}\,f(z_{1}/n,\cdots,z_{d-1}/n)^{n}}\enspace.$

Let $\displaystyle{Z=\left\\{(\zeta_{1},\ldots,\zeta_{d-1}):\zeta_{i}\geq
0\mbox{ for all $i$, and }\sum\zeta_{i}\leq 1\right\\}}$. Define $g$ on $Z$ as

$g(\zeta_{1},\ldots,\zeta_{d-1})=\frac{f(\zeta_{1},\ldots,\zeta_{d-1})}{\zeta_{1}^{\zeta_{1}}\cdots\zeta_{d-1}^{\zeta_{d-1}}\,(1-\zeta_{1}-\cdots-\zeta_{d-1})^{1-\zeta_{1}-\cdots-\zeta_{d-1}}}\enspace.$

If i) there exists $\vec{\zeta}_{\max}$ in the interior of $Z$ such that for
all $\vec{\zeta}\in Z$ with $\vec{\zeta}\neq\vec{\zeta}_{\max}$, we have
$g(\vec{\zeta}_{\max})\equiv g_{\max}>g(\vec{\zeta})$, and ii) the determinant
of the $(d-1)\times(d-1)$ matrix of second derivatives of $g$ is nonzero at
$\vec{\zeta}_{\max}$, then there exists a constant $D$ such that for all
sufficiently large $n$

$S\,<\,D\times g_{\max}^{n}\enspace.$

Applying Lemma 5 to the sum in Lemma 4 we see that we need to maximize

$g_{r}(\alpha,\beta,\gamma)=\frac{\displaystyle{\left(\frac{p(\alpha,\beta,\gamma)}{q^{2}}\right)^{r}}}{4\,\gamma^{\gamma}\,(\alpha-\gamma)^{\alpha-\gamma}\,(\beta-\gamma)^{\beta-\gamma}\,(1-\alpha-\beta+\gamma)^{1-\alpha-\beta+\gamma}}\enspace,$
(4)

where for convenience we defined $g_{r}$ in terms of $\alpha,\beta,\gamma$
instead of $\zeta_{1},\zeta_{2},\zeta_{3}$. We will show that $g_{r}$ has a
unique maximum at

$\vec{\zeta}^{*}=(1/2,1/2,1/4)\enspace.$

###### Lemma 6

For every $\epsilon>0$, and all $k\geq k_{0}(\epsilon)$ if

$r\leq 2^{k-1}\ln 2-\frac{\ln 2}{2}-\frac{1+\epsilon}{2}$

then $g_{r}(\vec{\zeta}^{*})=1>g_{r}(\vec{\zeta})$ for all $\vec{\zeta}\in Z$
with $\vec{\zeta}\neq\vec{\zeta}^{*}$. Moreover, the determinant of the matrix
of second derivatives of $g_{r}$ at $\vec{\zeta}^{*}$ is nonzero.

Therefore, for all $r,k,\epsilon$ as in Lemma 6

$\frac{{\bf E}[X^{2}]}{{\bf E}[X]^{2}}\;<\;\frac{D}{A^{2}}\times
g_{r}(\vec{\zeta}^{*})^{n}\;=\;D/A^{2}\enspace,$

completing the proof of Theorem 1.1 modulo Lemmata 4, 5 and 6.

The proof of Lemma 4 is a straightforward probabilistic calculation. The proof
of Lemma 5 is somewhat technical but follows standard asymptotic methods. To
prove Lemma 6 we will rely very heavily on Lemma 3. In particular, we will
show that all local maxima of $g_{r}$ occur within a one-dimensional subspace,
in which $g_{r}$ coincides with the function $\ell_{r}$ of (3). Specifically,
we prove

###### Lemma 7

If $(\alpha,\beta,\gamma)$ is a local extremum of $g_{r}$, then
$\alpha=\beta=1/2$.

This reduces our problem to the one-dimensional maximization for NAE $k$-SAT,
allowing us to easily prove Lemma 6.

Proof of Lemma 6. Observe that

$\displaystyle g_{r}(1/2,1/2,\gamma)=\ell_{r}(2\gamma)\enspace,$

where $\ell_{r}$ is the function defined in (3) for NAE $k$-SAT. Thus, the
inequality $g_{r}(\vec{\zeta}^{*})>g_{r}(\vec{\zeta})$ for
$\vec{\zeta}\neq\vec{\zeta}^{*}$ follows readily from Lemma 3, giving the
first part of the lemma.

To prove the condition on the determinant of the $3\times 3$ matrix of second
derivatives, a little arithmetic shows that at $\vec{\zeta}^{*}$ it is equal
to

$\frac{256\,(2^{k}-2-2kr+2k^{2}r)^{2}}{2^{4k}q^{4}}\,(4k(k-1)\,r-2^{2k}\,q^{2})\enspace.$

Thus, the determinant is negative whenever

$4k(k-1)\,r<2^{2k}\,q^{2}\enspace.$

For $k=3,4$ this is true for $r<3/2$ and $r<49/12$ respectively, while for
$k\geq 5$ it is true for all $r<\ln 2\times 2^{k-1}$. ∎

##  4 Proof of Lemma 4

Recall that $X$ denotes the number of 2-colorings of $H_{k}(n,m=rn)$.

###  4.1 First moment

Recall that

$q=1-2^{1-k}\enspace.$

The probability that a 2-coloring with $z=\alpha n$ black vertices and
$n-z=(1-\alpha)n$ white vertices makes a random hyperedge of size $k$
bichromatic is

$s(\alpha)=1-\alpha^{k}-(1-\alpha)^{k}\leq q\enspace.$

Summing over the $2^{n}$ colorings gives

${\bf E}[X]\,=\,\sum_{z=0}^{n}{n\choose z}\,s(z/n)^{rn}\enspace.$

To bound this sum from below we apply the lower bound of Lemma 2 with
$f(\alpha)=s(\alpha)^{r}$. In particular, it is easy to see that for all $r>0$

$g(\alpha)=\frac{s(\alpha)^{r}}{\alpha^{\alpha}(1-\alpha)^{1-\alpha}}=\frac{(1-\alpha^{k}-(1-\alpha)^{k})^{r}}{\alpha^{\alpha}(1-\alpha)^{1-\alpha}}$

is maximized at $\alpha=1/2$ and that $g(1/2)=2q^{r}$. Moreover, for any $k>1$

$g^{\prime\prime}(1/2)=-8\,(1-2^{1-k})^{r-1}\,\bigl{(}1+2^{1-k}(k(k-1)r-1)\bigr{)}\,<\,0\enspace.$

Therefore, we see that there exists a constant $A$ such that

${\bf E}[X]\,\geq\,A\times\left(2q^{r}\right)^{n}\enspace.$ (5)

###  4.2 Second moment

We first observe that ${\bf E}[X^{2}]$ equals the expected number of ordered
pairs $S,T$ of 2-partitions of the vertices such that both $S$ and $T$ are
2-colorings. Suppose that $S$ and $T$ have $\alpha n$ and $\beta n$ black
vertices respectively, while $\gamma n$ vertices are black in both. By
inclusion-exclusion a random hyperedge of size $k$ is bichromatic under both
$S$ and $T$ with probability $p(\alpha,\beta,\gamma)$, i.e.

$1-\alpha^{k}-(1-\alpha)^{k}-\beta^{k}-(1-\beta)^{k}+\gamma^{k}+(\alpha-\gamma)^{k}+(\beta-\gamma)^{k}+(1-\alpha-\beta+\gamma)^{k}\enspace.$

The negative terms above represent the probability that the hyperedge is
monochromatic under either $S$ or $T$, while the positive terms represent the
probability that it is monochromatic under both (potentially with different
colors). Since the $m=rn$ hyperedges are chosen independently and with
replacement, the probability that all $m=rn$ hyperedges are bichromatic is
$p(\alpha,\beta,\gamma)^{rn}$.

If $z_{1}$, $z_{2}$, $z_{3}$ and $z_{4}$ vertices are respectively black in
both assignments, black in $S$ and white in $T$, white in $S$ and black in
$T$, and white in both, then $\alpha=(z_{1}+z_{2})/n$, $\beta=(z_{1}+z_{3})/n$
and $\gamma=z_{1}/n$. Thus,

${\bf E}[X^{2}]=\sum{n\choose
z_{1},z_{2},z_{3},z_{4}}\,p\\!\left(\frac{z_{1}+z_{2}}{n},\,\frac{z_{1}+z_{3}}{n}\,\frac{z_{3}}{n}\right)^{rn}\enspace.$
(6)

∎

##  5 Proof of Lemma 7

We wish to show that at any extremum of $g_{r}$ we have $\alpha=\beta=1/2$. We
start by proving that at any such extremum $\alpha=\beta$. Note that since, by
symmetry, we are free to flip either or both colorings, we can restrict
ourselves to the case where $\alpha\leq 1/2$ and $\gamma\leq\alpha/2$.

Let $h(x_{1},x_{2},x_{3},x_{4})=\sum_{i}x_{i}\ln x_{i}$ denote the entropy
function, and let us define the shorthand $(\partial/\partial
x-\partial/\partial y)f$ for $\partial f/\partial x-\partial f/\partial y$.
Also, recall that $q=1-2^{1-k}$ and that $p(\alpha,\beta,\gamma)\equiv p$ is

$1-\alpha^{k}-(1-\alpha)^{k}-\beta^{k}-(1-\beta)^{k}+\gamma^{k}+(\alpha-\gamma)^{k}+(\beta-\gamma)^{k}+(1-\alpha-\beta+\gamma)^{k}.$

We will consider the gradient of $\ln g_{r}$ along a vector that increases
$\alpha$ while decreasing $\beta$. We see

$\displaystyle\left(\frac{\partial}{\partial\alpha}-\frac{\partial}{\partial\beta}\right)\ln
g_{r}(\alpha,\beta,\gamma)$ (7) $\displaystyle=$
$\displaystyle\left(\frac{\partial}{\partial\alpha}-\frac{\partial}{\partial\beta}\right)\Bigl{(}h\bigl{(}\gamma,\alpha-\gamma,\beta-\gamma,1-\alpha-\beta+\gamma\bigr{)}\,-\,\ln
4\,+\,r\,(\ln p-2\ln q)\Bigr{)}$ $\displaystyle=$
$\displaystyle-\ln(\alpha-\gamma)+\ln(\beta-\gamma)$
$\displaystyle+\frac{kr}{p}\,\Bigl{(}-\alpha^{k-1}+(1-\alpha)^{k-1}+(\alpha-\gamma)^{k-1}+\beta^{k-1}-(1-\beta)^{k-1}-(\beta-\gamma)^{k-1}\Bigr{)}$
$\displaystyle\equiv$ $\displaystyle\phi(\alpha)-\phi(\beta)\enspace,$

where

$\phi(x)\,=\,-\ln(x-\gamma)+\frac{kr}{p}\,\left(-x^{k-1}+(1-x)^{k-1}+(x-\gamma)^{k-1}\right)\enspace.$

Here we regard $p$ as a constant in the definition of $\phi(x)$.

Observe now that if $(\alpha,\beta,\gamma)$ is an extremum of $g_{r}$ then it
is also an extremum of $\ln g_{r}$. Therefore, it must be that
$(\partial/\partial\alpha)\ln g_{r}=(\partial/\partial\beta)\ln g_{r}=0$ at
$(\alpha,\beta,\gamma)$ which, by (7), implies $\phi(\alpha)=\phi(\beta)$.
This, in turn, implies $\alpha=\beta$ since $\phi(x)$ is monotonically
decreasing in the interval $\gamma<x<1$:

$\frac{{\rm d}\phi}{{\rm
d}x}\,=\,-\frac{1}{x-\gamma}-\frac{k(k-1)r}{p}\,\Bigl{(}x^{k-2}+(1-x)^{k-2}-(x-\gamma)^{k-2}\Bigr{)}\,<\,0\enspace.$

Next we wish to show that in fact $\alpha=\beta=1/2$. Setting $\alpha=\beta$,
we consider the gradient of $\ln g_{r}$ along a vector that increases $\alpha$
and $\gamma$ simultaneously (using a similar shorthand for $\partial
g/\partial\alpha+\partial g/\partial\gamma$):

$\displaystyle\left(\frac{\partial}{\partial\alpha}+\frac{\partial}{\partial\gamma}\right)\ln
g_{r}(\alpha,\alpha,\gamma)$ $\displaystyle=$
$\displaystyle\left(\frac{\partial}{\partial\alpha}+\frac{\partial}{\partial\gamma}\right)\Bigl{(}h\bigl{(}\gamma,\alpha-\gamma,\alpha-\gamma,1-2\alpha+\gamma\bigr{)}\,-\,\ln
4\,+\,r\,(\ln p-2\ln q)\Bigr{)}$ $\displaystyle=$
$\displaystyle-\ln\gamma+\ln(1-2\alpha+\gamma)$
$\displaystyle+\frac{kr}{p}\left(-2\alpha^{k-1}+2(1-\alpha)^{k-1}-(1-2\alpha+\gamma)^{k-1}+\gamma^{k-1}\right)$
$\displaystyle\equiv$ $\displaystyle\psi(\alpha)\enspace.$

Clearly, $\psi(\alpha)=0$ when $\alpha=1/2$. To show that $1/2$ is the only
such $\alpha$, we show that $\psi$ decreases monotonically with $\alpha$ by
showing that if $0<\alpha<1/2$ and $\gamma\leq\alpha/2$, all three terms below
are negative for $k\geq 3$.

$\displaystyle\frac{\partial\psi}{\partial\alpha}$ $\displaystyle=$
$\displaystyle\,-\frac{2}{1-2\alpha+\gamma}$
$\displaystyle+\,\frac{2k(k-1)r}{p}\left(-\alpha^{k-2}-(1-\alpha)^{k-2}+(1-2\alpha+\gamma)^{k-2}\right)$
$\displaystyle+\,\frac{2k^{2}r}{p^{2}}\times\left(\gamma^{k-1}-(1-2\alpha+\gamma)^{k-1}-2\alpha^{k-1}+2(1-\alpha)^{k-1}\right)$
$\displaystyle\;\;\;\;\times\,\left(-(\alpha-\gamma)^{k-1}+\alpha^{k-1}+(1-2\alpha+\gamma)^{k-1}-(1-\alpha)^{k-1}\right)\enspace.$

The first and second terms are negative since $1-\alpha>1-2\alpha+\gamma>0$,
implying $(1-\alpha)^{k-2}>(1-2\alpha+\gamma)^{k-2}$. The second factor of the
third term is positive since $f(z)=z^{k-1}$ is convex and
$(1-\alpha)-\alpha=(1-2\alpha+\gamma)-\gamma$ (the factor of 2 on the last two
terms only helps us since $1-\alpha\geq\alpha$). Similarly, the third factor
is negative since
$(1-\alpha)-(1-2\alpha+\gamma)=\alpha-\gamma\geq\alpha-(\alpha-\gamma)=\gamma$.

Thus, $\partial\psi/\partial\alpha<0$ and $\alpha=1/2$ is the unique solution
to $\psi(\alpha)=0$. Therefore, if $(\alpha,\alpha,\gamma)$ is an extremum of
$g_{r}$ we must have $\alpha=1/2$. ∎

##  6 Proof of Lemma 5

###  6.1 Preliminaries

We will use the following form of Stirling’s approximation for $n!$, valid for
$n>0$

$\sqrt{2\pi n}\,n^{n}\,{\rm
e}^{-n}\,\left(1+\frac{1}{12n}\right)\,<\,n!\,<\,\sqrt{2\pi n}\,n^{n}\,{\rm
e}^{-n}\,\left(1+\frac{1}{6n}\right)\enspace.$ (8)

We will also use the following crude lower bound for $n!$, valid for $n\geq 0$

$n!\geq(n/{\rm e})^{n}\enspace,$ (9)

using the convention $0^{0}\equiv 1$.

Let $z_{1},\ldots,z_{d}$ be such that $\sum_{i=1}^{d}z_{i}=n$. Let
$\zeta_{i}=z_{i}/n$. Let $\vec{\zeta}=(\zeta_{1},\ldots,\zeta_{d-1})$.

  * •

If $z_{i}>0$ for all $i$, then using the upper and lower bounds of (8) for
$n!$ and $z_{i}!$ respectively, and reducing the denominator further by
changing the factor $1+1/(12z_{i})$ to $1$ for $i\neq 1$, we get

$\displaystyle{n\choose z_{1},\cdots,z_{d}}$ $\displaystyle<$
$\displaystyle(2\pi
n)^{-(d-1)/2}\,\left(\prod_{i=1}^{d}\zeta_{i}^{\,-1/2}\right)\,\left(\prod_{i=1}^{d}\zeta_{i}^{\,-\zeta_{i}}\right)^{n}\times\frac{1+1/(6n)}{1+1/(12z_{1})}$
(10) $\displaystyle\leq$ $\displaystyle(2\pi
n)^{-(d-1)/2}\,\left(\prod_{i=1}^{d}\zeta_{i}^{\,-1/2}\right)\,\left(\prod_{i=1}^{d}\zeta_{i}^{\,-\zeta_{i}}\right)^{n}\enspace,$

where for (10) we assumed w.l.o.g. that $z_{1}\leq n/2$. Thus,

${n\choose z_{1},\cdots,z_{d}}\,f(z_{1}/n,\ldots,z_{d-1}/n)^{n}\,\leq\,(2\pi
n)^{-(d-1)/2}\,\left(\prod_{i=1}^{d}\zeta_{i}^{\,-1/2}\right)\,g(\vec{\zeta})^{n}\enspace.$
(11)

  * •

For any $z_{i}\geq 0$, the upper bound of (8) and (9) give

${n\choose z_{1},\cdots,z_{d}}\,<\,\frac{7}{6}\,\sqrt{2\pi
n}\,\left(\prod_{i=1}^{d}\zeta_{i}^{\,-\zeta_{i}}\right)^{n}\enspace,$

implying a cruder bound

${n\choose
z_{1},\cdots,z_{d}}\,f(z_{1}/n,\ldots,z_{d-1}/n)^{n}\,\leq\,\frac{7}{6}\,\sqrt{2\pi
n}\,g(\vec{\zeta})^{n}\enspace.$ (12)

###  6.2 The main proof

Our approach is a crude form of the Laplace method for asymptotic integrals
[9] which amounts to approximating functions near their peak as Gaussians.

We wish to approximate $g(\vec{\zeta})$ in the vicinity of
$\vec{\zeta}_{\max}$. We will do this by Taylor expanding $\ln g$, which is
analytic since $g$ is analytic and positive. Since $\ln g$ increases
monotonically with $g$, both $g$ and $\ln g$ are maximized at
$\vec{\zeta}_{\max}$. Furthermore, at $\vec{\zeta}_{\max}$ the matrix of
second derivatives of $\ln g$ is that of $g$ divided by a constant, since

$\left.\frac{\partial^{2}\ln
g}{\partial\zeta_{i}\,\partial\zeta_{j}}\right|_{\vec{\zeta}=\vec{\zeta}_{\max}}\\!\\!\\!=\,\frac{1}{g_{\max}}\frac{\partial^{2}g}{\partial\zeta_{i}\,\partial\zeta_{j}}-\frac{1}{g_{\max}^{2}}\frac{\partial
g}{\partial\zeta_{i}}\frac{\partial g}{\partial\zeta_{j}}$

and at $\vec{\zeta}_{\max}$ the first derivatives of $g$ are all zero.
Therefore, if the matrix of second derivatives of $g$ at $\vec{\zeta}_{\max}$
has nonzero determinant, so does the matrix of the second derivatives of $\ln
g$.

Note now that since the matrix of second derivatives is by definition
symmetric, it can be diagonalized, and its determinant is the product of its
eigenvalues. Therefore, if its determinant is nonzero, all its eigenvalues are
smaller than some $\lambda_{\max}<0$. Thus, Taylor expansion around
$\vec{\zeta}_{\max}$ gives

$\ln g(\vec{\zeta})\,\leq\,\ln
g_{\max}\,+\,\frac{1}{2}\lambda_{\max}\,|\vec{\zeta}-\vec{\zeta}_{\max}|^{2}\,+\,O(|\vec{\zeta}-\vec{\zeta}_{\max}|^{3})$

or, exponentiating to obtain $g$,

$g(\vec{\zeta})\,\leq\,g_{\max}\,\exp\left(\frac{1}{2}\lambda_{\max}|\vec{\zeta}-\vec{\zeta}_{\max}|^{2}\right)\times\bigl{(}1+O(|\vec{\zeta}-\vec{\zeta}_{\max}|^{3})\bigr{)}\enspace.$

Therefore, there is a ball of radius $\rho>0$ around $\vec{\zeta}_{\max}$ and
constants $Y>0$ and $g_{*}<g_{\max}$ such that

$\displaystyle\mbox{If
}|\vec{\zeta}-\vec{\zeta}_{\max}|\leq\rho,\;\;g(\vec{\zeta})$
$\displaystyle\leq$ $\displaystyle
g_{\max}\,\exp\left(-Y\,|\vec{\zeta}-\vec{\zeta}_{\max}|^{2}\right)\enspace,$
(13) $\displaystyle\mbox{If
}|\vec{\zeta}-\vec{\zeta}_{\max}|>\rho,\;\;g(\vec{\zeta})$ $\displaystyle\leq$
$\displaystyle g_{*}\enspace.$ (14)

We will separate $S$ into two sums, one inside the ball and one outside:

$\sum_{\vec{\zeta}\in
Z:\,|\vec{\zeta}-\vec{\zeta}_{\max}|\leq\rho}{n\choose\zeta_{1}n,\cdots,\zeta_{d}n}\,f(\vec{\zeta})^{n}\;\;\;+\;\sum_{\vec{\zeta}\in
Z:\,|\vec{\zeta}-\vec{\zeta}_{\max}|>\rho}{n\choose\zeta_{1}n,\cdots,\zeta_{d}n}\,f(\vec{\zeta})^{n}\enspace.$

For the terms inside the ball, first note that if
$|\vec{\zeta}-\vec{\zeta}_{\max}|\leq\rho$ then

$\prod_{i=1}^{d}\zeta_{i}^{\,-1/2}\leq W\mbox{ where
}W=\left(\min_{i}\zeta_{\max,i}-\rho\right)^{-d/2}\enspace.$

Then, since
$|\zeta-\vec{\zeta}_{\max}|^{2}=\sum_{i=1}^{d-1}(\zeta_{i}-\zeta_{\max,i})^{2}$,
using (11) and (13) we have

$\displaystyle\sum_{\vec{\zeta}\in
Z:\,|\vec{\zeta}-\vec{\zeta}_{\max}|\leq\rho}{n\choose\zeta_{1}n,\cdots,\zeta_{d}n}\,f(\vec{\zeta})^{n}$
$\displaystyle\leq$ $\displaystyle(2\pi
n)^{-(d-1)/2}\,W\,g_{\max}^{n}\quad\times\sum_{z_{1},\cdots,z_{d-1}=-\infty}^{\infty}\exp\left(-nY\,\sum_{i=1}^{d-1}(\zeta_{i}-\zeta_{\max,i})^{2}\right)$
$\displaystyle=$ $\displaystyle(2\pi
n)^{-(d-1)/2}\,W\,g_{\max}^{n}\times\;\prod_{i=1}^{d-1}\left(\sum_{z_{i}=-\infty}^{\infty}\exp\left(-nY\,(z_{i}/n-\zeta_{\max,i})^{2}\right)\right)\enspace.$

Now if a function $\phi(z)$ has a single peak, on either side of which it is
monotonic, we can replace its sum with its integral with an additive error at
most twice its largest term:

$\left|\sum_{z=-\infty}^{\infty}\phi(z)-\int_{-\infty}^{\infty}\phi(z)\,{\rm
d}z\,\right|\leq 2\max_{z}\phi(z)$

and so

$\displaystyle\sum_{z_{i}=-\infty}^{\infty}\exp\left(-nY(z_{i}/n-\zeta_{\max,i})^{2}\right)$
$\displaystyle\leq$ $\displaystyle
2+\int_{-\infty}^{\infty}\exp\left(-nY(z_{i}/n-\zeta_{\max,i})^{2}\right)\,{\rm
d}z$ $\displaystyle=$ $\displaystyle\sqrt{\pi n/Y}+2<\sqrt{2\pi n/Y}$

where the last inequality holds for sufficiently large $n$. Multiplying these
$d-1$ sums together gives

$\sum_{\vec{\zeta}\in
Z:\,|\vec{\zeta}-\vec{\zeta}_{\max}|\leq\rho}{n\choose\zeta_{1}n,\cdots,\zeta_{d}n}\,f(\vec{\zeta})^{n}\,\leq\,WY^{-(d-1)/2}\,g_{\max}^{n}\enspace.$
(15)

Outside the ball, we use (12), (14) and the fact that the entire sum has at
most $n^{d-1}$ terms to write

$\sum_{\vec{\zeta}\in
Z:\,|\vec{\zeta}-\vec{\zeta}_{\max}|>\rho}{n\choose\zeta_{1}n,\cdots,\zeta_{d}n}\,f(\vec{\zeta})^{n}\,\leq\,n^{d-1}\times\frac{7}{6}\,\sqrt{2\pi
n}\,g_{*}^{n}\,<\,g_{\max}^{n}$ (16)

where the last inequality holds for sufficiently large $n$. Combining (16) and
(15) gives

$S<(WY^{-(d-1)/2}+1)\,g_{\max}^{n}\equiv D\times g_{\max}^{n}$

which completes the proof. (We note that the constant $D$ can be optimized by
replacing our sums by integrals and using Laplace’s method [2, 9].)

##  7 Conclusions

We have shown that the second moment method yields a very sharp estimate of
the threshold for hypergraph 2-colorability. It allows us not only to close
the asymptotic gap between the previously known bounds but, in fact, to get
the threshold within a small additive constant. Yet:

$\bullet$ While the second moment method tells us that w.h.p. an exponential
number of 2-colorings exist for $r=\mathrm{\Theta}(2^{k})$, it tells us
nothing about how to find a single one of them efficiently. The possibility
that such colorings actually cannot be found efficiently is extremely
intriguing.

$\bullet$ While we have shown that the second moment method works really well,
we’d be hard pressed to say why. In particular, we do not have a criterion for
determining a constraint satisfaction problem’s amenability to the method. The
fact that the method fails spectacularly for random $k$-SAT suggests that,
perhaps, rather subtle forces are at play.

Naturally, one can always view the success of the second moment method in a
particular problem as an aposteriori indication that the satisfying solutions
of the problem are “largely uncorrelated”. This viewpoint, though, is hardly
predictive. (Yet, it might prove useful to the algorithmic question above).

The solution-symmetry shared by NAE $k$-SAT and hypergraph 2-colorability but
not by $k$-SAT, i.e. the property that the complement of a solution is also a
solution, explains why the method gives a nonzero lower bound for these two
problems (and why it fails for $k$-SAT). Yet symmetry alone does not explain
why the bound becomes essentially tight as $k$ grows. In any case, we hope
(and, worse, consider it natural) that an appropriate notion of symmetry is
present in many more problems.

## References

  * [1] Achlioptas, D., Kim, J.H., Krivelevich, M. and Tetali, P. Two-coloring random hypergraphs.  Random Structures Algorithms, 20(2):249–259, 2002.
  * [2] Achlioptas, D., and Moore, C. The asymptotic order of the random $k$-SAT threshold.  To appear in 43rd Annual Symposium on Foundations of Computer Science (Vancouver, BC, 2002).
  * [3] Alon, N., and Spencer, J. A note on coloring random $k$-sets.  Unpublished manuscript.
  * [4] Beck, J. On $3$-chromatic hypergraphs.  Discrete Math. 24, 2 (1978), 127–137.
  * [5] Berge, C. Hypergraphs.  North-Holland Publishing Co., Amsterdam, 1989.  Combinatorics of finite sets, Translated from the French.
  * [6] Bernstein, F. Zur theorie der trigonometrische reihen.  Leipz. Ber. 60, (1908), 325–328.
  * [7] Chao, M.-T., and Franco, J. Probabilistic analysis of a generalization of the unit-clause literal selection heuristics for the $k$-satisfiability problem.  Inform. Sci. 51, 3 (1990), 289–314.
  * [8] Chvátal, V., and Reed, B. Mick gets some (the odds are on his side).  In 33rd Annual Symposium on Foundations of Computer Science (Pittsburgh, PA, 1992). IEEE Comput. Soc. Press, Los Alamitos, CA, 1992, pp. 620–627.
  * [9] de Bruijn, N.G. Asymptotic Methods in Analysis. North-Holland, 1958.
  * [10] Erdős, P. On a combinatorial problem.  Nordisk Mat. Tidskr. 11 (1963), 5–10, 40.
  * [11] Erdős, P., and Lovász, L. Problems and results on $3$-chromatic hypergraphs and some related questions.  609–627. Colloq. Math. Soc. János Bolyai, Vol. 10.
  * [12] Flajolet, P., Knuth, D. E., and Pittel, B. The first cycles in an evolving graph.  Discrete Math. 75, 1-3 (1989), 167–215.  Graph theory and combinatorics (Cambridge, 1988).
  * [13] Friedgut, E. Necessary and sufficient conditions for sharp thresholds of graph properties, and the $k$-SAT problem.  J. Amer. Math. Soc. 12 (1999), 1017–1054.
  * [14] Frieze, A. and McDiarmid, C. Algorithmic theory of random graphs.  Random Struct. Alg. 10, (1997), 5–42.
  * [15] Karoński, M. and Łuczak, T. Random hypergraphs.  Combinatorics, Paul Erdős is eighty, Vol. 2 (Kesztheley, 1993), 283–293, Bolyai Soc. Math. Stud., 2, Janos Bolyai Math. Soc., Budapest, 1996\.
  * [16] Krivelevich, M. and Sudakov, B. The chromatic numbers of random hypergraphs.  Random Struct. Alg. 12, (1998), 381–403.
  * [17] Lovász, L. Coverings and coloring of hypergraphs.  Proceedings of the Fourth Southeastern Conference on Combinatorics, Graph Theory, and Computing (Boca Raton, Florida, 1973), 3–12.
  * [18] Motwani, R. and Raghavan, P. Randomized Algorithms.  Cambridge University Press, Cambridge, 1995.
  * [19] Radhakrishnan, J., and Srinivasan, A. Improved bounds and algorithms for hypergraph $2$-coloring.  Random Structures Algorithms 16, 1 (2000), 4–32.
  * [20] Schmidt-Pruzan, J., and Shamir, E., and Upfal, E. Random hypergraph coloring algorithms and the weak chromatic number.  Journal of Graph Theory 8, (1985), 347–362.