# Wigner description of massless Dirac plasmas

J. L. Figueiredo<EMAIL_ADDRESS>Instituto de Plasmas
e Fusão Nuclear, Lisboa, Portugal  Instituto Superior Técnico, Lisboa,
Portugal  J. P. Bizarro  Instituto de Plasmas e Fusão Nuclear, Lisboa,
Portugal  Instituto Superior Técnico, Lisboa, Portugal  H. Terças  Instituto
de Plasmas e Fusão Nuclear, Lisboa, Portugal  Instituto Superior Técnico,
Lisboa, Portugal

###### Abstract

We derive a quantum kinetic model describing the dynamics of graphene
electrons in phase space based on the Wigner-Moyal procedure. In order to take
into account the quantum nature of the carriers, we start from the low-energy
Schrödinger equation describing the mean-field wavefunction for both the
conduction and valence electrons. The equation of motion for the Wigner matrix
is established, with the Coulomb interaction being introduced self-
consistently (in the Hartree approximation). The long wavelength limit for the
plasmon dispersion relation is obtained, for both ungated and gated
situations. As an application, we derive the corresponding hydrodynamical
equations and add discuss the correct value of the effective hydrodynamic mass
of the carriers from first principles, an issue that is crucial in the
establishment of the correct hydrodynamics of Dirac electrons, thus paving the
stage towards a more comprehensive description of graphene plasmonics.

xxx ; xxx ; xxx ; xxx

##  I Introduction

In the recent years, graphene has been extensively studied due to its
spectacular optical, electronic and mechanical properties Castro Neto  _et
al._ (2009). In addiction to its two-dimensional (2D) nature, the elementary
electronic excitations are described by a Dirac-like dispersion in the low-
energy limit Hwang and Das Sarma (2007); Wallace (1947). The relativistic
nature of graphene electrons, resulting from the cone-like dispersion relation
near the Dirac points, makes it useful for transparent electronic devices,
ultra-sensitive photodetectors and other high-performance optoelectronic
structures Rodrigo  _et al._ (2015); Chen  _et al._ (2012); Zeng  _et al._
(2015). Furthermore, graphene possesses extremely high quantum efficiency for
light-matter interactions Koppens  _et al._ (2011). Moreover, the collective
oscillations of the electron and hole densities lead to the formation of
plasmons (or plasma waves) Agarwal  _et al._ (2018). In fact, plasmonic of 2D
materials, such as graphene and transition metal dichalcogenides (TMDCs) and
hexagonal boron nitride (hBN) Dai  _et al._ (2014); Lin  _et al._ (2009), is
nowadays a very prominent field of research Li  _et al._ (2017). Graphene-
based plasmonics finds a variety of applications, as the versatility of
graphene enables the manufacture of optical devices working in different
frequency ranges, namely in the terahertz (THz) and the infra-red domains
Grigorenko  _et al._ (2012). While metal plasmonics exhibit large Ohmic
losses, which limits their applicability to optical processing devices, doped
graphene emerges as an alternative. Its large conductivity, in part due to the
zero-mass character of the carriers, encloses a wide range of potential
applications, such as high-frequency nanoelectronics, nanomechanics,
transparent electrodes, and composite materials Geim and Novoselov (2007). For
this reason, the possibility of electric gating has been extensively studied
in graphene, allowing for the manipulation of the Fermi level Yang  _et al._
(2016). Recently, gating with a solid electrolyte allowed carrier
concentrations as large as $10^{14}\ \text{cm${}^{-2}$}$ to be achieved, which
results in a Fermi energy of $E_{F}\simeq 1$ eV, such that a modulation of
optical transmission in the visible spectrum is possible Wang  _et al._
(2008); Laboratory  _et al._ (2008). The potentiality of THz emission was also
pointed out Cosme and Terças (2020) as a possible application recently, by
making use of a graphene field-effect transistor (gFET), and controlling the
applied gate voltage and injected current.

From the theoretical point of view, a variety of techniques have been
developed to establish the dynamics of Dirac electrons and holes in graphene,
ranging from semiclassical hydrodynamical models Ryzhii  _et al._ (2007);
Bistritzer and MacDonald (2009); Svintsov  _et al._ (2012) to quantum
formalisms that involve collective Green’s functions, such as the time-
dependent Hartree-Fock approximation Schwengelbeck  _et al._ (2000); Jasiak
_et al._ (2009), or the time-dependent density functional theory Teperik  _et
al._ (2013). In one hand, it is often the case that, when going towards a
complete quantum description, cumbersome equations crop up, which are of very
reduced utility; on the other, the semiclassical approach is based on the
Vlasov equation, which is adequate in the case of high temperatures and low
electron-hole densities, but may fail to describe other important quantum
phenomena Dugaev and Katsnelson (2013).

In this work, we establish a kinetic formalism based on Wigner’s formulation
of quantum mechanics to study the Dirac electrons and holes in phase space
Wigner (1932). We start by constructing a Schrödinger equation for the
conduction and valence electrons, derived from a microscopic tight-binding
model, for the low-energy electrons. After moving to the electron-hole basis,
we derive a kinetic equation for the Wigner matrix components, which correctly
incorporates the pseudo-spin degrees of freedom. The interaction is introduced
self-consistently via the Hartree approximation, which obeys the Poisson
equation. With this in hand, we are able to obtain the plasmon dispersion
relation, recovering the usual result based on the random-phase approximation
(RPA). By performing averages over the phase-space distributions (more
precisely, by taking the moments of the Wigner equation), hydrodynamical
equations are obtained allowing for a fluid description of the Dirac particles
in graphene. As a consequence, we are able to derive microscopically the
effective mass of a fluid particle and relate it to the Drude mass, thus
contributing to the understanding of an important question in graphene
hydrodynamics. We also show that the classical limit corresponds to previous
results based on the Vlasov equation for the classical distribution function.

##  II Graphene preliminaries

The electronic dynamics can be captured starting with a general form for the
Hamiltonian, written in a second quantized fashion,

$\hat{H}=\sum_{s,s^{\prime}}\sum_{\bm{R},\bm{R^{\prime}}}\hat{u}_{s}^{{\dagger}}(\bm{R})\matrixelement{\hat{u}_{s},\bm{R}}{\hat{H}}{\hat{u}_{s^{\prime}},\bm{R^{\prime}}}\hat{u}_{s^{\prime}}(\bm{R^{\prime}}),$
(1)

where $\bm{R}/\bm{R^{\prime}}$ run over the real lattice, and $s/s^{\prime}$
over the two sublattices $A$ and $B$. Additionally,
$\hat{u}_{s}^{\dagger}(\bm{R})$ and $\hat{u}_{s}(\bm{R})$ denote the creation
and annihilation operators, respectively, for each lattice point $\bm{R}$ and
sublattice $s$. Moreover, the tight-binding approximation can be settled with
a proper restriction on the matrix elements
$\matrixelement{\hat{u}_{s},\bm{R}}{\hat{H}}{\hat{u}_{s^{\prime}},\bm{R^{\prime}}}$.
By allowing hopping only between nearest neighbors, we can set all matrix
elements to zero except the cases of
$\matrixelement{\hat{u}_{s},\bm{R}}{\hat{H}}{\hat{u}_{s^{\prime}},\bm{R}+\bm{\delta}_{i}}=-t(1-\delta_{ss^{\prime}})$,
where t $\simeq$ 2.97 eV is the hopping integral and $\bm{\delta}_{i}$ are the
nearest-neighbor vectors Wallace (1947). In momentum space, the Hamiltonian
reduces to

$\hat{H}=\sum_{\bm{q}}\ \bm{\varphi}_{\bm{q}}^{{\dagger}}\
\left(\begin{array}[]{cc}0&-t\Delta\\\ -t\Delta^{\ast}&0\end{array}\right)\
\bm{\varphi}_{\bm{q}},$ (2)

where $\bm{q}=(q_{x},q_{y})$ is the wave-vector,
$\bm{\varphi}_{\bm{q}}=(\hat{u}_{A\bm{q}},\hat{u}_{B\bm{q}})^{T}$ and
$\Delta(\bm{q})=\sum_{i}e^{-i\bm{q}\cdot\bm{\delta}_{i}}$. We can diagonalize
Eq. (2) by the means of a similarity transformation
$\bm{\Phi}(\bm{q},t)=S\bm{\varphi}(\bm{q},t)$, where
$\bm{\Phi}^{T}(\bm{q},t)=(\hat{c}_{\bm{q}},\hat{v}_{\bm{q}})$ and
$\hat{c}_{\bm{q}}$ and $\hat{v}_{\bm{q}}$ respectively denote conduction and
valence electron annihilation operators. Plugging in Eq. (2), we obtain

$\hat{H}=\sum_{\bm{q}}\ \bm{\Phi}_{\bm{q}}^{{\dagger}}\
\left(\begin{array}[]{cc}\hbar\omega(\bm{q})&0\\\
0&-\hbar\omega(\bm{q})\end{array}\right)\ \bm{\Phi}_{\bm{q}},$ (3)

where

$S=\frac{1}{\sqrt{2}}\left(\begin{array}[]{cc}-e^{-i\theta(\bm{q})}&1\\\
e^{i\theta(\bm{q})}&1\end{array}\right),$ (4)

$e^{i\theta(\bm{q})}=\sqrt{\Delta(\bm{q})/\Delta^{\ast}(\bm{q})}$, and

$\displaystyle\hbar\omega(\bm{q})=$ $\displaystyle\
t\sqrt{\cos\Big(\sqrt{3}dq_{x}/2\Big{missing})\cos\Big(3dq_{y}/2\Big{missing})+2\cos\Big(\sqrt{3}dq_{x}\Big{missing})+3}\
,$ (5)

is the single-particle dispersion relation. Above,
$d\simeq$1.5\text{\,}\mathrm{\SIUnitSymbolAngstrom}$$ is the carbon-carbon
distance, related to the lattice parameter $a$ as $a=\sqrt{3}d$. We can expand
Eq. (3) around the Dirac point $\bm{K}=(4\pi/3\sqrt{3}d,0)$ and obtain
Catarina  _et al._ (2019)

$\omega(\bm{q})\simeq v_{F}\absolutevalue{\bm{q}},$ (6)

where $v_{F}=3ta/(2h)\simeq 10^{6}\text{ ms${}^{-1}$}$ is the Fermi velocity.
As such, the low-energy Hamiltonian in Fourier space yields

$\hat{H}=\sum_{\bm{q}}\xi(\bm{q})\big{(}\hat{c}_{\bm{q}}^{\dagger}\hat{c}_{\bm{q}}-\hat{v}_{\bm{q}}^{\dagger}\hat{v}_{\bm{q}}\big{)},$
(7)

where $\xi(\bm{q})=\hbar v_{F}\absolutevalue{\bm{q}}$ is the single-particle
dispersion relation near the Dirac point.

By inspecting Eq. (7), we still face the problem of treating negative-energy
excitations. To circumvent it, we move to the electron-hole basis by defining
the (fermionic) hole creation and annihilation operators,
$\hat{h}^{\dagger}_{\bm{q}}=\hat{v}_{-\bm{q}}$ and
$\hat{h}_{\bm{q}}=\hat{v}_{-\bm{q}}^{\dagger}$, and by setting the vacuum
(Fermi sea) as

$|\Omega_{\rm F}\rangle=\prod_{q}\hat{v}_{q}^{\dagger}|0\rangle.$ (8)

The latter is related to the Fermi level as
$\hat{H}|\Omega_{F}\rangle=E_{F}|\Omega_{F}\rangle$, where $E_{F}$ is the
Fermi energy. In this new basis, the Hamiltonian in Eq. (7) reads

$\hat{H}=E_{F}\
+\sum_{\bm{q}}\xi(\bm{q})\big{(}\hat{c}_{\bm{q}}^{\dagger}\hat{c}_{\bm{q}}+\hat{h}_{\bm{q}}^{\dagger}\hat{h}_{\bm{q}}\big{)}.$
(9)

We now write the equation of motion for both electrons and holes,
$i\hbar\partial_{t}\bm{\psi}=\hat{H}\bm{\psi}$, where
$\bm{\psi}^{T}(\bm{q},t)=(\hat{c}_{\bm{q}},\hat{h}_{\bm{q}})$. In what
follows, we shall drop the constant $E_{F}$, for simplicity. Then, the
Hamiltonian $\hat{H}$ comprises kinetic and potential terms, accounting for
the Coulomb interaction, which leads to

$\displaystyle i\hbar\frac{\partial}{\partial t}\psi^{\alpha}(\bm{q},t)=\ $
$\displaystyle\xi(\bm{q})\psi^{\alpha}(\bm{q},t)$ $\displaystyle+\int
d\bm{q}^{\prime}\
V^{\alpha}(\bm{q^{\prime}},t)\psi^{\alpha}(\bm{q}-\bm{q^{\prime}},t).$ (10)

Here, $V^{\alpha}(\bm{q},t)=\mathcal{Q}^{\alpha}\phi(\bm{q},t)$,
$\mathcal{Q}^{\alpha}$ is the charge of the carries of band $\alpha$ and
$\phi(\bm{q},t)$ is the Fourier transform of the electrostatic potential. In
real space, $\phi$ satisfies the Poisson equation

$\bm{\nabla}^{2}\phi=-\frac{1}{\varepsilon_{0}\varepsilon_{r}}\sum_{\beta}\mathcal{Q}^{\beta}n^{\beta},$
(11)

with $n^{\beta}(\bm{r},t)=\absolutevalue{\psi^{\beta}(\bm{r},t)}^{2}$ denoting
the electron ($\beta=1$) and hole ($\beta=2$) densities. We perform our
calculations within the Hartree (mean-field) approximation, thus neglecting
exchange contributions (correlations). This procedure is valid for small
values of the coupling parameter $r_{s}\sim\alpha_{s}/\varepsilon_{r}$ (ratio
of the average potential energy to the average kinetic energy), where
$\alpha_{s}=e^{2}/(4\pi\varepsilon_{0}\hbar v_{F})\simeq 2.2$ is the graphene
structure constant. As a result, the Hartree approximation is reliable as long
as $\varepsilon_{r}\gtrsim 2.2$. A formal solution to Eq. (11) can be written
as

$\phi(\bm{r},t)=\frac{1}{4\pi\varepsilon_{0}\varepsilon_{r}}\int
d\bm{r}^{\prime}\
\sum\displaylimits_{\beta}\frac{\mathcal{Q}^{\beta}n^{\beta}(\bm{r^{\prime}},t)}{\absolutevalue{\bm{r}-\bm{r^{\prime}}}},$
(12)

thus providing a closure relation to the system of Eq. (10). In other words,
we have

$\phi(\bm{q},t)=\mathcal{U}(\bm{q})\sum\displaylimits_{\beta}\
s^{\beta}n^{\beta}(\bm{q},t),$ (13)

where
$\mathcal{U}(\bm{q})=e/2\varepsilon_{0}\varepsilon_{r}\absolutevalue{\bm{q}}$
is the 2D Fourier transform of the Coulomb potential, $e=\absolutevalue{e}$ is
the elementary charge and $s^{\beta}=2\beta-3$ is the sign of the
$\beta-$charge.

##  III Wigner formalism

A very handy way to treat the electronic system of (10) is by using the Wigner
picture of quantum mechanics Wigner (1932), which allows for a fully phase-
space description, in close analogy with the classical case. In the classical
limit, the Wigner function denotes the probability density of finding a
particle in a given infinitesimal phase space volume $d\bm{r}\ d\bm{k}$
centred in $(\bm{r},\bm{k})$. In the quantum case, due to the commutation
relation between $\bm{r}$ and $\bm{q}$, the Heisenberg uncertainty principle
prevents particles to localize in a specific phase-space point, and a proper
distribution function is not possible to construct Moyal (1949). However, we
can still construct the Wigner function $W(\bm{r},\bm{k},t)$ whose properties
are similar to those of a classical distribution function. It is a function of
both a spatial coordinate $\bm{r}$ and a wave-vector coordinate $\bm{k}$.
Although $\bm{r}$ and $\bm{k}$ are not conjugated to each other, they both
give information about the spatial and momentum distributions of the system,
which is described by a wave-function $\psi(\bm{r},t)$. The phase-space
regions where $W(\bm{r},\bm{k},t)$ takes negative values are purely quantum,
having no classical analogue as it cannot be univocally defined as a particle
phase-space density. For that reason, the Wigner function is often refereed to
as a quasi-distribution function. The mathematical definition of the Wigner
function is given by the Weyl transform Weyl (2002) of the density operator
$\hat{\rho}$. The expectation values of any operator can be computed by
integrating the Weyl transform of the same operator multiplied by
$W(\bm{r},\bm{k},t)$, very similarly to the classical case Hillery  _et al._
(1997). For the present case, it is given by a tensor
$W^{\alpha\gamma}(\bm{r},\bm{k},t)$, accounting for the pseudo-spin degrees of
freedom,

$\displaystyle W^{\alpha\gamma}(\bm{r},\bm{k},t)$
$\displaystyle=\int\frac{d\bm{s}}{(2\pi)^{2}}\ e^{i\bm{k}\cdot\bm{s}}$
$\displaystyle\times{\psi}^{\alpha}(\bm{r}-\bm{s}/2,t)\
{\psi^{\ast}}^{\gamma}(\bm{r}+\bm{s}/2,t).$ (14)

Although the diagonal is real, the off-diagonal elements can have non-zero
imaginary parts (despite being complex-conjugated to each other). Those
elements represent the density correlations between the two populations.
Notwithstanding, in the present mean-field approximation, the off-diagonal
elements decouple, and will thus be neglected. The potential in Eq. (13) can
be calculated self-consistently from the diagonal components as (no summation)

$\displaystyle n^{\alpha}(\bm{r},t)$ $\displaystyle=\int d\bm{k}\
W^{\alpha\alpha}(\bm{r},\bm{k},t)\ ,$ (15)

while in momentum space Eq. (14) transforms to

$W^{\alpha\gamma}(\bm{q},\bm{k},t)={\psi}^{\alpha}(\bm{k}+\bm{q}/2,t)\
{\psi^{\ast}}^{\gamma}(\bm{k}-\bm{q}/2,t).$ (16)

To construct an equation for the Wigner matrix, we start by writing Eq. (10)
and its hermitian conjugate for two independent momentum coordinates,
$\bm{q}_{1}$ and $\bm{q}_{2}$, and distinct pseudo-spin indices, $\alpha$ and
$\gamma$. By multiplying the fist by ${\psi^{\ast}}^{\gamma}(\bm{q}_{2},t)$
and the second by ${\psi}^{\alpha}(\bm{q}_{1},t)$ and subtracting the
equations, we obtain

$\displaystyle
i\hbar\frac{\partial}{\partial_{t}}\Big{[}{\psi}^{\alpha}(\bm{q}_{1},t){\psi^{\ast}}^{\gamma}(\bm{q}_{2},t)\Big{]}$
$\displaystyle=\Big{[}\xi(\bm{q}_{1})-\xi(\bm{q}_{2})\Big{]}{\psi}^{\alpha}(\bm{q}_{1},t){\psi^{\ast}}^{\gamma}(\bm{q}_{2},t)$
$\displaystyle+\int
d\bm{q}^{\prime}\Big{[}V^{\alpha}(\bm{q}^{\prime},t)\psi^{\alpha}(\bm{q}_{1}-\bm{q}^{\prime},t){\psi^{\ast}}^{\gamma}(\bm{q}_{2},t)-V^{\gamma}(\bm{q}^{\prime},t)\psi^{\alpha}(\bm{q}_{1},t){\psi^{\ast}}^{\gamma}(\bm{q}_{2}-\bm{q}^{\prime},t)\Big{]},$
(17)

where we used
$V^{\alpha}(\bm{q^{\prime}},t)^{\ast}=V^{\alpha}(-\bm{q^{\prime}},t)$.
Recalling the identity
$f(\bm{k}+\bm{q})=e^{\bm{q}\cdot\bm{\nabla}_{\bm{k}}}f(\bm{k})$, and
introducing the coordinate transformation $\bm{q}_{1}=\bm{k}+\bm{q}/2$,
$\bm{q}_{2}=\bm{k}-\bm{q}/2$ into Eq. (17), it follows, after some
straightforward algebra

$\displaystyle\Big{[}i\hbar\ \frac{\partial}{\partial
t}-\Delta\xi^{-}(\bm{q},\bm{k})\Big{]}W^{\alpha\gamma}(\bm{q},$
$\displaystyle\bm{k},t)=e\int d\bm{q}^{\prime}\ \phi(\bm{q}^{\prime},t)$
$\displaystyle\times\Delta W^{\alpha\gamma}(\bm{q},\bm{k},\bm{q^{\prime}},t),$
(18)

where $\Delta
W^{\alpha\gamma}(\bm{q},\bm{k},\bm{q^{\prime}},t)=s^{\alpha}W^{\alpha\gamma}(\bm{q}-\bm{q^{\prime}},\bm{k}-\bm{q^{\prime}}/2,t)-s^{\gamma}W^{\alpha\gamma}(\bm{q}-\bm{q^{\prime}},\bm{k}+\bm{q^{\prime}}/2,t)$
and $\Delta\xi^{-}(\bm{q},\bm{k})=\xi(\bm{k}+\bm{q}/2)-\xi(\bm{k}-\bm{q}/2)$.
Equations (18) and (13) define the Wigner-Poisson model for 2D Dirac particles
in Fourier space. For the diagonal components
$W^{\alpha\alpha}(\bm{r},\bm{k},t)\doteq W^{\alpha}$, the real space version
of Eq. (18) reads

$\displaystyle i\hbar\frac{\partial}{\partial
t}W^{\alpha}+i\hbar\mathcal{K}\\{W^{\alpha}\\}=\int d\bm{q}\
e^{i\bm{q}\cdot\bm{r}}\ \Big{(}W^{\alpha}_{-}-W^{\alpha}_{+}\Big{)}$
$\displaystyle\times V^{\alpha}(\bm{q},t),$ (19)

where $W^{\alpha}_{\pm}\doteq W^{\alpha}(\bm{r},\bm{k}\pm\bm{q}/2,t)$ and

$\mathcal{K}\\{W^{\alpha}\\}=v_{F}\int\frac{d\bm{r}^{\prime}}{(2\pi)^{2}}\frac{\sin\big(2\bm{k}\cdot\bm{r^{\prime}}\big{missing})}{\absolutevalue{\bm{r^{\prime}}}^{3}}\
W^{\alpha}(\bm{r}-\bm{r^{\prime}},\bm{k},t),$ (20)

represents the kinetic operator. Its nonlocal nature is profoundly related to
the relativistic nature of the Dirac electrons, and strongly differs from the
term appearing in semi-classical equations proposed in the literature,
$\mathcal{K}\\{W^{\alpha}\\}=\bm{v}_{\bm{k}}\cdot\bm{\nabla}W^{\alpha}$, with
$\hbar\bm{v}_{\bm{k}}=\partial_{\bm{k}}\xi(\bm{k})=\hbar
v_{F}\bm{k}/\absolutevalue{\bm{k}}$ (see Ref. Stauber  _et al._ (2007) and
references therein). It is convenient to compare it with the conventional
parabolic dispersion relation, $\xi(\bm{k})=\hbar^{2}\bm{k}^{2}/(2m)$, for
which the kinetic operator reduces to the usual (local) convective derivative
Hillery  _et al._ (1997)

$\mathcal{K}\\{W^{\alpha}\\}=\frac{\hbar\bm{k}}{m}\cdot\bm{\nabla}W^{\alpha}.$
(21)

Moreover, in the limit $\hbar\rightarrow 0$, we can approximate
$W_{\pm}^{\alpha}$ as

$W_{\pm}^{\alpha}\simeq
W^{\alpha}\pm\frac{\hbar\bm{q}}{2}\cdot\bm{\nabla}_{\bm{p}}W^{\alpha},$ (22)

with $\bm{p}=\hbar\bm{k}$, which allows us to write the ray-tracing version of
Eq. (19) as (no summation)

$\frac{\partial}{\partial
t}W^{\alpha}+\frac{\bm{p}}{m}\cdot\bm{\nabla}W^{\alpha}-\bm{\nabla}V^{\alpha}\cdot\bm{\nabla}_{\bm{p}}W^{\alpha}\simeq
0.$ (23)

This means that, for parabolic particles, the classical limit of the Wigner
equation corresponds to the Vlasov equation. However, no classical limit of
Eq. (20) exists (i.e. it does not converge as $\hbar\rightarrow 0$), which
then confers to Dirac particles a purely quantum nature.

We now consider small perturbations around an equilibrium configuration,
keeping the lowest order contributions to the Wigner components. As such, we
can write

$W^{\alpha\gamma}(\bm{r},\bm{k},t)\simeq
W_{0}^{\alpha}(\bm{k})\delta^{\alpha\gamma}+\tilde{W}^{\alpha\gamma}(\bm{r},\bm{k},t),$
(24)

with $|\tilde{W}^{\alpha\gamma}|\ll W^{\alpha}_{0}$, which in momentum space
reads $W^{\alpha\gamma}(\bm{q},\bm{k},t)\simeq
W_{0}^{\alpha}(\bm{k})\delta(\bm{q})\delta^{\alpha\gamma}+\tilde{W}^{\alpha\gamma}(\bm{q},\bm{k},t)$.
Similarly, the density and electrostatic potential will be perturbed as

$\displaystyle n^{\alpha}(\bm{q},t)$ $\displaystyle\simeq
n_{0}^{\alpha}\delta(\bm{q})+\tilde{n}^{\alpha}(\bm{q},t),$ (25)
$\displaystyle\phi(\bm{q},t)$
$\displaystyle\simeq\mathcal{U}(\bm{q})\sum\displaylimits_{\beta}\left[n_{0}\delta(\bm{q})+\tilde{n}^{\beta}(\bm{q},t)\right],$
(26)

where $n_{0}^{\alpha}=\int d\bm{k}W_{0}^{\alpha}(\bm{k})$ is the equilibrium
density of each band. We can show that the lowest order contributions of Eq.
(18) vanish. Moreover, the first order terms provide

$\displaystyle i\hbar\frac{\partial}{\partial
t}\tilde{W}^{\alpha\gamma}(\bm{q},\bm{k},t)$
$\displaystyle=\Delta\xi^{-}(\bm{q},\bm{k})\tilde{W}^{\alpha\gamma}(\bm{q},\bm{k},t)+\sum\displaylimits_{\beta}\tilde{n}^{\beta}(\bm{q},t)$
$\displaystyle\times\mathcal{Q}^{\alpha}\mathcal{U}(\bm{q})\Delta
W_{0}^{\alpha}(\bm{k},\bm{q})\delta^{\alpha\gamma},$ (27)

where $\Delta W_{0}^{\alpha}(\bm{k},\bm{q})=\
W_{0}^{\alpha}(\bm{k}-\bm{q}/2)-W_{0}^{\alpha}(\bm{k}+\bm{q}/2)$. Equation
(27) is formally equivalent to Kubo’s formula for the linear response of a
many-body system Fetter and Walecka (1971) and reproduces the features
contained in RPA Das Sarma and Hwang (2009); Wunsch  _et al._ (2006a). This
formalism is specially advantageous to describe the dynamics of electrons that
are far from equilibrium, such as the case of plasma instabilities, with the
configuration being solely defined by $W_{0}$.

In what follows, we consider the case of negatively doped graphene, with the
conduction band filled up to the Fermi level $E_{F}>0$. In momentum-energy
space, this is defined by the cone $E(\bm{k})=\hbar
v_{F}\absolutevalue{\bm{k}}\Theta(k_{F}-\absolutevalue{\bm{k}})$, where
$\Theta(x)$ is the Heaviside step function, and $k_{F}=E_{F}/(\hbar v_{F})$ is
the Fermi wave-number. The latter is related with the doping density $n_{0}$
by

$k_{F}=\sqrt{\frac{4\pi n_{0}}{g_{s}g_{v}}},$ (28)

where $g_{s}g_{v}$ accounts for the spin ($g_{s}=2$) and valley ($g_{v}=2$)
degeneracy. The first results from the degeneracy of the spin populations in
each energy band, which we have neglected in our treatment so far, and the
latter should be incorporated to consistently include the two minima in the
first Brillouin zone (BZ) Das Sarma  _et al._ (2009). Typical experimental
values of $n_{0}$ between $10^{9}-5\times 10^{12}~{}\text{cm}^{-2}$ are
achievable in graphene Adam  _et al._ (2007). In the case $E_{F}\gg k_{B}T$,
the presence of holes is negligible, which allows us to set $\tilde{W}^{22}=0$
and solve Eq. (18) for the Wigner component $\tilde{W}^{11}$ (notice that the
off-diagonal terms $W^{12}$ and $W^{21}$ are also not relevant in that case;
we also simplify the notation by dropping the superscripts, $W\doteq W^{11}$
and $n\doteq n_{1}$). After Fourier transforming the time coordinate, we can
recast Eq. (27) into a more familiar form,

$\tilde{W}(\bm{q},\bm{k},\omega)=\tilde{n}(\bm{q},t)\ \mathcal{U}(\bm{q})\
\frac{\Delta
W_{0}(\bm{k},\bm{q})}{\hbar\omega-\xi(\bm{k}+\bm{q}/2)+\xi(\bm{k}-\bm{q}/2)}.$
(29)

Upon integrating both sides in $\bm{k}$, we obtain the dielectric function

$\epsilon(\bm{q},\omega)=1-\mathcal{U}(\bm{q})\Pi(\bm{q},\omega),$ (30)

where the polarizability function $\Pi(\bm{q},\omega)$ is defined as

$\Pi(\bm{q},\omega)=\int d\bm{k}\
\frac{W_{0}(\bm{k})-W_{0}(\bm{k}+\bm{q})}{\hbar\omega+\xi(\bm{k})-\xi(\bm{k}+\bm{q})}.$
(31)

The plasmon dispersion relation is given by the zeros of
$\epsilon(\bm{q},\omega)$. We note that, in the absence of streaming,
$\epsilon(\mathbf{q},\omega)$ and $\omega(\mathbf{q})$ depend only on
$q\doteq|\mathbf{q}|$. The formal result of Eq. (30) is equivalent to the RPA
Das Sarma and Hwang (2009). In the long wavelength limit $q\rightarrow 0$, the
plasmon frequency can be computed using the noninteracting irreducible
polarizability, and Eq. (30) is recovered. We assume the equilibrium
configuration to be given by

$W_{0}(\bm{k})=n_{0}\Theta(k_{F}-\absolutevalue{\bm{k}})/\pi k_{F}^{2},$ (32)

and expanding Eq. (30) around $q=0$, keeping only terms up to
$\mathcal{O}(q^{2})$, the plasmon dispersion relation is obtained

$\omega^{2}=\omega_{p}^{2}\frac{q}{k_{F}}+\frac{3}{4}v_{F}^{2}q^{2},$ (33)

where $\omega_{p}$ is the characteristic plasmon frequency

$\omega_{p}=\Bigg{(}\frac{e^{2}n_{0}v_{F}}{2\hbar\varepsilon_{0}\varepsilon_{r}}\Bigg{)}^{1/2}.$
(34)

For the typical experimental values $\varepsilon_{r}=2.5$ and
$n_{0}\in[5\times 10^{-3},1]\times 10^{12}\ $\mathrm{c}\mathrm{m}^{-2}$$,
$\omega_{p}$ lies in the $\mathrm{THz}$-region, $\omega_{p}\in[2.6,37.3]\
$\mathrm{THz}$$. The first term $\omega\sim\sqrt{q}$ describes the long
wavelength signature of plasmons in two-dimensional electron gases (2DEG)
Kittel (1963); Eliasson and Liu (2018); Wunsch  _et al._ (2006b); Liu  _et
al._ (2008). The most notable difference, when compared to the characteristic
plasmon frequency in the 3-dimensional parabolic case,
$\omega_{p}^{3D}=\sqrt{e^{2}n_{0}/(\varepsilon_{0}\varepsilon_{r}m)}$, is the
appearance of $\hbar$ in leading order, revealing its pure quantum nature.
Therefore, no classical counterpart exists for the 2D Dirac plasma, which
confirms our previous observation.

The same kinetic approach can be used to describe graphene electrons in a
field effect transistor (FET) structure, i.e., placed between two metallic
contacts, source and drain, and controlled by a gate. The gate voltage is
related to the carrier density by Zhu  _et al._ (2009)

$U(\mathbf{r},t)=en(\mathbf{r},t)\left(\frac{1}{C_{g}}+\frac{1}{C_{q}}\right),$
(35)

where $C_{g}$/$C_{q}$ are, respectively, the gate and quantum capacitance. The
gate capacitance is given by $C_{g}=\varepsilon_{0}\varepsilon_{r}/d_{0}$,
where $d_{0}$ is the gate separation. The quantum capacitance $C_{q}$ reflects
the change in the potential with the band occupancy and is defined as
$C_{q}=e^{2}D(E)$, where
$D(E)=g_{s}g_{v}\absolutevalue{E}/(2\pi\hbar^{2}v_{F}^{2})$ is the density of
states. For the typical carrier densities we are interested in, $C_{g}\ll
C_{q}$, and so the second term in Eq. (35) can be neglected. To include this
effect in our model, one just needs to add $-eU(\mathbf{r},t)$ to the
effective potential, which amounts to the substitution
$\mathcal{U}(q)\rightarrow\mathcal{U}(q)+e^{2}d_{0}/(\varepsilon_{0}\varepsilon_{r})$.
Thus, the dispersion relation becomes

$\omega^{2}=\omega_{p}^{2}\frac{q}{k_{F}}+\left(S^{2}+\frac{3}{4}v_{F}^{2}\right)q^{2}+\mathcal{O}(q^{3}),$
(36)

where
$S^{2}=2e^{2}v_{F}d_{0}\sqrt{n_{0}}/(\hbar\varepsilon_{0}\varepsilon_{r}\sqrt{\pi})$
is the sound velocity of the electron fluid Cosme and Terças (2020). Equation
(36) is plotted in Fig. 1, alongside with Eqs. (6) and (33) for comparison. It
is patent that the gate induces a linear term $\omega\sim q$, which rapidly
dominates over the $\sqrt{q}$ term, given a typical value for $k_{F}d_{0}$ of
the order of unit. This is a consequence of the (static) screening produced by
the gate potential.

Figure 1: Positive branch of the plasmon dispersion relation in ungated (blue)
and gated (dashed red) configurations, along with the electron/hole dispersion
relation $\omega=v_{F}q$ (dotted brown) for $\varepsilon_{r}=2.5$ and
$d_{0}k_{F}=1$. Here, $q\doteq\absolutevalue{\bm{q}}$ is the wave vector
module.

##  IV Hydrodynamical model

One of the major advantages of the present description is the possibility of
calculating the moments out of Eq. (18), and subsequently construct
hydrodynamical models. Since we are concerned with the diagonal elements only,
we set $W^{\alpha}\doteq W^{\alpha\alpha}$. Similarly to the classical case,
we define the average value of an operator $\hat{\mathcal{G}}$ as

$\overline{\mathcal{G}}^{\alpha}=\frac{1}{n^{\alpha}(\bm{r},t)}\int d\bm{k}\
\mathcal{T}\\{\hat{\mathcal{G}}\\}W^{\alpha}(\bm{r},\bm{k},t),$ (37)

where $\bm{k}=(k_{x},k_{y})$ and $\mathcal{T}\\{\hat{\mathcal{G}}\\}$ denotes
the Weyl transform of the operator $\hat{\mathcal{G}}$ Weyl (2002),

$\mathcal{T}\\{\hat{\mathcal{G}}\\}=\int\frac{d\bm{s}}{(2\pi\hbar)^{d}}\
e^{i\bm{s}\cdot\bm{p}/\hbar}\
\matrixelement{\bm{r}-\bm{s}/2}{\hat{\mathcal{G}}}{\bm{r}+\bm{s}/2},$ (38)

The cases $\hat{\mathcal{G}}=1$ and $\hat{\bm{\mathcal{G}}}=\hat{\bm{k}}$
define the density and momentum variables,

$\displaystyle n^{\alpha}(\bm{r},t)$ $\displaystyle=\int d\bm{k}\
W^{\alpha}(\bm{r},\bm{k},t),$ (39)
$\displaystyle\overline{\bm{p}}^{\alpha}(\bm{r},t)$
$\displaystyle=\frac{1}{n^{\alpha}(\bm{r},t)}\int d\bm{k}\ \hbar\bm{k}\
W^{\alpha}(\bm{r},\bm{k},t).$ (40)

Above, we used $\mathcal{T}\\{\hat{\bm{k}}\\}=\bm{k}$ c (1). Equations (39)
and (40) are valid under the assumptions that the relevant quantities are
slow-varying in both space and time, which is the natural requirement to go
from a microscopic field theory to a macroscopic hydrodynamical model (coarse-
graining condition). The fluid description is valid if changes in the
macroscopic quantities take place on small spacial and temporal scales, i.e.
if $q\ll k_{F}$ and $\omega\ll\omega_{p}$, respectively. This requirement is
fulfilled if the characteristic time of the kinematic processes, $t\sim
1/\omega$, is much longer than the inverse collision frequency, $1/\nu_{c}$,
and the typical length $L$ is much greater than the mean free path $l\sim
v_{F}/\nu_{c}$, so that the plasma can be regarded locally as in a quasi-
equilibrium configuration. However, due to the quantum nature of the model, we
were able to capture the relativistic Dirac structure in a rigorous way.
Actually, in the context of 2D quantum plasmas achieved in semiconductor
structures, the De Broglie wavelength is replaced by the Thomas-Fermi
screening length $\lambda_{\rm
TF}=\sqrt{g_{s}g_{v}n_{0}}e^{2}/(\sqrt{4\pi}\epsilon\hbar v_{F})$. This is the
analogue of the classical Debye length in plasmas, and differs from the Fermi
wavelength defined above, $\lambda_{F}=2\pi/k_{F}$ (ranging as $\sim 10~{}{\rm
nm}-10\mu{\rm m}$ for typical graphene parameters). Consequently, the
hydrodynamical limit is expected to be valid provided the condition

$q\lambda_{\rm TF}\ll 1,\quad q\lambda_{F}\ll 1.$ (41)

In graphene, $\lambda_{F}/\lambda_{\rm TF}\simeq 3.2$, such that the later
condition is the most determinant of validity. Typical values of $k_{F}$ are
found between $10^{3}$ and $10^{6}$ $\text{cm}^{-1}$. At smaller wavelengths,
the microscopic structure becomes important, and the hydrodynamical
approximation no longer holds.

The time evolution equations for the quantities in Eqs. (39) and (40) are
derived in Appendix A, and provide

$\frac{\partial}{\partial t}n^{\alpha}=-\frac{2}{\hbar}\int d\bm{k}\ d\bm{q}\
ie^{i\bm{r}\cdot\bm{q}}\
\Bigg{[}\sinh(\frac{\bm{q}}{2}\cdot\bm{\nabla}_{\bm{k}})\xi(\bm{k})\Bigg{]}\
W^{\alpha}(\bm{q},\bm{k},t),$ (42) $\frac{\partial}{\partial
t}(n^{\alpha}\overline{\bm{p}}^{\alpha})=\ -2\int d\bm{k}\ d\bm{q}\
ie^{i\bm{r}\cdot\bm{q}}\ \bm{p}\
\Bigg{[}\sinh(\frac{\bm{q}}{2}\cdot\bm{\nabla}_{\bm{k}})\xi(\bm{k})\Bigg{]}\
W^{\alpha}(\bm{q},\bm{k},t)-\mathcal{Q}^{\alpha}n^{\alpha}\bm{\nabla}\phi.$
(43)

###  IV.1 Classical and semi-classical limits

In order to clarify the meaning of some of the terms in Eqs. (42) and (43), it
is instructive to derive the classical limit $\hbar\rightarrow 0$ (the semi-
classical limits are obtained by keeping higher orders of $\hbar$). Performing
the Taylor expansion of the $\sinh(\cdot)$ operator, we get

$\displaystyle\frac{\partial}{\partial
t}n^{\alpha}+\bm{\nabla}\cdot\overline{\bm{j}}^{\alpha}$ $\displaystyle=\
\hbar\int d\bm{k}\ \mathcal{N}\\{W^{\alpha}\big{\\}},$ (44)
$\displaystyle\frac{\partial}{\partial
t}(n^{\alpha}\overline{\bm{p}}^{\alpha}+\bm{\nabla}P^{\alpha}+\mathcal{Q}^{\alpha}n^{\alpha}\bm{\nabla}\phi=$
$\displaystyle\ \hbar\int d\bm{k}\ \bm{p}\ \mathcal{N}\\{W^{\alpha}\\},$ (45)

where $\overline{\bm{j}}^{\alpha}=n^{\alpha}\overline{\bm{v}}^{\alpha}$ is the
density current,

$\overline{\bm{v}}^{\alpha}=\frac{v_{F}}{n^{\alpha}}\int d\bm{k}\
\frac{\bm{k}}{\absolutevalue{\bm{k}}}\ W^{\alpha}(\bm{r},\bm{k},t)$ (46)

is the velocity field and $P^{\alpha}$ is the pressure-tensor

$P^{\alpha}(\bm{r},t)=v_{F}\int d\bm{k}\
\frac{1}{\absolutevalue{\bm{p}}}\left(\begin{array}[]{cc}p_{x}^{2}&p_{x}p_{y}\\\
p_{y}p_{x}&p_{y}^{2}\\\ \end{array}\right)W^{\alpha}(\bm{r},\bm{k},t).$ (47)

The kernel $\mathcal{N}$ introduces quantum mechanical contributions

$\displaystyle\mathcal{N}\\{W^{\alpha}\big{\\}}=\sum_{n=1}^{+\infty}\frac{2\hbar^{2n-1}}{i(2n+1)!}$
$\displaystyle\int d\bm{q}\ e^{i\bm{q}\cdot\bm{r}}\
W^{\alpha}(\bm{q},\bm{k},t)$
$\displaystyle\times\Big{(}\frac{\bm{q}}{2}\cdot\bm{\nabla}_{\bm{p}}\Big{)}^{2n+1}\xi(\bm{p}).$
(48)

Note that, in the case of systems with parabolic dispersion,
$\xi(\bm{p})=\bm{p}^{2}/(2m)$, $\mathcal{N}$ vanishes. However, in the
massless case, all-order derivatives of $\xi(\bm{p})$ exist. By letting
$\hbar\rightarrow 0$, we obtain the classical limit of Eqs. (44) and (45).

With the aim of obtaining semi-classical models explicitly, we consider the
first-order terms of $\mathcal{N}\\{W^{\alpha}\\}$. Casting
$\order{\hbar^{2}}$ terms in Eq. (48), we obtain the modified hydrodynamical
equations

$\displaystyle\frac{\partial}{\partial
t}n^{\alpha}+\bm{\nabla}\cdot\overline{\bm{j}}^{\alpha}=\frac{\hbar^{2}}{24}\Bigg{[}\frac{\partial^{3}}{\partial
x^{3}}\Big{(}n^{\alpha}\overline{J_{xxx}}^{\alpha}\Big{)}+\frac{\partial^{3}}{\partial
y^{3}}\Big{(}n^{\alpha}\overline{J_{yyy}}^{\alpha}\Big{)}$
$\displaystyle+3\frac{\partial^{2}}{\partial x^{2}}\frac{\partial}{\partial
y}\Big{(}n^{\alpha}\overline{J_{xxy}}^{\alpha}\Big{)}+3\frac{\partial^{2}}{\partial
y^{2}}\frac{\partial}{\partial
x}\Big{(}n^{\alpha}\overline{J_{yyx}}^{\alpha}\Big{)}\Bigg{]},$ (49)
$\displaystyle\frac{\partial}{\partial
t}(n^{\alpha}\overline{\bm{p}}^{\alpha})+\bm{\nabla}P^{\alpha}+\mathcal{Q}^{\alpha}n^{\alpha}\bm{\nabla}\phi=\frac{\hbar^{2}}{24}\Bigg{[}\frac{\partial^{3}}{\partial
x^{3}}\Big{(}n^{\alpha}\overline{\bm{T}_{xxx}}^{\alpha}\Big{)}$
$\displaystyle+\frac{\partial^{3}}{\partial
y^{3}}\Big{(}n^{\alpha}\overline{\bm{T}_{yyy}}^{\alpha}\Big{)}+3\frac{\partial^{2}}{\partial
x^{2}}\frac{\partial}{\partial
y}\Big{(}n^{\alpha}\overline{\bm{T}_{xxy}}^{\alpha}\Big{)}$
$\displaystyle+3\frac{\partial^{2}}{\partial y^{2}}\frac{\partial}{\partial
x}\Big{(}n^{\alpha}\overline{\bm{T}_{yyx}}^{\alpha}\Big{)}\Bigg{]},$ (50)

where $J_{ijl}$ and $\bm{T}_{ijl}$ are the dispersive tensors

$\displaystyle J_{ijl}$
$\displaystyle=v_{F}\Bigg{(}\frac{3p_{i}p_{j}p_{l}}{\absolutevalue{\bm{p}}^{5}}-\frac{\delta_{ij}p_{l}+\delta_{jl}p_{i}+\delta_{li}p_{j}}{\absolutevalue{\bm{p}}^{3}}\Bigg{)},$
(51) $\displaystyle\bm{T}_{ijl}$ $\displaystyle=\bm{p}J_{ijl}.$ (52)

We will provide explicit expressions for these quantities below and discuss
how they modify the plasmon dispersion relation.

###  IV.2 Mass transport

Having set the relevant transport equations, we shall move now to a more
detailed discussion concerning the averaged momentum and velocity fields. Once
again, we are interested in the electron case, for which we shall drop the
band index $\alpha$, for the sake of notation. As one can readily observe by
comparing Eqs. (40) and (46), a proportionality relation of the form relation
$\overline{\bm{p}}=m\overline{\bm{v}}$ does not hold for massless particles
Chaves  _et al._ (2017). However, by allowing a space and time dependence on
the mass, we are able to define an effective “mass tensor” as

$\displaystyle m_{ij}(\bm{r},t)$
$\displaystyle=\frac{\overline{p}_{i}(\bm{r},t)}{\overline{v}_{i}(\bm{r},t)}\delta_{ij},$
$\displaystyle=\frac{\hbar}{v_{F}}\frac{\displaystyle\int d\bm{k}\ k_{i}\
W(\bm{r},\bm{k},t)}{\displaystyle\int d\bm{k}\ k_{i}\
W(\bm{r},\bm{k},t)/\absolutevalue{\bm{k}}}\delta_{ij},$ (53)

and the corresponding “mass density tensor” $\rho_{ij}=nm_{ij}$. The meaning
of such fields should be clear: although the carriers have no mass (Dirac
particles), the fluid velocity and momentum fields can be related via an
effective mass, providing a measure of the inertia of a fluid element. Due to
rotational symmetry, we must have $m_{ij}(\bm{r},t)=m(\bm{r},t)\delta_{ij}$.

We now modify the equilibrium used in the last section in order to describe
adiabatic (hydrodynamic) changes in the system,

$W(\bm{r},\bm{k},t)=\frac{n_{0}}{\pi k_{F}^{2}}\Theta\left(k_{F}-n_{0}\
\absolutevalue{\frac{\bm{k}-\overline{\bm{k}}(\bm{r},t)}{n(\bm{r},t)}}\right).$
(54)

This particular form of the equilibrium is consistent with the definitions
(39) and (40). Without loss of generality, we can assume $\overline{k}_{y}=0$.
Plugging Eq. (54) into the $x-$component of Eq. (46), it leads to

$\displaystyle\overline{v}_{x}(\bm{r},t)$
$\displaystyle=\frac{\overline{p}_{x}(\bm{r},t)}{\mathcal{M}\gamma(\bm{r},t)},$
(55)

where $\mathcal{M}=\hbar k_{F}/v_{F}$ is the Drude mass and $\gamma$ is the
Lorentz-like factor,

$\displaystyle\gamma^{-1}(n,\overline{p}_{x})$ $\displaystyle=\
\frac{8}{\pi}\int\displaylimits_{0}^{1}dy\
\frac{\sqrt{1-y^{2}}}{\sqrt{f_{+}}+\sqrt{f_{-}}}\ ,$ (56) $\displaystyle
f_{\pm}=\Big{(}p^{\prime}\pm$ $\displaystyle
n^{\prime}\sqrt{1-y^{2}}\Big{)}^{2}+y^{2},\
p^{\prime}=\frac{\overline{p}_{x}}{p_{F}},\ n^{\prime}=\frac{n}{n_{0}}.$

As we can see, $\gamma$ depends on the position and time merely through the
hydrodynamical variables $n(\bm{r},t)$ and
$\overline{p}_{x}(\bm{r},t)=\hbar\overline{k}_{x}(\bm{r},t)$.

Although the integral in Eq. (56) has no analytic solution in general, it
admits amenable forms in the limiting cases of small and large momenta. We can
verify that
$\left(\partial\gamma/\partial\overline{p}_{x}\right)|_{\overline{p}_{x}=0}=0$,
thus for $\overline{p}_{x}/p_{F}\ll 1$, $\gamma$ becomes

$\gamma(n)\simeq\frac{\pi}{4}\
\frac{1-n^{2}/n_{0}^{2}}{K_{1}(1-n^{2}/n_{0}^{2})-K_{2}(1-n^{2}/n_{0}^{2})},$
(57)

where $K_{1}(x)=\int_{0}^{\pi/2}d\theta\ (1-x^{2}\sin^{2}\theta)^{-1/2}$ and
$K_{2}(x)=\int_{0}^{\pi/2}d\theta\ (1-x^{2}\sin^{2}\theta)^{1/2}$ are the
complete elliptic integrals of the first and second kind, respectively.
Accordingly, the quantities $\overline{p}_{x}$ and $\overline{v}_{x}$ become
proportional to each other, ensuring a relation of the form
$\overline{p}_{x}=m(n)\overline{v}_{x}$. In the case $n=n_{0}$,
$m=\mathcal{M}$, as $\gamma(n_{0})=1$. For general out-of-equilibrium
conditions with small values of $\overline{p}_{x}/p_{F}$, the mass differs
from the Drude mass, as $m(n)=\gamma(n)\mathcal{M}$, with $\gamma(n)$ depicted
in Fig. 2. We also plot the commonly adopted expression
$\gamma(n)=\sqrt{n/n_{0}}$, for comparison. This form of $\gamma$ leads to the
effective mass $m=\mathcal{M}\sqrt{n^{\prime}}$, which is based on the local
approximation for the Fermi wave-vector, $k_{F}=\sqrt{\pi
n_{0}}\rightarrow\sqrt{\pi n(\bm{r},t)}$ Castro Neto  _et al._ (2009);
Svintsov  _et al._ (2013); Chaves  _et al._ (2017). Equation (57) results in
significant changes, when compared to the latter, with a functional dependence
on $n$ which strongly differs from proposed $\sqrt{n}$ behaviour, specially in
the case of $n^{\prime}\gg 1$ (see Fig. 2).

Figure 2: Graphical representation of Eq. (57) (solid blue), as a function of
the normalised density $n^{\prime}=n/n_{0}$, evaluated in the limit
$\overline{p}_{x}/p_{F}\ll 1$, which verifies $\gamma(n^{\prime})\sim\pi
n^{\prime}/4$ for $n^{\prime}\gg 1$; $\gamma(n^{\prime})=\sqrt{n^{\prime}}$
(dashed red), which has been adopted in the literature.  Figure 3: Exact
relation between average velocity and momentum fields, for several values of
the normalised density $n^{\prime}=n/n_{0}$. Near the origin, we find a linear
relation between both variables, while for large momentum, the velocity
converges, in absolute value, to the Fermi velocity, becoming only a function
of the sign of $\overline{p}_{x}$.

For $\overline{p}_{x}/p_{F}\gg 1$, we find
$\gamma\simeq\absolutevalue{\overline{p}_{x}}/(v_{F}\mathcal{M})$. This
implies $\overline{v}_{x}\simeq{\rm sign}(\overline{p}_{x})v_{F}$, meaning
that, for large fluid momentum, the fluid velocity approaches the Fermi
velocity. In Fig. 3, we can observe that a linear relation between momentum
and velocity holds for small velocities only.

###  IV.3 Semi-classical dispersion relation

Next, we provide analytical expressions for the quantities appearing in Eqs.
(47), (51) and (52) - in terms of the previous hydrodynamical variables
$n^{\alpha}$ and $\overline{\bm{p}}^{\alpha}$ \- via Eq. (54). We restrict,
yet again, the variations to the $x-$direction, such that the relevant
quantities are $\overline{J_{xxx}}$, $\overline{T^{x}_{xxx}}$ and $P_{xx}$.
Considering

$n(x,t)\simeq
n_{0}+\tilde{n}(x,t),\quad\overline{p}_{x}(x,t)\simeq\overline{p}_{x,0}+\tilde{\overline{p}}_{x}(x,t).$
(58)

and neglecting second order terms, the relevant components featuring the
hydrodynamical equations read (see Appendix B)

$\displaystyle P_{xx}$ $\displaystyle\simeq
p_{F}v_{F}n_{0}\left(\frac{1}{3}+\frac{3}{4}\frac{\tilde{n}}{n_{0}}\right),$
(59) $\displaystyle\overline{J_{xxx}}$
$\displaystyle\simeq-\frac{3v_{F}}{4p_{F}^{3}}\tilde{\overline{p}}_{x},$ (60)
$\displaystyle\overline{T^{x}_{xxx}}$
$\displaystyle\simeq\frac{3v_{F}}{p_{F}}\left(\frac{1}{4}-\frac{1}{8}\frac{\tilde{n}}{n_{0}}\right),$
(61)

Keeping only the $\order{\tilde{n}}$ terms in Eq. (57), we find

$\gamma(\tilde{n})\simeq 1+\frac{3}{4}\frac{\tilde{n}}{n_{0}},$ (62)

which leads to linearized velocity

$\overline{v}_{x}\simeq\frac{\overline{p}_{x}}{\mathcal{M}}.$ (63)

By combining the latter results, we obtain

$\frac{\partial}{\partial t}\tilde{n}+\frac{\partial}{\partial
x}\Big{(}n_{0}\frac{\tilde{\overline{p}}_{x}}{\mathcal{M}}\Big{)}=-\frac{\hbar^{2}}{24}\frac{\partial^{3}}{\partial
x^{3}}\Big{(}n_{0}\frac{3v_{F}}{4p_{F}^{3}}\tilde{\overline{p}}_{x}\Big{)},$
(64) $\frac{\partial}{\partial
t}(n_{0}\tilde{\overline{p}}_{x})+\frac{\partial}{\partial
x}\Big{(}\frac{3}{4}p_{F}v_{F}\tilde{n}\Big{)}+\mathcal{Q}n_{0}\frac{\partial}{\partial
x}\tilde{\phi}=\frac{\hbar^{2}}{24}\frac{\partial^{3}}{\partial
x^{3}}\Big{(}\frac{3v_{F}}{8p_{F}}\tilde{n}\Big{)},$ (65)

where $\tilde{\phi}(\bm{r},t)$ is the perturbed electrostatic potential,

$\tilde{\phi}(\bm{r},t)=\frac{\mathcal{Q}}{4\pi\varepsilon_{0}\varepsilon_{r}}\int
d\bm{r}^{\prime}\
\frac{\tilde{n}(x,t)}{\absolutevalue{\bm{r}-\bm{r^{\prime}}}},$ (66)

and $\mathcal{Q}=e$. Fourier analysis of Eqs. (64) and (65) leads to the
modified plasmon dispersion

$\omega^{2}=\omega_{p}^{2}\frac{q}{k_{F}}+\frac{3}{4}v_{F}^{2}q^{2}-\frac{\omega_{p}^{2}}{32}\frac{q^{3}}{k_{F}^{3}}+\nu\frac{\hbar^{2}q^{4}}{4\mathcal{M}^{2}},$
(67)

where $\nu=-1/32$ is a numerical factor resulting from the contributions of
both $\overline{J_{xxx}}$ and $\overline{T^{x}_{xxx}}$. The two first terms in
the r.h.s are in agreement with Eq. (33), and the third terms is a classical
($\hbar$-independent) contribution to the plasmon dispersion. The fourth one
accounts for a quantum correction ($\sim\hbar^{2}$), and is a consequence of
the Bohm potential Lee and Jung (2017), which plays the role of a quantum
pressure. For Dirac plasmas, we find a negative sign in the quantum correction
(the Bohm pressure softens the plasmon mode), contrary to the what has been
reported for the case of parabolic particles, $\hbar^{2}q^{4}/(4m^{2})$
Akbari-Moghanjoughi (2012); Manfredi and Haas (2001). Since the absolute value
of $\nu$ is rather small, such correction only becomes important for
intermediate wave-number values, as can be seen in Fig. 4. The plasmon
dispersion relation of Eq. (67) indicates that the present quantum
hydrodynamical model is superior to the RPA calculation performed in previous
chapters, as the latter hides the effect of Bohm dispersion, which is known to
play an important role in the case of dense plasmas. Akbari-Moghanjoughi
(2012); Li (2010).

Figure 4: Positive branch of the semi-classical plasmon dispersion relation in
Eq. (67) (orange dashed), together with its classical counterpart, given by
Eq. (33) (blue), with $\varepsilon_{r}=2.5$.

##  V Conclusions and future work

Using the Wigner-Moyal description, we have derived a quantum hydrodynamical
model for a Dirac electron-hole plasma in single-layer graphene. The massless
nature of the quasi-particles was captured by the low-energy limit of the
many-body Hamiltonian, which served as the basis for the construction of the
Wigner matrix. Then, the equation of motion for the Wigner matrix elements in
the phase-space was established, with the interacting potential being given by
the Hartree approximation. The latter turns out to be an excellent
approximation for graphene, as the coupling constant is very small,
$r_{s}\sim\alpha_{s}/\varepsilon_{r}\ll 1$, and independent of the electron
density. Moreover, a closed set of generalised quantum hydrodynamical
equations was derived by taking the moments of the quantum kinetic equation.
We found an infinite number of $\hbar-$dependent terms arising in both the
continuity and force equations, as a consequence of the linear dispersion
relation of the carriers. As it is shown, in the conventional case of
parabolic systems (i.e. for which the single-particle dispersion relation is
$\sim p^{2}$), these contributions vanish. To the best of our knowledge, such
corrections have not been discussed so far in the literature. By neglecting
the quantum corrections, our model reduces to existing semi-classical models
for Dirac fermions Svintsov  _et al._ (2012, 2013); Ryzhii  _et al._ (2007).
Moreover, a nonlinear relation between the fluid velocity and momentum fields
have been put forward in order to construct a hydrodynamical mass. We found a
relation of the form
$\overline{\bm{p}}=\gamma(\bm{r},t)\mathcal{M}\overline{\bm{v}}$, where
$\mathcal{M}$ is the Drude mass and $\gamma(\bm{r},t)$ is a local Lorentz-like
factor. In the limiting cases of small and large averaged momentum,
$\gamma(\bm{r},t)$ is given analytically. Up to first order in both the
density and momentum fluctuations, we find
$\overline{\bm{p}}\simeq\mathcal{M}\overline{\bm{v}}$. We stress the fact that
a quantitative description for the hydrodynamical mass is a major issue in the
current research, therefore being one important result of the current kinetic
description. We have also computed higher-order corrections to the plasmon
dispersion relation ($\sim q^{4}$) by taking into account the first quantum
corrections to the hydrodynamical equations.

We expect our formalism to be particularly suited to describe out-of-
equilibrium quantum dynamics of Dirac electrons in a more systematic manner.
While the linear response can be indeed described within the framework RPA
(which we recover upon linearization of the Wigner equation), the nonlinear
regimes deserve a more complete description, which are captured by the kinetic
equation that we have proposed here. This is crucial in the description of
saturation in the late stages of dynamical instabilities, for example.
Moreover, (quantum) diffusive processes may be captured by adapting the
strategy of quasi-linear theories accounting for particle-wave interactions
Nicholson (1983); by working out the collision integral at different levels of
approximation, we may obtain quantum versions of the Fokker-Planck equation
for relativistic-like particles. Another interesting aspect of our work is the
contribution to the hydrodynamical modelling of graphene plasmas. With the ab-
initio establishment of a hydrodynamical mass, we are able to simulate
experimentally relevant configurations, such as hydrodynamical instabilities
in field-effect transistor THz emitters Cosme and Terças (2020). Moreover, at
a phenomenological level, our hydrodynamical model is particularly suited to
accommodate both shear and odd (or Hall) viscosity terms (while a more
microscopic, kinetic theory is more demanding, involving the details of the
processes leading to dissipation). The latter has been identified to play a
special role in the topology of acoustic waves Souslov  _et al._ (2019), and
we hope that plasmons may also feature such topological transitions, opening
the venue for a range of applications.

###### Acknowledgements.

H. T. acknowledges Fundação para a Ciência e Tecnologia (FCT-Portugal) through
the Contract No. CEECIND/00401/2018 and J. P. B. acnkowledges the financial
support from the Quantum Flagship Grant PhoQuS (Grant No. 820392) of the
European Union.

##  Appendix A Derivation of the hydrodynamical set of equations

In this section, we derive the hydrodynamical set of Eqs. (42) and (43). To do
that, we start by time differentiating Eqs. (39) and (40) (timed
$n^{\alpha}$), and replace $\partial_{t}W^{\alpha}(\bm{r},\bm{k},t)$ by $\int
d\bm{q}\ e^{i\bm{q}\cdot\bm{r}}\ \partial_{t}W^{\alpha}(\bm{q},\bm{k},t)$,
given that Eq. (18) only provides time derivatives of the Wigner matrix
components in Fourier space. We get

$\displaystyle\frac{\partial}{\partial t}n^{\alpha}=$
$\displaystyle-\frac{1}{\hbar}\int d\bm{k}\ d\bm{q}\ ie^{i\bm{r}\cdot\bm{q}}\
\Delta\xi^{-}(\bm{q},\bm{k})\ W^{\alpha}(\bm{q},\bm{k},t)-e\int d\bm{q}\
d\bm{q}^{\prime}\ ie^{i\bm{q}\cdot\bm{r}}\phi(\bm{q^{\prime}},t)\int d\bm{k}\
\Delta W^{\alpha\alpha}(\bm{q},\bm{k},\bm{q^{\prime}},t),$ (68)
$\displaystyle\frac{\partial}{\partial
t}\big{(}n^{\alpha}\overline{\bm{p}}^{\alpha}\big{)}=$ $\displaystyle-\int
d\bm{k}\ d\bm{q}\ ie^{i\bm{r}\cdot\bm{q}}\ \bm{k}\
\Delta\xi^{-}(\bm{q},\bm{k})\ W^{\alpha}(\bm{q},\bm{k},t)-e\int d\bm{q}\
d\bm{q}^{\prime}\ ie^{i\bm{q}\cdot\bm{r}}\ \phi(\bm{q^{\prime}},t)\int
d\bm{k}\ \bm{k}\ \Delta W^{\alpha\alpha}(\bm{q},\bm{k},\bm{q^{\prime}},t),$
(69)

where $n^{\alpha}\doteq n^{\alpha}(\bm{r},t)$ and
$\overline{\bm{p}}^{\alpha}\doteq\overline{\bm{p}}^{\alpha}(\bm{r},t)$. Now,
we use the following relations

$\displaystyle\Delta W^{\alpha\alpha}(\bm{q},\bm{k},\bm{q^{\prime}},t)=$
$\displaystyle-2s^{\alpha}e^{-\bm{q^{\prime}}\cdot\bm{\nabla}_{\bm{q}}}\sinh(\frac{\bm{q^{\prime}}}{2}\cdot\bm{\nabla}_{\bm{k}})W^{\alpha}(\bm{q},\bm{k},t),$
(70) $\displaystyle\Delta\xi^{-}(\bm{q},\bm{k})=$ $\displaystyle\
2\sinh(\frac{\bm{q}}{2}\cdot\bm{\nabla}_{\bm{k}})\xi(\bm{k}),$ (71)

into Eqs. (68) and (69), which leads to

$\displaystyle\frac{\partial}{\partial t}n^{\alpha}=$
$\displaystyle-\frac{2}{\hbar}\int d\bm{k}\ d\bm{q}\ ie^{i\bm{r}\cdot\bm{q}}\
\Bigg{[}\sinh(\frac{\bm{q}}{2}\cdot\bm{\nabla}_{\bm{k}})\xi(\bm{k})\Bigg{]}W^{\alpha}(\bm{q},\bm{k},t)$
$\displaystyle+2\mathcal{Q}^{\alpha}\int d\bm{q}\ d\bm{q}^{\prime}\
ie^{i\bm{q}\cdot\bm{r}}\
\phi(\bm{q^{\prime}},t)e^{-\bm{q^{\prime}}\cdot\bm{\nabla}_{\bm{q}}}\ \int
d\bm{k}\
\sinh(\frac{\bm{q^{\prime}}}{2}\cdot\bm{\nabla}_{\bm{k}})W^{\alpha}(\bm{q},\bm{k},t),$
(72) $\displaystyle\frac{\partial}{\partial
t}\big{(}n^{\alpha}\overline{\bm{p}}^{\alpha}\big{)}=$ $\displaystyle-2\int
d\bm{k}\ d\bm{q}\ ie^{i\bm{r}\cdot\bm{q}}\ \bm{k}\
\Bigg{[}\sinh(\frac{\bm{q}}{2}\cdot\bm{\nabla}_{\bm{k}})\xi(\bm{k})\Bigg{]}\
W^{\alpha}(\bm{q},\bm{k},t)$ $\displaystyle+2\mathcal{Q}^{\alpha}\int d\bm{q}\
d\bm{q}^{\prime}\ ie^{i\bm{q}\cdot\bm{r}}\
\phi(\bm{q^{\prime}},t)e^{-\bm{q^{\prime}}\cdot\bm{\nabla}_{\bm{q}}}\int
d\bm{k}\ \bm{k}\
\sinh(\frac{\bm{q^{\prime}}}{2}\cdot\bm{\nabla}_{\bm{k}})W^{\alpha}(\bm{q},\bm{k},t).$
(73)

In 2D, the differential operator
$\sinh(\frac{\bm{q^{\prime}}}{2}\cdot\bm{\nabla}_{\bm{k}})$ reads

$\displaystyle\sinh(\frac{\bm{q^{\prime}}}{2}\cdot\bm{\nabla}_{\bm{k}})$
$\displaystyle=\sum_{n=0}^{+\infty}\ \sum_{m=0}^{2n+1}\ \theta_{nm}\
\Big{(}\frac{q^{\prime}_{x}}{2}\frac{\partial}{\partial_{k_{x}}}\Big{)}^{m}\Big{(}\frac{q^{\prime}_{y}}{2}\frac{\partial}{\partial_{k_{y}}}\Big{)}^{2n+1-m}=\frac{\bm{q}^{\prime}}{2}\cdot\bm{\nabla}_{\bm{k}}+\mathcal{O}\left(\frac{\partial^{3}}{\partial
k_{i}^{3}}\right),$ (74)

where $\theta_{nm}=1/[m!(2n+1-m)!]$. By plugging Eq. (74) into the second term
on the r.h.s. of Eq. (72), we get a sum of terms, all of them at least linear
in the derivatives in $\bm{k}$. Once the remaining integrand is
$\bm{k}-$independent, we can perform partial integration in $\bm{k}$, which
reduces to surface contributions. Assuming that the Winger matrix elements and
all order $\bm{k}-$derivatives go to zero as $k_{i}\rightarrow\pm\infty$, we
readly obtain Eq. (42). By repeating the same procedure, a similar argument
allows to simplify the last term on the r.h.s. of Eq. (73). However, in this
case, the remaining integrand is linear in $\bm{k}$, which gives a non-trivial
contribution, coming from the first term of Eq. (74). Therefore, we obtain

$\displaystyle\frac{\partial}{\partial
t}\big{(}n^{\alpha}\overline{\bm{p}}^{\alpha}\big{)}=$ $\displaystyle\ -2\int
d\bm{k}\ d\bm{q}\ ie^{i\bm{r}\cdot\bm{q}}\ \bm{k}\
\Bigg{[}\sinh(\frac{\bm{q}}{2}\cdot\bm{\nabla}_{\bm{k}})\xi(\bm{k})\Bigg{]}W^{\alpha}(\bm{q},\bm{k},t)$
$\displaystyle-\mathcal{Q}^{\alpha}\int d\bm{q}\ d\bm{q}^{\prime}\
d\bm{r}^{\prime}\ ie^{i\bm{q}\cdot\bm{r}}\
\bm{q^{\prime}}\phi(\bm{q^{\prime}},t)\
e^{-\bm{q^{\prime}}\cdot\bm{\nabla}_{\bm{q}}}\
\frac{e^{-i\bm{q}\cdot\bm{r}^{\prime}}}{(2\pi)^{2}}.\
n^{\alpha}(\bm{r}^{\prime},t).$ (75)

Above, we used the property $\int d\bm{k}\
W^{\alpha}(\bm{q},\bm{k},t)=n^{\alpha}(\bm{q},t)$, which follows from Fourier
transforming Eq. (15), and further replaced $n^{\alpha}(\bm{q},t)$ by $\int
d\bm{r}^{\prime}\
e^{-i\bm{q}^{\prime}\cdot\bm{r}^{\prime}}n^{\alpha}(\bm{r}^{\prime},t)\
(2\pi)^{-2}$. Using the identities
$e^{-\bm{q^{\prime}}\cdot\bm{\nabla}_{\bm{q}}}e^{-i\bm{q}\cdot\bm{r}^{\prime}}=e^{i(\bm{q}^{\prime}-\bm{q})\bm{r}^{\prime}}$,
$i\bm{q}e^{i\bm{q}\cdot\bm{r}}=\bm{\nabla}e^{i\bm{q}\cdot\bm{r}}$ and $\int
d\bm{q}\
e^{i(\bm{r}-\bm{r}^{\prime})\bm{q}}(2\pi)^{-2}=\delta(\bm{r}-\bm{r}^{\prime})$,
we finally arrive at Eq. (43).

##  Appendix B Calculation of the dispersive tensors

Plugging Eq. (54) into the relevant components of Eqs. (47), (51) and (52)
gives

$\displaystyle P_{xx}$ $\displaystyle=\frac{\hbar v_{F}n_{0}}{\pi
k_{F}^{2}}\int\displaylimits_{0}^{k_{F}}dk_{y}\left[k_{+}(k_{+}^{2}+k_{y}^{2})^{1/2}-k_{-}(k_{-}^{2}+k_{y}^{2})^{1/2}-k_{y}^{2}\log(\frac{k_{+}(k_{-}^{2}+k_{y}^{2})^{1/2}}{k_{-}(k_{+}^{2}+k_{y}^{2})^{1/2}})\right].$
(76) $\displaystyle\overline{J_{xxx}}$ $\displaystyle=\frac{2v_{F}n_{0}}{\pi
p_{F}^{2}n}\int\displaylimits_{0}^{k_{F}}dk_{y}\
k_{y}^{2}\left[\frac{1}{(k_{+}^{2}+k_{y}^{2})^{1/2}}-\frac{1}{(k_{-}^{2}+k_{y}^{2})^{1/2}}\right],$
(77) $\displaystyle\overline{T^{x}_{xxx}}$ $\displaystyle=\frac{2\hbar
v_{F}n_{0}}{\pi
p_{F}^{2}n}\int\displaylimits_{0}^{k_{F}}dk_{y}\left[\frac{k_{-}^{3}}{(k_{-}^{2}+k_{y}^{2})^{3/2}}-\frac{k_{+}^{3}}{(k_{+}^{2}+k_{y}^{2})^{1/2}}\right],$
(78)

where $k_{\pm}=\overline{k}_{x}\pm n^{\prime}\sqrt{k_{F}^{2}-k_{y}^{2}}$ and
$n^{\prime}=n/n_{0}$. For the equilibrium configuration, we expect to have
$\overline{p}_{x,0}=0$, such that only the first order term contributes.
Expanding Eqs. (76)$-$(78) and neglecting
$\order{\tilde{\overline{p}}_{x}^{2}/p_{F}^{2}}$ terms, it leads to

$\displaystyle\overline{J_{xxx}}$ $\displaystyle\simeq\frac{4v_{F}}{\pi
p_{F}^{3}}\frac{n_{0}/n}{(n^{2}/n_{0}^{2}-1)^{2}}\Bigg{[}2K_{1}\Big{(}1-\frac{n_{0}^{2}}{n^{2}}\Big{)}-\Big{(}1+\frac{n^{2}}{n_{0}^{2}}\Big{)}K_{2}\Big{(}1-\frac{n_{0}^{2}}{n^{2}}\Big{)}\Bigg{]}\tilde{\overline{p}}_{x},$
(79) $\displaystyle\overline{T^{x}_{xxx}}$
$\displaystyle\simeq\frac{4v_{F}}{\pi
p_{F}}\frac{n/n_{0}}{(n^{2}/n_{0}^{2}-1)^{2}}\Bigg{[}2K_{1}\Big{(}1-\frac{n_{0}^{2}}{n^{2}}\Big{)}-\Big{(}1+\frac{n^{2}}{n_{0}^{2}}\Big{)}K_{2}\Big{(}1-\frac{n_{0}^{2}}{n^{2}}\Big{)}\Bigg{]},$
(80) $\displaystyle P_{xx}$
$\displaystyle\simeq\frac{4p_{F}v_{F}n_{0}}{3\pi}\frac{1}{n^{2}/n_{0}^{2}-1}\Bigg{[}\frac{n^{4}}{n_{0}^{4}}K_{2}\Big{(}1-\frac{n_{0}^{2}}{n^{2}}\Big{)}-\frac{n^{2}}{n_{0}^{2}}K_{1}\Big{(}1-\frac{n_{0}^{2}}{n^{2}}\Big{)}\Bigg{]}.$
(81)

Keeping only the $\order{n}$ terms in Eqs. (79)-(81), it yields Eqs.
(59)-(61).

## References

  * Castro Neto  _et al._ (2009) A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 (2009).
  * Hwang and Das Sarma (2007) E. H. Hwang and S. Das Sarma, Phys. Rev. B 75, 205418 (2007).
  * Wallace (1947) P. R. Wallace, Phys. Rev. 71, 622 (1947).
  * Rodrigo  _et al._ (2015) D. Rodrigo, O. Limaj, D. Janner, D. Etezadi, F. J. Garcia de Abajo, V. Pruneri, and H. Altug, Science 349, 165–168 (2015).
  * Chen  _et al._ (2012) J. Chen, M. Badioli, P. Alonso-González, S. Thongrattanasiri, F. Huth, J. Osmond, M. Spasenović, A. Centeno, A. Pesquera, P. Godignon, and et al., Nature 487, 77–81 (2012).
  * Zeng  _et al._ (2015) S. Zeng, K. V. Sreekanth, J. Shang, T. Yu, C.-K. Chen, F. Yin, D. Baillargeat, P. Coquet, H.-P. Ho, A. V. Kabashin, and K.-T. Yong, Advanced Materials 27, 6163 (2015).
  * Koppens  _et al._ (2011) F. H. L. Koppens, D. E. Chang, and F. J. García de Abajo, Nano Letters 11, 3370 (2011), pMID: 21766812.
  * Agarwal  _et al._ (2018) A. Agarwal, M. Vitiello, L. Viti, A. Cupolillo, and A. Politano, Nanoscale 10 (2018), 10.1039/C8NR01395K.
  * Dai  _et al._ (2014) S. Dai, Z. Fei, Q. Ma, A. S. Rodin, M. Wagner, A. S. McLeod, M. K. Liu, W. Gannett, W. Regan, K. Watanabe, T. Taniguchi, M. Thiemens, G. Dominguez, A. H. C. Neto, A. Zettl, F. Keilmann, P. Jarillo-Herrero, M. M. Fogler, and D. N. Basov, Science 343, 1125 (2014).
  * Lin  _et al._ (2009) Y. Lin, T. V. Williams, and J. W. Connell, The Journal of Physical Chemistry Letters 1, 277 (2009).
  * Li  _et al._ (2017) Y. Li, Z. Li, C. Chi, H. Shan, L. Zheng, and Z. Fang, Advanced Science 4, 1600430 (2017).
  * Grigorenko  _et al._ (2012) A. N. Grigorenko, M. Polini, and K. S. Novoselov, Nature Photonics 6, 749 (2012).
  * Geim and Novoselov (2007) A. Geim and K. Novoselov, Nature materials 6, 183 (2007).
  * Yang  _et al._ (2016) Y. Yang, Z. Shi, J. Li, and Z.-Y. Li, Photonics Research 4, 65 (2016).
  * Wang  _et al._ (2008) F. Wang, Y. Zhang, C. Tian, C. Girit, A. Zettl, M. Crommie, and Y. R. Shen, Science 320, 206 (2008).
  * Laboratory  _et al._ (2008) L. B. N. Laboratory, U. S. D. of Energy. Office of Scientific, and T. Information,  _Dirac Charge Dynamics in Graphene by Infrared Spectroscopy_ (Lawrence Berkeley National Laboratory, 2008).
  * Cosme and Terças (2020) P. Cosme and H. Terças, ACS Photonics 7, 1375 (2020).
  * Ryzhii  _et al._ (2007) V. Ryzhii, A. Satou, and T. Otsuji, Journal of Applied Physics - J APPL PHYS 101 (2007), 10.1063/1.2426904.
  * Bistritzer and MacDonald (2009) R. Bistritzer and A. H. MacDonald, Phys. Rev. B 80, 085109 (2009).
  * Svintsov  _et al._ (2012) D. Svintsov, V. Vyurkov, S. Yurchenko, T. Otsuji, and V. Ryzhii, Journal of Applied Physics 111, 083715 (2012).
  * Schwengelbeck  _et al._ (2000) U. Schwengelbeck, L. Plaja, L. Roso, and E. C. Jarque, Journal of Physics B: Atomic, Molecular and Optical Physics 33, 1653 (2000).
  * Jasiak  _et al._ (2009) R. Jasiak, G. Manfredi, P.-A. Hervieux, and M. Haefele, New Journal of Physics 11, 063042 (2009).
  * Teperik  _et al._ (2013) T. V. Teperik, P. Nordlander, J. Aizpurua, and A. G. Borisov, Phys. Rev. Lett. 110, 263901 (2013).
  * Dugaev and Katsnelson (2013) V. K. Dugaev and M. I. Katsnelson, Phys. Rev. B 88, 235432 (2013).
  * Wigner (1932) E. Wigner, Phys. Rev. 40, 749 (1932).
  * Catarina  _et al._ (2019) G. Catarina, B. Amorim, E. V. Castro, J. M. V. P. Lopes, J. M. V. P. Lopes, and N. Peres, “Twisted bilayer graphene: Low-energy physics, electronic and optical properties,” in  _Handbook of Graphene Set_ (John Wiley & Sons, Ltd, 2019) Chap. 6, pp. 177–231.
  * Moyal (1949) J. E. Moyal, Proc. Cambridge Phil. Soc. 45, 99 (1949).
  * Weyl (2002) H. Weyl, Z. Phys. 46 (2002).
  * Hillery  _et al._ (1997) M. Hillery, R. F. O’Connell, M. O. Scully, and E. P. Wigner, “Distribution functions in physics: Fundamentals,” in  _Part I: Physical Chemistry. Part II: Solid State Physics_ (Springer Berlin Heidelberg, Berlin, Heidelberg, 1997) pp. 273–317.
  * Stauber  _et al._ (2007) T. Stauber, N. M. R. Peres, and F. Guinea, Phys. Rev. B 76, 205423 (2007).
  * Fetter and Walecka (1971) A. L. Fetter and J. D. Walecka,  _Quantum Theory of Many-Particle Systems_ (McGraw-Hill, Boston, 1971).
  * Das Sarma and Hwang (2009) S. Das Sarma and E. H. Hwang, Physical Review Letters 102 (2009), 10.1103/physrevlett.102.206412.
  * Wunsch  _et al._ (2006a) B. Wunsch, T. Stauber, F. Sols, and F. Guinea, New Journal of Physics 8, 318 (2006a).
  * Das Sarma  _et al._ (2009) S. Das Sarma, E. Hwang, and Q. Li, Phys. Rev. B 80 (2009), 10.1103/PhysRevB.80.121303.
  * Adam  _et al._ (2007) S. Adam, E. H. Hwang, V. M. Galitski, and S. Das Sarma, Proceedings of the National Academy of Sciences 104, 18392–18397 (2007).
  * Kittel (1963) C. Kittel,  _Quantum Theory of Solids_ (John Wiley and Sons, Inc, 1963).
  * Eliasson and Liu (2018) B. Eliasson and C. S. Liu, Physics of Plasmas 25, 012105 (2018).
  * Wunsch  _et al._ (2006b) B. Wunsch, T. Stauber, F. Sols, and F. Guinea, New Journal of Physics 8, 318 (2006b).
  * Liu  _et al._ (2008) Y. Liu, R. Willis, K. Emtsev, and T. Seyller, Phys. Rev. B 78 (2008), 10.1103/PhysRevB.78.201403.
  * Zhu  _et al._ (2009) W. Zhu, V. Perebeinos, M. Freitag, and P. Avouris, Phys. Rev. B 80, 235402 (2009).
  * c (1) _This can be proved directly from the definition in ( 38), after performing the coordinate transformation $\bm{s}=\bm{r}-\bm{z}/2$, together with the relation $\bra{\bm{r}}\hat{\bm{k}}\ket{\bm{r^{\prime}}}=-i\bm{\nabla}\ \delta(\bm{r}-\bm{r^{\prime}})$_.
  * Chaves  _et al._ (2017) A. J. Chaves, N. M. R. Peres, G. Smirnov, and N. A. Mortensen, Physical Review B 96 (2017), 10.1103/physrevb.96.195438.
  * Svintsov  _et al._ (2013) D. Svintsov, V. Vyurkov, V. Ryzhii, and T. Otsuji, Physical Review B 88 (2013), 10.1103/physrevb.88.245444.
  * Lee and Jung (2017) M.-J. Lee and Y.-D. Jung, Physics Letters A 381, 636 (2017).
  * Akbari-Moghanjoughi (2012) M. Akbari-Moghanjoughi, Physics of Plasmas 19, 042701 (2012), https://doi.org/10.1063/1.3699535 .
  * Manfredi and Haas (2001) G. Manfredi and F. Haas, Phys. Rev. B 64, 075316 (2001).
  * Li (2010) S.-C. Li, Physics of Plasmas 17, 082307 (2010).
  * Nicholson (1983) D. Nicholson,  _Introduction to Plasma Theory_ (Wiley, 1983).
  * Souslov  _et al._ (2019) A. Souslov, K. Dasbiswas, M. Fruchart, S. Vaikuntanathan, and V. Vitelli, Phys. Rev. Lett. 122, 128001 (2019).


# Obliquely propagating ion-acoustic shock waves in degenerate quantum plasma

M.K. Islam1,∗, S. Biswas1,∗∗, N.A. Chowdhury2,∗∗∗, A. Mannan1,†, M.
Salahuddin1,‡, and A.A. Mamun1,§ 1 Department of Physics, Jahangirnagar
University, Savar, Dhaka-1342, Bangladesh
2Plasma Physics Division, Atomic Energy Centre, Dhaka-1000, Bangladesh
e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

A theoretical investigation has been carried out on the propagation of
nonlinear ion-acoustic shock waves (IASHWs) in a collsionless magnetized
degenerate quantum plasma system composed of inertial non-relativistic
positively charged light and heavy ions, inertialess ultra-relativistically
degenerate electrons and positrons. The reductive perturbation method has been
employed to drive the Burgers’ equation. It has been observed that under
consideration, our plasma model supports only positive potential shock
structure. It is also found that the amplitude and steepness of the IASHWs
have been significantly modified by the variation of ion kinematic viscosity,
oblique angle, number density, and charge state of the plasma species. The
results of our present investigation will be helpful for understanding the
propagation of IASHWs in white dwarfs and neutron stars.

###### keywords:

Shock waves; Degenerate quantum plasma; White dwarfs; Neutron stars.

††journal: “Contributions to Plasma Physics”

##  1 Introduction

The existence of light (viz., ${}^{1}_{1}{\mbox{H}}$ [1, 2, 3, 4],
${}^{4}_{2}{{\mbox{He}}}$ [1, 5, 6, 7, 8], ${}^{12}_{6}{\mbox{C}}$ [7, 8], and
${}^{16}_{8}{\mbox{O}}$ [7, 8, 9], etc.) and heavy (viz.,
${}^{56}_{26}{\mbox{Fe}}$ [10, 11], ${}^{85}_{37}{\mbox{Rb}}$ [11, 12], and
${}^{96}_{42}{\mbox{Mo}}$ [11, 12], etc.) ions in astrophysical compact
objects (viz., white dwarfs and neutron stars) has received a substantial
attention to investigate ion-acoustic (IA) waves (IAWs) in degenerate quantum
plasma (DQP). The particle number density in the white dwarfs (i.e., $6\times
10^{29}\mbox{cm}^{-3}$) and neutron stars (i.e., $6\times
10^{36}\mbox{cm}^{-3}$) is extremely high [1, 9, 11, 13]. This extreme number
density of these degenerate particles dictates to follow the Heisenberg’s
uncertainty principle, and according to the uncertainty principle when the
position of these particles is confined then the momentum of these particles
tends to very large. This excessive momentum leads to generate extreme outward
degenerate pressure which is counter-balanced by the inward gravitational
compression. The dynamics of these particles is mathematically modeled under
two categories, namely, non-relativistic and ultra-relativistic limits by
Chandrasekhar [1]. The pressure of non-relativistic light and heavy ions is
expressed as $P_{j}=K_{j}N_{j}^{\alpha}$; where $\alpha=5/3$,
$K_{j}=3\pi\hbar^{2}/5m_{j}$, $j=l$ for light ion, $j=h$ for heavy ion, and
$\hbar$ is the Planck constant [14, 15, 16]. The degenerate pressure of ultra-
relativistic electrons and positrons can be expresses as
$P_{s}=K_{s}N_{s}^{\gamma_{s}}$; where $s=e$ for electron, $s=p$ for positron,
$\gamma_{s}=4/3$, $K_{s}=3\hbar c/4$, and $c$ is the speed of light [14, 15,
16].

The presence of positrons in white dwarfs and neutron stars has been
extensively discussed in the Refs. [17, 18, 19, 20]. Sultana and Schlickeiser
[16] investigated IA solitary waves in DQP containing degenerate electrons,
light ions, and inertial mobile non-degenerate heavy ions. Gill et al. [21]
investigated the IA shock waves (IASHWs) in relativistic DQP composed of
electrons, positrons and ions, and found that the height of the potential is
maximum for the lower positron density. Ata-ur-Rahman et al. [22] considered
an unmagnetized DQP containing inertial ions, and inertialess electrons and
positrons, and demonstrated that the amplitude of the IAWs decreases with the
increase of positron number density. Hossen et al. [23] studied IASHWs in a
four-component plasma system having inertialess electrons and positrons, and
inertial heavy and light ions, and reported that the amplitude of the shock
profile decreases with positron number density. Mamun and Shukla [24]
investigated electrostatic solitary waves propagating in ultra-relativistic
plasma medium consisting of degenerate electrons and cold mobile ions, and
reported that the wave amplitude increases with the increase of ion number
density.

The strong magnetic field (i.e., about 1 Mega Gauss) in white dwarfs was
predicted by Blackett [25] and observed by Zeeman spectroscopy [26, 27]. El-
Taibany et al. [28] analyzed solitary waves in a magnetized degenerate
electron-positron plasma, and found that the wave amplitude increases with the
oblique angle which is the angle between the external magnetic field and the
direction of wave propagation. Shaukat [29] investigated IA solitary waves in
the presence of magnetic field, and observed that the solitary wave amplitude
increases with increasing oblique angle.

The energy dissipation of the shock wave, which is governed by the Burgers’
equation [30], may arise due to the kinematic viscosity of the medium. Hafez
et al. [31] studied IASHWs in weakly relativistic plasma containing electrons,
positrons, and ions, and noticed that the steepness of the IASHWs decreases
with the increase in the value of viscosity of plasma species but the
amplitude of the IASHWs does not change. Abdelwahed et al. [32] analyzed
IASHWs in a pair-ion plasma, and also found that the shock steepness decreases
with increasing ion viscosity.

Recently, Saini et al. [33] investigated heavy nucleus acoustic periodic waves
in DQP. Haider [34] examined the shock profiles in the presence of degenerate
inertial ions, and inertialess electrons and positrons. Atteya et al. [35]
studied IASHWs in a DQP which contains ion fluids, degenerate electrons, and
stationary heavy ions. To the best of authors’ knowledge, still no one
investigated IASHWs in a magnetized DQP having inertialess ultra-
relativistically degenerate electrons and positrons, and inertial positively
charged non-relativistic light and heavy ions. In this manuscript, we will
derive the Burgers’ equation, and will also use associated solution to examine
the basic features of IASHW in DQP.

The manuscript is organized in the preceding way: The governing equations are
described in section 2. The Burgers’ equation and associated shock solution
are presented in section 3. The results and discussion are presented in
section 4. A brief conclusion is presented in section 5.

##  2 Governing Equations

We consider a magnetized DQP system consisting of inertial positively charged
light ion (mass $m_{l}$; charge $eZ_{l}$; number density $N_{l}$), positively
charged heavy ion (mass $m_{h}$; charge $eZ_{h}$, number density $N_{h}$),
inertialess electron (mass $m_{e}$; charge $-e$, number density $N_{e}$), and
positron (mass $m_{p}$; charge $e$; number density $N_{p}$); where $Z_{l}$
($Z_{h}$) is the charge state of the light (heavy) ion. An uniform external
magnetic field $\mathbf{B}$ is existing in the direction of $z$-axis
($\mathbf{B}={B_{0}}\hat{z}$ and $\hat{z}$ is the unit vector). The
propagation of IAWs in DQP system is governed by the following equations:

$\displaystyle\frac{\partial N_{h}}{\partial
T}+\tilde{\nabla}\cdot(N_{h}{U}_{h})=0,$ (1) $\displaystyle\frac{\partial
U_{h}}{\partial
T}+({U}_{h}\cdot\tilde{\nabla}){U}_{h}=-\frac{Z_{h}e}{m_{h}}\nabla\tilde{\Phi}+\frac{Z_{h}eB_{0}}{m_{h}}({U}_{h}\times\hat{z})$
$\displaystyle\hskip
48.36958pt-\frac{1}{m_{h}N_{h}}\tilde{\nabla}P_{h}+\tilde{\eta}_{h}\tilde{\nabla}^{2}{U}_{h},$
(2) $\displaystyle\frac{\partial N_{l}}{\partial
T}+\tilde{\nabla}\cdot(N_{l}{U}_{l})=0,$ (3)
$\displaystyle\frac{\partial{U}_{l}}{\partial
T}+({U}_{l}\cdot\tilde{\nabla}){U}_{l}=-\frac{Z_{l}e}{m_{l}}\tilde{\nabla}\tilde{\Phi}+\frac{Z_{l}eB_{0}}{m_{l}}({U}_{l}\times\hat{z})$
$\displaystyle\hskip
43.2483pt-\frac{1}{m_{l}N_{l}}\tilde{\nabla}P_{l}+\tilde{\eta}_{l}\tilde{\nabla}^{2}{U}_{l},$
(4) $\displaystyle\tilde{\nabla}^{2}\tilde{\Phi}=4\pi
e(N_{e}-N_{p}-Z_{l}N_{l}-Z_{h}N_{h}),\ $ (5)

where ${U}_{l}$ (${U}_{h}$) is the fluid speed of light (heavy) ion;
$\tilde{\Phi}$ is the electrostatic wave potential; $P_{l}$ ($P_{h}$) is the
pressure for light (heavy) ion; $\tilde{\eta}_{l}=\mu/m_{l}N_{l}$
($\tilde{\eta}_{h}=\mu/m_{h}N_{h}$) is the kinematic viscosity for light
(heavy) ion. The degenerate pressure equations for electrons and positrons can
be expressed, respectively, as

$\displaystyle\tilde{\nabla}\tilde{\Phi}-\frac{1}{eN_{e}}\tilde{\nabla}P_{e}=0,$
(6)
$\displaystyle\tilde{\nabla}\tilde{\Phi}+\frac{1}{eN_{p}}\tilde{\nabla}P_{p}=0.\
$ (7)

Now, we have introduced the normalizing parameters: $n_{h}\rightarrow
N_{h}/n_{h0}$; $n_{l}\rightarrow N_{l}/n_{l0}$; $n_{e}\rightarrow
N_{e}/n_{e0}$; $n_{p}\rightarrow N_{p}/n_{p0}$; $u_{h}\rightarrow
U_{h}/C_{h}$; $u_{l}\rightarrow U_{l}/C_{h}$; $\phi\rightarrow
e\tilde{\Phi}/m_{e}c^{2}$; $t\rightarrow\omega_{ph}T$;
$\nabla\rightarrow\lambda_{Dh}^{-1}\tilde{\nabla}$;
$\eta\rightarrow\tilde{\eta}/\omega_{ph}\lambda_{Dh}^{2}$ [where IAWs speed
$C_{h}=(Z_{h}m_{e}c^{2}/m_{h})^{1/2}$; plasma frequency $\omega_{ph}=(4\pi
Z_{h}^{2}e^{2}n_{h0}/m_{h})^{1/2}$; the Debye length
$\lambda_{Dh}=(m_{e}c^{2}/4\pi Z_{h}e^{2}n_{h0})^{1/2}$, and for simplicity we
have considered $\eta=\eta_{l}=\eta_{h}$]. At equilibrium, the charge
neutrality condition can be written as
$n_{e0}=n_{p0}+Z_{l}n_{l0}+Z_{h}n_{h0}$. By using these normalizing
parameters, Eqs. (1)-(5) can be expressed in the normalized form

$\displaystyle\frac{\partial n_{h}}{\partial t}+\nabla\cdot(n_{h}u_{h})=0,$
(8) $\displaystyle\frac{\partial}{\partial
t}(u_{h})+(u_{h}\cdot\nabla)u_{h}=-\nabla\phi+\Omega_{ch}(u_{h}\times\hat{z})$
$\displaystyle\hskip 51.21504pt-\frac{\mu_{1}K_{1}}{n_{h}}\nabla
n_{h}^{\alpha}+\eta\nabla^{2}{u_{h}},$ (9) $\displaystyle\frac{\partial
n_{l}}{\partial t}+\nabla\cdot(n_{l}u_{l})=0,$ (10)
$\displaystyle\frac{\partial}{\partial
t}(u_{l})+(u_{l}\cdot\nabla)u_{l}=-\mu_{2}\nabla\phi+\mu_{2}\Omega_{ch}(u_{l}\times\hat{z})$
$\displaystyle\hskip 45.52458pt-\frac{\mu_{3}K_{2}}{n_{l}}\nabla
n_{l}^{\alpha}+\eta\nabla^{2}{u_{l}},$ (11)
$\displaystyle\nabla^{2}\phi=(1+\mu_{4}+\mu_{5})n_{e}-\mu_{4}n_{p}-\mu_{5}n_{l}-n_{h},\
$ (12)

where $\Omega_{ch}=\omega_{ch}/\omega_{ph}$, $\mu_{1}=1/Z_{h}$,
$\mu_{2}=Z_{l}m_{h}/Z_{h}m_{l}$, $\mu_{3}=m_{h}/Z_{h}m_{e}$,
$\mu_{4}=n_{p0}/Z_{h}n_{h0}$, $\mu_{5}=Z_{l}n_{l0}/Z_{h}n_{h0}$,
$K_{1}={n_{h0}}^{\alpha-1}K_{h}/m_{e}c^{2}$, and
$K_{2}={n_{l0}}^{\alpha-1}K_{l}/m_{e}c^{2}$. Now, by normalizing and
integrating Eqs. (6) and (7), the number densities of the inertialess
electrons and positrons can be obtained in terms of electrostatic potential
$\phi$, respectively, as

$\displaystyle
n_{e}=\bigg{[}1+\frac{\gamma_{e}-1}{\gamma_{e}K_{3}}\phi\bigg{]}^{\frac{1}{\gamma_{e}-1}},$
(13) $\displaystyle
n_{p}=\bigg{[}1-\frac{\gamma_{p}-1}{\gamma_{p}K_{4}}\phi\bigg{]}^{\frac{1}{\gamma_{p}-1}},\
$ (14)

where $K_{3}={n_{e0}}^{\gamma_{e}-1}K_{e}/m_{e}c^{2}$ and
$K_{4}={n_{p0}}^{\gamma_{p}-1}K_{p}/m_{e}c^{2}$. By expanding the right hand
side of Eqs. (13) and (14) up to second order in $\phi$, and substituting in
Eq. (12), we get

$\displaystyle\nabla^{2}\phi+\mu_{5}n_{l}+n_{h}=1+\mu_{5}+\beta_{1}\phi+\beta_{2}\phi^{2}+\cdot\cdot\cdot,\
$ (15)

where

$\displaystyle\beta_{1}=\frac{[\gamma_{p}K_{4}(1+\mu_{4}+\mu_{5})+\mu_{4}\gamma_{e}K_{3}]}{\gamma_{e}\gamma_{p}K_{3}K_{4}},$
$\displaystyle\beta_{2}=\frac{[(\gamma_{p}K_{4})^{2}(1+\mu_{4}+\mu_{5})(2-\gamma_{e})-\mu_{4}(\gamma_{e}K_{3})^{2}(2-\gamma_{p})]}{2(\gamma_{e}\gamma_{p}K_{3}K_{4})^{2}}.\
$

##  3 Derivation of the Burgers’ Equation

To study IASHWs, we derive Burgers’ equation by employing the reductive
perturbation method (RPM) [36, 37], and the stretched coordinates for
independent variables can be written as [38, 39, 40]

$\displaystyle\xi=\epsilon(l_{x}x+l_{y}y+l_{z}z-v_{p}t),$ (16)
$\displaystyle\tau=\epsilon^{2}t,\ $ (17)

where $v_{p}$ is the phase speed and $\epsilon$ is a smallness parameter
measuring the weakness of the dissipation (0$<$$\epsilon$$<$1). The $l_{x}$,
$l_{y}$, and $l_{z}$ are the directional cosines of $\boldsymbol{k}$ (wave
vector) along $x$, $y$, and $z$-axes, respectively (i.e.,
$l_{x}^{2}+l_{y}^{2}+l_{z}^{2}=1$). Then, the dependent variables can be
expressed in power series of $\epsilon$ as [39]

$\displaystyle n_{h}=1+\epsilon
n_{h}^{(1)}+\epsilon^{2}n_{h}^{(2)}+\cdot\cdot\cdot,$ (18) $\displaystyle
n_{l}=1+\epsilon n_{l}^{(1)}+\epsilon^{2}n_{l}^{(2)}+\cdot\cdot\cdot,$ (19)
$\displaystyle
u_{hx,y}=\epsilon^{2}u_{hx,y}^{(1)}+\epsilon^{3}u_{hx,y}^{(2)}+\cdot\cdot\cdot,$
(20) $\displaystyle
u_{lx,y}=\epsilon^{2}u_{lx,y}^{(1)}+\epsilon^{3}u_{lx,y}^{(2)}+\cdot\cdot\cdot,$
(21) $\displaystyle u_{hz}=\epsilon
u_{hz}^{(1)}+\epsilon^{2}u_{hz}^{(2)}+\cdot\cdot\cdot,$ (22) $\displaystyle
u_{lz}=\epsilon u_{lz}^{(1)}+\epsilon^{2}u_{lz}^{(2)}+\cdot\cdot\cdot,$ (23)
$\displaystyle\phi=\epsilon\phi^{(1)}+\epsilon^{2}\phi^{(2)}+\cdot\cdot\cdot.\
$ (24)

Now, by substituting Eqs. (16)-(24) into Eqs. (8)-(11) and (15), and
collecting the terms containing $\epsilon$, the first-order equations reduce
to

$\displaystyle
n_{h}^{(1)}=\frac{l_{z}^{2}}{\big{(}v_{p}^{2}-\alpha\mu_{1}l_{z}^{2}K_{1}\big{)}}\phi^{(1)},$
(25) $\displaystyle
u_{hz}^{(1)}=\frac{v_{p}l_{z}}{\big{(}v_{p}^{2}-\alpha\mu_{1}l_{z}^{2}K_{1}\big{)}}\phi^{(1)},$
(26) $\displaystyle
n_{l}^{(1)}=\frac{\mu_{2}l_{z}^{2}}{\big{(}v_{p}^{2}-\alpha\mu_{3}l_{z}^{2}K_{2}\big{)}}\phi^{(1)},$
(27) $\displaystyle
u_{lz}^{(1)}=\frac{\mu_{2}v_{p}l_{z}}{\big{(}v_{p}^{2}-\alpha\mu_{3}l_{z}^{2}K_{2}\big{)}}\phi^{(1)}.\
$ (28)

Now, the phase speed of IASHWs can be written as

$\displaystyle
v_{p}=v_{p+}=l_{z}\sqrt{\frac{a+\sqrt{a^{2}-4\beta_{1}b}}{2\beta_{1}}},$ (29)
$\displaystyle
v_{p}=v_{p-}=l_{z}\sqrt{\frac{a-\sqrt{a^{2}-4\beta_{1}b}}{2\beta_{1}}},\ $
(30)

where
$a=\alpha\mu_{1}\beta_{1}K_{1}+\alpha\mu_{3}\beta_{1}K_{2}+\mu_{2}\mu_{5}+1$
and
$b=\alpha^{2}\mu_{1}\mu_{3}\beta_{1}K_{1}K_{2}-\alpha\mu_{1}\mu_{2}\mu_{5}K_{1}-\alpha\mu_{3}K_{2}$.
The $x$ and $y$-components of the first-order momentum equations can be
manifested as

$\displaystyle
u_{hx}^{(1)}=-\frac{l_{y}v_{p}^{2}}{\Omega_{ch}(v_{p}^{2}-\alpha\mu_{1}l_{z}^{2}K_{1})}\frac{\partial\phi^{(1)}}{\partial\xi},$
(31) $\displaystyle
u_{hy}^{(1)}=\frac{l_{x}v_{p}^{2}}{\Omega_{ch}(v_{p}^{2}-\alpha\mu_{1}l_{z}^{2}K_{1})}\frac{\partial\phi^{(1)}}{\partial\xi},$
(32) $\displaystyle
u_{lx}^{(1)}=-\frac{l_{y}v_{p}^{2}}{\Omega_{ch}(v_{p}^{2}-\alpha\mu_{3}l_{z}^{2}K_{2})}\frac{\partial\phi^{(1)}}{\partial\xi},$
(33) $\displaystyle
u_{ly}^{(1)}=\frac{l_{x}v_{p}^{2}}{\Omega_{ch}(v_{p}^{2}-\alpha\mu_{3}l_{z}^{2}K_{2})}\frac{\partial\phi^{(1)}}{\partial\xi}.\
$ (34)

Now, by taking the next higher-order terms, the equation of continuity,
momentum equation, and Poisson’s equation can be written as

$\displaystyle\frac{\partial n_{h}^{(1)}}{\partial\tau}-v_{p}\frac{\partial
n_{h}^{(2)}}{\partial\xi}+l_{x}\frac{\partial
u_{hx}^{(1)}}{\partial\xi}+l_{y}\frac{\partial
u_{hy}^{(1)}}{\partial\xi}+l_{z}\frac{\partial u_{hz}^{(2)}}{\partial\xi}$
$\displaystyle+l_{z}\frac{\partial}{\partial\xi}\Big{(}n_{h}^{(1)}u_{hz}^{(1)}\Big{)}=0,$
(35) $\displaystyle\frac{\partial
u_{hz}^{(1)}}{\partial\tau}-v_{p}\frac{\partial
u_{hz}^{(2)}}{\partial\xi}+l_{z}u_{hz}^{(1)}\frac{\partial
u_{hz}^{(1)}}{\partial\xi}+\mu_{1}l_{z}\frac{\partial\phi^{(2)}}{\partial\xi}$
$\displaystyle+\alpha\mu_{1}l_{z}K_{1}\Biggr{[}\frac{\partial
n_{h}^{(2)}}{\partial\xi}+\frac{(\alpha-2)}{2}\frac{\partial
n_{h}^{(1)^{2}}}{\partial\xi}\Biggr{]}-\eta\frac{\partial^{2}u_{hz}^{(1)}}{\partial\xi^{2}}=0,$
(36) $\displaystyle\frac{\partial
n_{l}^{(1)}}{\partial\tau}-v_{p}\frac{\partial
n_{l}^{(2)}}{\partial\xi}+l_{x}\frac{\partial
u_{lx}^{(1)}}{\partial\xi}+l_{y}\frac{\partial
u_{ly}^{(1)}}{\partial\xi}+l_{z}\frac{\partial u_{lz}^{(2)}}{\partial\xi}$
$\displaystyle+l_{z}\frac{\partial}{\partial\xi}\Big{(}n_{l}^{(1)}u_{lz}^{(1)}\Big{)}=0,$
(37) $\displaystyle\frac{\partial
u_{lz}^{(1)}}{\partial\tau}-v_{p}\frac{\partial
u_{lz}^{(2)}}{\partial\xi}+l_{z}u_{lz}^{(1)}\frac{\partial
u_{lz}^{(1)}}{\partial\xi}+\mu_{2}l_{z}\frac{\partial\phi^{(2)}}{\partial\xi}$
$\displaystyle+\alpha\mu_{3}l_{z}K_{2}\Biggr{[}\frac{\partial
n_{l}^{(2)}}{\partial\xi}+\frac{(\alpha-2)}{2}\frac{\partial
n_{l}^{(1)^{2}}}{\partial\xi}\Biggr{]}-\eta\frac{\partial^{2}u_{lz}^{(1)}}{\partial\xi^{2}}=0,$
(38)
$\displaystyle\mu_{5}n_{l}^{(2)}+n_{h}^{(2)}=\beta_{1}\phi^{(2)}+\beta_{2}\phi^{(1)^{2}}.\
$ (39)

Finally, the next higher-order terms of Eqs. (8)-(11) and (15), with the help
of Eqs. (25)-(39), can provide the Burgers’ equation

$\displaystyle\frac{\partial\Phi}{\partial\tau}+A\Phi\frac{\partial\Phi}{\partial\xi}=B\frac{\partial^{2}\Phi}{\partial\xi^{2}},\
$ (40)

where $\Phi=\phi^{(1)}$ for simplicity. In Eq. (32), the nonlinear coefficient
$A$ and dissipative coefficient $B$ are given

$\displaystyle A=A^{\prime}(A^{\prime\prime}-2\beta_{2}),$ (41) $\displaystyle
B=\frac{\eta}{2},\ $ (42)

where

$\displaystyle
A^{\prime}=\frac{v_{p}^{4}+\alpha^{2}\mu_{1}\mu_{3}l_{z}^{4}K_{1}K_{2}-\alpha
l_{z}^{2}v_{p}^{2}\big{(}\mu_{1}K_{1}+\mu_{3}K_{2}\big{)}}{2v_{p}l_{z}^{2}\big{[}(v_{p}^{2}-\alpha\mu_{3}l_{z}^{2}K_{2})^{2}+\mu_{2}\mu_{5}(v_{p}^{2}-\alpha\mu_{1}l_{z}^{2}K_{1})^{2}\big{]}},$
$\displaystyle
A^{\prime\prime}=\frac{l_{z}^{4}\big{[}\\{3v_{p}^{2}+\mu_{1}l_{z}^{2}K_{1}\alpha(\alpha-2)\\}+\mu_{2}^{2}\mu_{5}\\{3v_{p}^{2}+\mu_{3}l_{z}^{2}K_{2}\alpha(\alpha-2)\\}\big{]}}{\big{\\{}(v_{p}^{2}-\alpha\mu_{1}l_{z}^{2}K_{1})(v_{p}^{2}-\alpha\mu_{3}l_{z}^{2}K_{2})\big{\\}}^{3}},\
$

Now, we look for stationary shock wave solution of this Burgers’ equation by
considering $\zeta=\xi-U_{0}\tau$ (where $\zeta$ is a new space variable and
$U_{0}$ is the speed of the ion fluid). These allow us to write the stationary
shock wave solution as

$\displaystyle\Phi=\Phi_{0}\bigg{[}1-\tanh\bigg{(}\frac{\zeta}{\Delta}\bigg{)}\bigg{]},\
$ (43)

where the amplitude $\Phi_{0}$ and width $\Delta$ are, respectively, given by

$\displaystyle\Phi_{0}=\frac{U_{0}}{A},~{}~{}~{}\mbox{and}~{}~{}~{}\Delta=\frac{2B}{U_{0}}.\
$ (44)

It is clear from Eqs. (43) and (44) that the IASHWs exist, which are formed
due to the balance between nonlinearity and dissipation, because $B>0$ and the
IASHWs with $\Phi>0$ ($\Phi<0$) exist if $A>0$ ($A<0$) because $U_{0}>0$.

Figure 1: Plot of the nonlinear coefficient $A$ vs $\mu_{5}$ when
$\alpha=5/3$, $\gamma_{e}=4/3$, $\gamma_{p}=4/3$, $\delta=20^{\circ}$,
$\eta=0.3$, $Z_{h}=37$, $Z_{l}=6$, $n_{h0}=7\times 10^{29}\mbox{cm}^{-3}$,
$n_{l0}=5\times 10^{30}\mbox{cm}^{-3}$, $n_{p0}=10^{31}\mbox{cm}^{-3}$, and
$v_{p}=v_{p+}$. Figure 2: Plot of $\Phi$ vs $\zeta$ for different values of
$\delta$ when $\alpha=5/3$; $\gamma_{e}=4/3$, $\gamma_{p}=4/3$, $\eta=0.3$,
$Z_{h}=37$, $Z_{l}=6$, $n_{h0}=7\times 10^{29}\mbox{cm}^{-3}$, $n_{l0}=5\times
10^{30}\mbox{cm}^{-3}$, $n_{p0}=10^{31}\mbox{cm}^{-3}$, $U_{0}=0.05$, and
$v_{p}=v_{p+}$. Figure 3: Plot of $\Phi$ vs $\zeta$ for different values of
$\eta$ when $\alpha=5/3$; $\gamma_{e}=4/3$, $\gamma_{p}=4/3$,
$\delta=20^{\circ}$, $Z_{h}=37$, $Z_{l}=6$, $n_{h0}=7\times
10^{29}\mbox{cm}^{-3}$, $n_{l0}=5\times 10^{30}\mbox{cm}^{-3}$,
$n_{p0}=10^{31}\mbox{cm}^{-3}$, $U_{0}=0.05$, and $v_{p}=v_{p+}$. Figure 4:
Plot of $\Phi$ vs $\zeta$ for different values of $Z_{l}$ when $\alpha=5/3$;
$\gamma_{e}=4/3$, $\gamma_{p}=4/3$, $\delta=20^{\circ}$, $\eta=0.3$,
$Z_{h}=37$, $n_{h0}=7\times 10^{29}\mbox{cm}^{-3}$, $n_{l0}=5\times
10^{30}\mbox{cm}^{-3}$, $n_{p0}=10^{31}\mbox{cm}^{-3}$, $U_{0}=0.05$, and
$v_{p}=v_{p+}$. Figure 5: Plot of $\Phi$ vs $\zeta$ for different values of
$Z_{h}$ when $\alpha=5/3$; $\gamma_{e}=4/3$, $\gamma_{p}=4/3$,
$\delta=20^{\circ}$, $\eta=0.3$, $Z_{l}=6$, $n_{h0}=7\times
10^{29}\mbox{cm}^{-3}$, $n_{l0}=5\times 10^{30}\mbox{cm}^{-3}$,
$n_{p0}=10^{31}\mbox{cm}^{-3}$, $U_{0}=0.05$, and $v_{p}=v_{p+}$. Figure 6:
Plot of $\Phi$ vs $\zeta$ for different values $n_{l0}$ when $\alpha=5/3$;
$\gamma_{e}=4/3$, $\gamma_{p}=4/3$, $\delta=20^{\circ}$, $\eta=0.3$,
$Z_{h}=37$, $Z_{l}=6$, $n_{h0}=7\times 10^{29}\mbox{cm}^{-3}$,
$n_{p0}=10^{31}\mbox{cm}^{-3}$, $U_{0}=0.05$, and $v_{p}=v_{p+}$. Figure 7:
Plot of $\Phi$ vs $\zeta$ for different values of $n_{h0}$ when $\alpha=5/3$;
$\gamma_{e}=4/3$, $\gamma_{p}=4/3$, $\delta=20^{\circ}$, $\eta=0.3$,
$Z_{h}=37$, $Z_{l}=6$, $n_{l0}=5\times 10^{30}\mbox{cm}^{-3}$,
$n_{p0}=10^{31}\mbox{cm}^{-3}$, $U_{0}=0.05$, and $v_{p}=v_{p+}$. Figure 8:
Plot of $\Phi$ vs $\zeta$ for different values of $n_{p0}$ when $\alpha=5/3$;
$\gamma_{e}=4/3$, $\gamma_{p}=4/3$, $\delta=20^{\circ}$, $\eta=0.3$,
$Z_{h}=37$, $Z_{l}=6$, $n_{h0}=7\times 10^{29}\mbox{cm}^{-3}$, $n_{l0}=5\times
10^{30}\mbox{cm}^{-3}$, $U_{0}=0.05$, and $v_{p}=v_{p+}$.

##  4 Results and Discussion

Our investigation is reasonably valid for degenerate cold plasma systems
(viz., white dwarfs and neutron stars) in which light (viz.,
${}^{1}_{1}{\mbox{H}}$ [1, 2, 3, 4], ${}^{4}_{2}{{\mbox{He}}}$ [1, 5, 6, 7,
8], ${}^{12}_{6}{\mbox{C}}$ [7, 8], and ${}^{16}_{8}{\mbox{O}}$ [7, 8, 9],
etc.) and heavy (viz., ${}^{56}_{26}{\mbox{Fe}}$ [10, 11],
${}^{85}_{37}{\mbox{Rb}}$ [11, 12], and ${}^{96}_{42}{\mbox{Mo}}$ [11, 12],
etc.) ions can exist. For our numerical analysis, we have considered the range
of the plasma parameters as $Z_{l}=6\thicksim 13$, $Z_{h}=35\thicksim 65$,
$n_{l0}=3\times 10^{30}\mbox{cm}^{-3}\thicksim 7\times 10^{30}\mbox{cm}^{-3}$,
$n_{h0}=4\times 10^{29}\mbox{cm}^{-3}\thicksim 10\times
10^{29}\mbox{cm}^{-3}$, and $n_{p0}=1\times 10^{31}\mbox{cm}^{-3}\thicksim
10\times 10^{31}\mbox{cm}^{-3}$.

It is obvious from Eq. (40) that the wave potential becomes infinite when
nonlinear coefficient is equal to zero (i.e., $A=0$), and in that case, the
validity of reductive perturbation method breaks down. The positive (i.e.,
$\Phi>0$) and negative (i.e., $\Phi<0$) electrostatic potentials can exist
corresponding to the values of $A>0$ and $A<0$, respectively. It is clear from
Fig. 1 that our plasma model supports only positive potential shock structure
(i.e., $\Phi>0$) associated with $A>0$ under consideration of non-relativistic
light and heavy ions (i.e., $\alpha=5/3$), and ultra-relativistic degenerate
electrons and positrons (i.e., $\gamma_{e}=\gamma_{p}=4/3$).

The oblique angle $\delta$ indicates the angle between the direction of the
propagation of IASHW and the direction of existing external magnetic field
which is parallel with $z$-axes. The variation of electrostatic shock
potential structure (i.e., $\Phi>0$) associated with $A>0$ with $\delta$ under
consideration of non-relativistic light and heavy ions (i.e., $\alpha=5/3$)
and ultra-relativistic degenerate electrons and positrons (i.e.,
$\gamma_{e}=\gamma_{p}=4/3$) can be observed in Fig. 2. It is obvious from
this figure that the electrostatic shock potential (i.e., $\Phi>0$) associated
with $A>0$ increases with the increase of $\delta$. Physically, the
interaction between the electrostatic shock potential associated with $A>0$
and the external magnetic field is clearly increased with the increase in
$\delta$.

Figure 3 describes the effects of the kinematic viscosity of the positively
charged non-relativistic light and heavy ions (via $\eta$) on the
electrostatic shock potential structure (i.e., $\Phi>0$) associated with
$A>0$. We have noticed that the amplitude of the electrostatic shock profile
(i.e., $\Phi>0$) associated with $A>0$ is independent to the variation of the
ion kinematic viscosity but the steepness of the electrostatic shock profile
(i.e., $\Phi>0$) associated with $A>0$ is rigorously dependent to the
variation of ion kinematic viscosity. The steepness of the electrostatic shock
potential structure (i.e., $\Phi>0$) associated with $A>0$ decreases with
$\eta$, and this result agrees with the result of Refs. [31, 32].

The amplitude of the positive electrostatic shock structure (i.e., $\Phi>0$)
associated with $A>0$ is so much sensitive to the change of the charge state
of non-relativistic light and heavy ions. Figures 4 and 5 show that the
amplitude of the positive shock structure (i.e., $\Phi>0$) associated with
$A>0$ increases with the charge state of non-relativistic light and heavy
ions. Similarly, the increasing number of non-relativistic light and heavy
ions enhances the positive electrostatic shock structure (i.e., $\Phi>0$)
associated with $A>0$. It can easily demonstrate from Figs. 6 and 7 that as we
increase $n_{l0}$ and $n_{h0}$, the amplitude of the positive electrostatic
shock structure (i.e., $\Phi>0$) associated with $A>0$ increases. Physically,
the charge state and number density of non-relativistic light and heavy ions
can control the dynamics of the DQP system in the similar way.

We have studied the characteristics of IASHWs for different values of positron
number density (via $n_{p0}$) under consideration of non-relativistic light
and heavy ions (i.e., $\alpha=5/3$) and ultra-relativistic degenerate
electrons and positrons (i.e., $\gamma_{e}=\gamma_{p}=4/3$) in Fig. 8. It can
be highlighted from this figure that the amplitude of electrostatic shock
structure (i.e., $\Phi>0$) associated with $A>0$ decreases with positron
number density, and this finding is analogous to the result of Refs. [21, 22]
.

##  5 Conclusion

We have studied the basic characteristics of IASHWs in an extremely dense DQP
containing non-relativistic light and heavy ions, and inertialess ultra-
relativistic degenerate electrons and positrons in the presence of external
magnetic field. The RPM [41] has been utilized to derive Burgers’ equation.
The results that have been found from our present investigation can be
summarized as follows:

  * 1.

Our plasma model supports only positive potential shock structure (i.e.,
$\Phi>0$) associated with $A>0$ under consideration of non-relativistic light
and heavy ions (i.e., $\alpha=5/3$), and ultra-relativistic degenerate
electrons and positrons (i.e., $\gamma_{e}=\gamma_{p}=4/3$).

  * 2.

The electrostatic shock potential (i.e., $\Phi>0$) associated with $A>0$
increases with the increase in $\delta$.

  * 3.

The steepness of the positive potential shock structure (i.e., $\Phi>0$)
associated with $A>0$ decreases with $\eta$.

  * 4.

The amplitude of electrostatic shock structure (i.e., $\Phi>0$) associated
with $A>0$ is to be found to increase with the charge state and number density
of non-relativistic light and heavy ions.

  * 5.

The increasing positron number density decreases the height of the positive
shock profile.

It may be noted that the self-gravitational effect of the plasma species is
important to be considered in our governing equations but beyond the scope of
our present work. Overall, the outcomes from our present investigation will be
helpful to understand the IASHWs in white dwarfs and neutron stars.

## References

  * [1] S. Chandrasekhar, Astrophys. J. 74, 81 (1931).
  * [2] S. Chandrasekhar, The Observatory 57, 373 (1934).
  * [3] R.S. Fletcher, et al., Phys. Rev. Lett. 96, 105003 (2006).
  * [4] T.C. Killian, Nature (London) 441, 297 (2006).
  * [5] H.M. Van Horn, Science 252, 384 (1991).
  * [6] R.H. Fowler, J. Astrophys. Astron. 15, 105 (1994).
  * [7] D. Koester and G. Chanmugam, Rep. Prog. Phys. 53, 837 (1990).
  * [8] D. Koester, Astron. Astrophys. Rev. 11, 33 (2002).
  * [9] S. Chandrasekhar, Philos. Mag. 11, 592 (1931).
  * [10] A. Vanderburg, et al., Nature 526, 546 (2015).
  * [11] N. Chamel and P. Haensel, Living Rev. Relativ. 11, 10 (2008).
  * [12] A. Witze, Nature 510, 196 (1990).
  * [13] S. Chandrasekhar and R.F. Tooper, Astrophys. J. 139, 1396 (1964).
  * [14] P.K. Shukla, et al., Phys. Rev. E 84, 026405 (2011).
  * [15] W.F. El-Taibany, et al., Adv. space res. 50, 101 (2012).
  * [16] S. Sultana and R. Schlickeiser, Astrophys. Space Sci. 363, 103 (2018).
  * [17] B. Zhang and J. Gil, Astrophys. J. 631, 143 (2005).
  * [18] P.A. Sturrock, Astrophys. J. 164 (1971) 529.
  * [19] A.K. Harding and A. G. Muslimov, Astrophys. J. 508, 328 (1998).
  * [20] A.K. Harding and D. Lai, Rep. Prog. Phys. 69, 2631 (2006).
  * [21] T.S. Gill, et al., J. Phys. 208, 012040 (2010).
  * [22] Ata-Ur-Rahman, et al., J. Plasma Phys. 79, 817 (2013).
  * [23] M.R. Hossen, et al., Plasma Phys. Rep. 43, 1189 (2017).
  * [24] A.A. Mamun and P.K. Shukla, Phys. Plasmas 17, 104504 (2010).
  * [25] P.M.S. Blackett, Nature 159, 658 (1947).
  * [26] J. Liebert, et al., Astrophys. J. 214, 457 (1977).
  * [27] F. Euchner, et al., Astron. Astrophys. 390, 633 (2002).
  * [28] W.F. El-Taibany, et al., Adv. Space Res. 50, 101 (2012).
  * [29] M.I. Shaukat, Phys. Plasmas 24, 102301 (2017).
  * [30] J.M. Burgers, Elsevier 1, 171 (1948).
  * [31] M.G. Hafez, et al., Plasma Phys. Rep. 43, 499 (2017).
  * [32] H.G. Abdelwahed, et al., J. Exp. Theor. Phys. 122, 1111 (2016).
  * [33] N.S. Saini, et al., Wave Random Complex. DOI: 10.1080/17455030.2020.1798561.
  * [34] M. M. Haider, Z. Naturforsch. A 71, 1131 (2016).
  * [35] A. Atteya, et al., Eur. Phys. J. Plus 132, 1 (2017).
  * [36] M.H. Rahman,et al., Phys. Plasmas 25, 102118 (2018); N.A. Chowdhury, et al., Phys. plasmas 24, 113701 (2017); M.H. Rahman, et al., Chin. J. Phys. 56, 2061 (2018); N.A. Chowdhury, et al., Plasma Phys. Rep. 45, 459 (2019); S. Jahan, et al., Plasma Phys. Rep. 46, 90 (2020).
  * [37] S.K. Paul, et al., Pramana J. Phys. 94, 58 (2020); N.A. Chowdhury, et al., Chaos 27, 093105 (2017); N.A. Chowdhury, et al., Contrib. Plasma Phys. 58, 870 (2018); M. Hassan, et al., Commun. Theor. Phys. 71, 1017 (2019); D.M.S. Zaman, et al., High Temp. 58, 789 (2020).
  * [38] A. A. Mamun, Phys. Scr. 59, 454 (1999).
  * [39] M. M. Hossen, et al., High Energy Density Phys. 24, 9 (2017).
  * [40] H. Washimi, T. Taniuti, Phys. Rev. Lett. 17, 996 (1966).
  * [41] S. Jahan, et al., Commun. Theor. Phys. 71, 327 (2019); N. Ahmed, et al., Chaos 28, 123107 (2018); N.A. Chowdhury, et al., Vacuum 147, 31 (2018); R.K. Shikha, et al., Eur. Phys. J. D 73, 177 (2019); S. Banik, et al., Eur. Phys. J. D 75, 43 (2021).