# Effect of non-equilibrium thermochemistry
on Pitot pressure measurements in shock tunnels
(or: Is 0.92 really the magic number?)

Tamara Sopek 111Postdoctoral Research Fellow, School of Engineering.
University of Southern Queensland, Toowoomba, Queensland, 4350, Australia
Peter Jacobs 222Assoc. Professor, Institute for Advanced Engineering and Space
Sciences University of Southern Queensland, Toowoomba, Queensland, 4350,
Australia  Suria-Devi Subiah333DPhil student, Department of Engineering
Science University of Oxford, Oxford, OX1 2JD, UK  Peter Collen444Postdoctoral
Research Assistant University of Oxford, Oxford, OX1 2JD, UK  and Matthew
McGilvray555Professor, Department of Engineering Science. University of
Oxford, Oxford, OX1 2JD, UK

###### Abstract

Pitot pressure is the most common measurement in high total enthalpy shock
tunnels for test condition verification. Nozzle calculations using multi-
temperature non-equilibrium thermochemistry are needed in conjunction with
Pitot measurements to quantify freestream properties. Pitot pressure is
typically matched by tuning the boundary layer transition location in these
simulations. However, non-equilibrium thermochemistry effects on the Pitot
probe are commonly ignored. A computational study was undertaken to estimate
the effect of non-equilibrium thermochemistry on Pitot pressure and freestream
conditions. The test flow was produced by a Mach 7 nozzle in a reflected shock
tunnel for air at a relatively low total enthalpy of 2.67 MJ/kg. Three
different thermochemical models (equilibrium, finite-rate chemistry and two-
temperature thermochemistry) were employed to compute flow variables at the
nozzle exit and Pitot probe. Pitot pressures from these simulations were
compared against those obtained via experiments. The results show a departure
from the commonly utilized C of 0.92 in the reduced Rayleigh-Pitot equation
form for high Mach numbers. Additionally, calculations were done with a sweep
of free-stream conditions and resulting in C$\approx$0.94. These results show
that the influence of finite-rate thermochemistry should be taken into
account, even at relatively low flow enthalpies.

## Nomenclature

$C_{p,\,max}$ = | maximum value of the pressure coefficient
---|---
$C$ = | constant, $C$=0.5*$C_{p,\,max}$
$M$ = | Mach number
$p$ = | pressure, Pa
$T$ = | temperature, K
$U$ = | axial velocity, m/s
$\rho$ = | density
Subscripts
$tr$ = | translation and rotational
$ve$ = | vibrational and electronic
Abbreviations
ESTCN = | Equilibrium Shock Tube Conditions, with Nozzle
RST = | reflected shock tunnel
NENZF1d = | Non-Equilibrium Nozzle Flow, one dimension

##  1 Introduction

Shock tunnels are a class of impulse facilities used to investigate hypersonic
flows [1]. They present a fundamental instrument to replicate conditions
encountered during high-speed flight. Because these facilities can simulate a
range of temperatures and pressures and use different test gas compositions,
they are a useful device for the investigation of conditions that might be
difficult to achieve in other facilities. However, the flow state and
variables must be determined to a sufficiently high level of accuracy to allow
analysis of the obtained data with reasonably low uncertainty. Typically,
Pitot pressure measurements are performed to verify the condition and size of
the core flow at the nozzle exit. To estimate additional properties of the
test flow, a computational study is usually undertaken from stagnation
conditions as it is unpractical to measure all the relevant flow variables.
However, the computed variables are highly dependent on the type of the gas
model and boundary layer modeling. In hypersonic high temperature gases, the
thermochemistry changes substantially with the increase of flow speed. These
gas effects, such as excitation of vibrational energy modes, may affect the
assumed Pitot pressure, though this is typically ignored.

The Rayleigh-Pitot formula [2] shown in (Eq. 1) gives a relationship between
stagnation and static pressure as a function of Mach number and ratio of
specific heats. It is derived as the product of the normal shock and
isentropic relations for pressure at stagnation point behind the normal shock
wave, and, for an ideal gas can be written as:

$\frac{\textit{p}_{0,2}}{\textit{p}_{1}}=\left[\frac{(1+\gamma)M_{1}^{2}}{2}\right]^{\frac{\gamma}{\gamma-1}}\left[\frac{\gamma+1}{2\gamma
M^{2}_{1}-(\gamma-1)}\right]^{\frac{1}{\gamma-1}}$ (1)

which, in the limit of high Mach number can be reduced to:

$\textit{p}_{0,2}=\textit{p}_{Pitot}=\frac{1}{2}\rho\textit{U}^{2}\cdot
C_{p,\,max}=C\rho\textit{U}^{2}$ (2)

where $C_{p\,max}$ and $C=\frac{1}{2}C_{p,\,max}$ are constants. $C_{p,\,max}$
is the maximum value of the pressure coefficient, evaluated at the stagnation
point behind a normal shock wave. Full derivation of the maximum pressure
coefficient, $C_{p,\,max}$, is presented in the Appendix, in Eqs 5, 6, 7,
LABEL:Cpmax-infinite-Mach-2, 9 and 10.

It is shown in Eq.10 from the Appendix that in the limit of high hypersonic
flows $C_{p,\,max}$ reaches the asymptotic value $C_{p,\,max}=1.8394$ for
$\gamma=1.4$ [3]. C which is used to compute $\textit{p}_{Pitot}$ is then
obtained by halving that value, i.e. C$\,=0.9197$ or more commonly a rounded
value of 0.92 [4]. This value of 0.92 is for an ideal gas model, which implies
that this value should be validated when the effect of thermochemistry present
in real shock tunnel nozzles is in play. The practice of using the reduced
form Rayleigh-Pitot equation goes back decades as seen in the tables starting
at page 25 in Ref.[4].

Kantrowitz [5] conducted a study investigating the influence of non-
equilibrium thermochemistry on energy dissipation. He performed measurements
with CO2 where the gas was expanded from a chamber through a nozzle, and
compressed at the nose of a Pitot tube aligned with the jet from the nozzle.
Pitot pressure was measured over a range of chamber pressures. The results
show dependence of the resultant Pitot pressure on the relaxation time of the
gas when the gas is in non-equilibrium, the measured Pitot pressure will be
lower than the chamber pressure, and this “loss” of Pitot pressure increases
with increased chamber pressure. These initial findings were confirmed by a
subsequent study [6] conducted with a range of gases.

Work was done previously [7] to investigate the effects of hypersonic and high
temperature gas on Pitot pressure. A parametric study was completed by
altering the flow conditions to isolate the effects of individual flow
characteristics, for example by running simulations with thermal effects on,
followed by thermal effects off. Calculations were performed for a range of
conditions and two nozzles (M6 and M8) using both an inviscid MATLAB code and
US3D code [8]. The full-form for an ideal gas (Eq. 1) and reduced form (Eq. 2)
of Rayleigh-Pitot equation were used to assess their accuracy in predicting
Pitot pressure without numerical modeling of Pitot probe. In that study it was
suggested that a value of 0.94 be used for C in the reduced form of the
Rayleigh-Pitot equation.

The aim of this work was to investigate the effect of non-equilibrium
thermochemistry on Pitot pressure measurements in shock tunnels, as well as
the validity of the Rayleigh-Pitot equation in both general and reduced form
used in hypersonic flows for different thermochemical models. Numerical
calculations were performed and compared against experimental data for a low
enthalpy air condition in the Oxford T6 Stalker facility operated in RST mode.
The nozzle-exit flow values were then subsequently used as an inflow for Pitot
probe simulations. The reduced form Rayleigh-Pitot equation
($\textit{p}_{Pitot}$=C$\rho\textit{U}^{2}$) for hypersonic flows was
calculated from the results of the simulations performed using three
thermochemical models.

##  2 Numerical methodology

Numerical results were obtained with two types of calculations: the first
calculation was performed for just the nozzle, and then the second calculation
was done using outflow from the nozzle simulation as an inflow for the Pitot
probe simulation. Positioning of the nozzle and the Pitot probe was done to
emulate that in measurements, as shown in Fig. 2.

Figure 1: Placement of the nozzle and Pitot probe.

###  2.1 Eilmer Code

The computational study was undertaken using the Eilmer code [9], which is an
open-source 666Eilmer code is a part of the Gas Dynamics Toolkit and freely
available at https://gdtk.uqcloud.net Reynolds-averaged Navier-Stokes
computational fluid dynamics solver. This CFD solver is a collection of codes
for the numerical simulation of transient (time-accurate), compressible gas
flows in two and three dimensions. It has been validated on numerous high-
speed (supersonic and hypersonic) test cases, and it can be used to model
these highly energetic flows in chemical and thermal equilibrium and non-
equilibrium. This ability makes it ideal for high total enthalpy flows where
the flow behind the shock wave experiences very high temperatures (sometimes
in excess of 10000 K) where non-equilibrium thermochemistry should be
considered.

###  2.2 Thermochemistry modeling

The simulations were performed with three different thermochemical gas models
to investigate the influence of high-enthalpy chemical kinetics. These models
include equilibrium thermochemistry (based on 5-species air model in the CEA2
program [10]), finite-rate chemistry (based on Gupta’s 5-species air model
[11]) and two-temperature thermochemistry (based on Park’s 2-temperature
5-species air model [11]). These thermochemistry models are referred to,
respectively, in this section of the paper as the equilibrium, chemical non-
equilibrium and thermochemical non-equilibrium models. The mixture transport
properties, that is, mixture viscosity and thermal conductivity, were
calculated using mixing rules of Gordon and McBride [10] which are a variant
of Wilke’s original formulation [12].

###  2.3 Nozzle setup

(a)

(b)

Figure 2: Computational grid for a) Mach 7 nozzle and b) Pitot probe
simulations. Inflow is from the left-hand boundary.

For the axisymmetric nozzle simulations, the approach taken was identical to
that used in previous studies of similar facilities, e.g. [13], where the
nozzle is simulated from the throat onwards. It has been shown that the
differences when including the region upstream of the nozzle vs not including
it are very small, while the cost is significantly reduced when simulating
from the throat onwards. The stagnation properties (Table 2) are
isentropically expanded to a Mach number of 1 to give a uniform inflow at the
throat and that process is assumed to be in thermochemical equilibrium due to
subsonic and relatively slow transient to reaction times. The wall of the
nozzle was assumed to be at a constant temperature of 300 K. Shock tunnels
have short flow durations and the nozzle wall does not change in temperature
by more than a few degrees during the test flow time. A pseudo-space-marching
approach was used to reduce the computational time needed to obtain a steady-
state nozzle flow solution [14]. The turbulence model used in the simulations
was the $k-\omega$ turbulence model of Wilcox [15]. To further reduce
computational time, the wall functions of Nichols [16] were used to model the
turbulent boundary layers. The application of wall functions is appropriate
since boundary layer separation is not expected in the nozzle flow field. The
turbulence intensity of the inflow to the nozzle was set to 5%, with a ratio
of the turbulent-to-laminar viscosity of 100. To examine the influence of
boundary layer transition location on the core flow parameters, simulations
with different transition locations downstream of the nozzle throat were
performed. The boundary layer transition location was varied until simulations
with three different thermochemical models reached agreement with the
experimental Pitot pressure data that was within experimental uncertainty of
approximately 3%. The final transition location is 380 mm downstream of the
nozzle throat. More information about the experimental facility and the
process of obtaining Pitot pressure data can be found in [17].

Table 2:  Stagnation conditions for the nozzle simulation.  Stagnation property | Value | Unit
---|---|---
pst | 2.5469e+07 | Pa
Tst | 2389.63 | K
mass fraction [N2] | 0.76222 | -
mass fraction [O2] | 0.22734 | -
mass fraction [NO] | 0.01041 | -
mass fraction [O] | 3.1125e-05 | -
mass fraction [N] | 0.0 | -

(a)

(b)

Figure 3: Grid convergence analysis based on pressure for a) nozzle and b)
Pitot probe simulations.

(a)

(b)

(c)

Figure 4: Comparison of the equilibrium (top half) and the thermochemical non-
equilibrium (bottom half) flow over Pitot probe. Flow variables are a)
pressure, b) translational-rotational temperature, and c) density.

The computational grid which produced converged solutions for the flow
simulations had 600 cells in the axial direction and 80 cells in the radial
direction as shown in Fig.2a. Grid clustering was employed near the throat and
initial expansion regions to resolve the larger flow gradients expected in
these regions and near the nozzle wall to ensure adequate resolution of the
boundary layer. The grids near the nozzle wall were configured to be
orthogonal to the wall to ensure that the near-wall flow parameters were
properly approximated by the flow solver. The level of grid convergence is
shown in Fig.3a with a plot of mass-flow-weighted static pressure vs. the
inverse of total number of cells for each grid resolution. The mass-flow-
weighted static pressure value was obtained at the nozzle exit plane, with
averaging area corresponding to the fraction of nozzle radius that would be
equal to the radius of Pitot probe. The coarse grid had 300 axial cells and 40
radial cells, the medium grid had 600 axial cells and 80 radial cells, whilst
the fine grid had 900 axial cells and 120 radial cells. The mass-flow-weighted
pressure for the medium grid differed from the fine grid by only 0.7%. This
indicated that the medium grid produced a flow field at the exit of the nozzle
that was sufficiently grid-converged for the analysis conducted in the present
study at a cost of 1.56 hrs compute per flow solution. The output from the
nozzle simulation was further used as a freestream inflow for a separate
simulation of the Pitot probe, described in Sec 2.4. The position of the Pitot
probe relative to the nozzle is shown in Fig. 1.

###  2.4 Pitot probe

The Pitot probe was modeled using a simplified geometry, as illustrated in
Fig. 2b. This did not account for the internal geometry which can lead to flow
measurement fluctuations [18]. However, the internal geometry of the Pitot
probe is so small that it can be approximated that the sensor is near the
surface of the probe. The probe was placed at an axial distance of 0.0033 m
from the nozzle exit plane, approximating the position that would be used in
the measurements. The nozzle and the probe have matching axes of symmetry
(i.e. $y$ = 0). The Pitot probe near the stagnation point was modeled as
inviscid flow because the flow is stagnated. There is no radius on the
physical Pitot probe tip, or at worst, the radius is negligibly small, thus
the probe was modeled with sharp edges. Grid clustering was employed around
the wall. The grid independence was checked by plotting mass-flow weighted
static pressure at the probe face against the inverse of the total number of
cells for each grid resolution, as shown in Fig. 3b. The mass-flow-weighted
static pressure value was obtained at the probe face, with averaging area
corresponding to the fraction of the probe face radius that would experience
the flow. The coarse grid had the total of 37500, the medium grid 150000, and
the fine grid 337500 cells. It can be seen in Fig. 3b that the grid
convergence was achieved as the difference in pressure levels for all three
grid sizes is <1%. The results reported here were obtained using the medium
grid. The area with 1.5 mm radius from the axis of symmetry is used for Pitot
pressure averaging as this area is representative of the Pitot probe opening.

Figure 4 shows a comparison of the flowfield for the equilibrium (top half of
the image) and thermochemical non-equilibrium (bottom half of the image) flow
over the Pitot probe, to illustrate the effect thermochemical non-equilibrium
has on the freestream properties. Only equilibrium and thermochemical non-
equilibrium are shown here as the differences between these two flowfields are
the largest. Variables shown are pressure, temperature and density, while it
is also possible to observe the difference in the bow shock stand-off distance
for different thermochemistry. The largest difference in the flowfield values
is observed in the region between the Pitot probe face and the bow shock. In
this region, pressure plot demonstrates a slightly lower pressure for the
thermochemical non-equilibrium compared to the equilibrium flow. Temperature
is higher for the case of thermochemical non-equilibrium, while density is
lower. As these flow properties give most of the information about the
flowfield, it is clear from Fig. 4 that it is necessary to consider the effect
of thermochemical non-equilibrium on the flow.

##  3 Pitot pressure predictions

A comparison of experimental versus simulated Pitot pressure at the
measurement location is provided in Fig. 5. Numerical data are computed using
the reduced form of the Rayleigh-Pitot equation and flow variables from the
simulations with values taken at the centerline of the nozzle exit plane. The
experimental values shown are measured Pitot pressure averaged over the test
time. The results show good uniformity across the core flow region, which the
experimental data suggests has a diameter of at least 160 mm. In Fig. 5a, the
numerical $\textit{p}_{Pitot}$ values are computed via the reduced form of the
Rayleigh-Pitot formula and using the C values obtained from the respective
simulations at the centerline of the nozzle exit plane. Fig. 5b shows
numerical Pitot pressures obtained using Eq. 2 with a C value of 0.92 for all
three thermochemical models to illustrate the difference when this value is
applied to all cases, compared to Fig. 5a where appropriately tuned C is
applied. Fig. 5b was previously reported in Ref.[17], an overview of the T6
Stalker facility. Table 3 provides C values for all three thermochemical
models obtained from Pitot probe simulations, where tabulated C values were
obtained from the post-shock stagnation point conditions of the Pitot probe
simulation using a respective thermochemical model. These computed C values
show a dependence on the non-equilibrium relaxation effects such that C
values, and thus computed Pitot pressures, decrease with the increase in
relaxation times, a result that is consistent with the findings of Ref. [5,
6].

(a) Numerical Pitot pressure obtained with coefficient C obtained from
respective simulations using Rayleigh-Pitot equation and computed flow
variables: CEQ=0.94388, C${}_{\bf{CHEM.~{}NON-EQ}}$=0.94137,
C${}_{\bf{THERMOCHEM.~{}NON-EQ}}$=0.93002

(b) Numerical Pitot pressure obtained with coefficient C equal to 0.92 for all
cases.

Figure 5: Pitot pressure comparison between experimental and numerical data.
Numerical results for all three cases were obtained for a boundary layer
transition location set at 0.38 m downstream of the nozzle throat.

Fig. 5 demonstrates that overall, there is an excellent match between the
experiments and the simulations. However, Fig. 5a illustrates that within the
core flow region, the thermochemical non-equilibrium result shows slightly
better agreement with measurement data points. Additionally, it is shown that
by using the tuned coefficient values to compute Pitot pressure, the resulting
values for Pitot pressures for all three thermochemical models agree somewhat
better with measurement data points than when using C = 0.92 for all cases.
While these improvements in agreement might be considered small judging from
the total values only, they become more significant considering that the
experimental errors (as reported on the measurement data points) are already
quite low. It will be shown in the following sections what is the effect on
the total test flow properties.

Most experimenters calculate Pitot pressures using the Rayleigh-Pitot relation
(Eq. 1) or its reduced form (Eq. 2) for ideal air at high Mach number (with
the assumption of a constant ratio of specific heats of 1.4). However,
applying the Rayleigh-Pitot relation (Eq. 1) in the thermochemical non-
equilibrium case resulted in Pitot pressure values that were not in agreement
with the other two thermochemistry models, which is not surprising considering
that the Rayleigh-Pitot equation was derived for an ideal gas (both thermally
and calorically perfect).

Table 3:  Constant in Rayleigh-Pitot equation for high Mach numbers, obtained for different thermochemical calculations for both nozzle and Pitot probe.  Thermochemical model | Constant C
---|---
equilibrium | 0.94388
chemical non-equilibrium | 0.94137
thermochemical non-equilibrium | 0.93002

Table 4 presents Pitot pressure values obtained from simulations with
different combinations of thermochemistry models in the nozzle simulations and
the subsequent Pitot probe simulations, i.e., the thermochemistry model in
nozzle simulation for some cases differed from the thermochemistry model of
the Pitot probe simulation. This resulted in a maximum of 5% difference on the
Pitot pressures (as shown in the matching Pitot pressure levels in Fig. 5).
This is not surprising since the Pitot pressure value is driven by the
momentum of the gas which is not greatly influenced by the detailed
distribution of energy within the gas stream.

Table 4:  Pitot pressure averaged over area of the Pitot probe opening of $\phi$=3 mm. Expressed in percentage of difference off the nominal thermochemical non-equilibrium value of 337 kPa. | 0.92$\rho U^{2}$ | EQ | CHEM. NON-EQ | THERMOCHEM. NON-EQ
---|---|---|---|---
nozzle sim. | EQuilibrium | 2.28% | 4.96% | 4.12% | 3.29%
CHEM. NON-EQuilibrium | 1.09% | 3.38% | 3.46% | 2.81%
THERMOCHEM. NON-EQuilibrium | -1.09% | 0.66% | 0.87% | 0.00%

##  4 Effect on freestream conditions

While different thermochemistry appears to have little effect on Pitot
pressure, it is expected to have a more significant effect on other flow
properties. To explain this, simulations are undertaken for the thermochemical
non-equilibrium case. In this work, we are assuming the correct value of the
Pitot pressure is the one obtained through the same chemistry assumption for
both nozzle and Pitot probe. The correct method is to tune to the Pitot probe
measurement to obtain correct freestream values. This is done by changing the
boundary layer transition location in the nozzle. Again, we are following the
typical approach of an experimenter where reduced form of the Rayleigh-Pitot
equation with 0.92 is used to compute numerical Pitot pressure that we are
trying to match to the measured Pitot pressure.

The results of this process are listed in Table 5, where the baseline case is
the thermochemical non-equilibrium nozzle and Pitot probe, resulting in Pitot
pressure of 337 kPa. The percentage differences in columns 2 and 3 in Table 5
were calculated using the following equation:

$\Delta V_{ESTCN}[\%]=\frac{V_{ESTCN}-V_{NON-EQ}}{V_{NON-EQ}}\cdot 100$ (3)

where $V$ is a freestream property, such as density, velocity, static
pressure, etc. Subscript ESTCN stands for calculations performed using ESTCN
with either non-equilibrium or equilibrium gas model, while subscript NON-EQ
relates to the result of the non-equilibrium nozzle and Pitot probe
simulations (i.e., the baseline case).

The difference in freestream values (column 3) is not massive, however, the
tuned transition location is only a few millimetres away from the initial one,
and so the values would not be significantly different between the initial
(column 2) and the transition-location-tuned simulation (column 3). The reason
for this is that the transition location was already tuned to match
experimental Pitot pressure, as explained in Sec. 2.4 and shown in Fig. 5.
Thus, the process explained here was the fine tuning of the transition
location and Pitot pressure. The resulting transition location is only 3 mm
upstream from the previous location of 380 mm, i.e., the resulting transition
location is now 377 mm. If less effort was put initially in tuning the
transition location, then naturally the process in this section would be
longer and the obtained differences in freestream properties would be larger.

Table 5:  Discrepancies in freestream values for the non-equilibrium and equilibrium nozzle.  Property | NON-EQ | $\Delta$NON-EQ nozzle | $\Delta$EQ nozzle
---|---|---|---
($p_{Pitot}$=0.92$\rho U^{2}$) | ($p_{Pitot}$=0.92$\rho U^{2}$)
$\rho$ [kg/m3] | 0.074 | 0.8% | -2.2%
$U$ [m/s] | 2214 | -0.02% | 1.6%
$p$ [kPa] | 5390 | 1.1% | 11.9%
$T_{tr}$ [K] | 253 | 0.3% | 14.4%
$T_{ve}$ [K] | 1271 | 0.01% | —
$M$ [-] | 7.1 | -0.2% | -7.7%

Table 5 also includes equilibrium nozzle calculations. As the experimenters
typically assume equilibrium flow properties at low temperatures, the
freestream conditions are generally assessed assuming the measured Pitot
pressure and tuning an equilibrium isentropic nozzle calculation to match that
value. Experimenters in reflected-shock tunnels usually perform this
equilibrium nozzle calculation using a state-to-state calculation program,
such as ESTCN [19], for estimating flow conditions. ESTCN simulations were run
targeting the Pitot pressure of the baseline case (non-equilibrium
thermochemistry simulations). Pitot pressure was computed using the reduced
form of the Rayleigh-Pitot equation with constant C equal to 0.92. It is
evident from the values in Table 5 that using the equilibrium calculations
with C = 0.92 to obtain Pitot pressure results in higher freestream
discrepancies compared to the baseline case than when the non-equilibrium
nozzle simulations with tuning the transition location are performed.

A sweep of free-stream conditions was performed, varying only enthalpy and
using the same Mach 7 nozzle. These calculations were done for air using the
following approach: first, one-dimensional nozzle simulations were done using
NENZF1d 777NENZF1d code is a part of the Gas Dynamics Toolkit and freely
available at https://gdtk.uqcloud.net. Second, using the nozzle simulation
result as an inflow, axisymmetric simulations were performed for a sphere, and
third, C values were obtained directly from the observed simulated pressure in
the cells against the stagnation point on the sphere as in the Eq. 4:

$C=\frac{\textit{p}_{simulation}}{\rho\textit{U}^{2}},$ (4)

where $\rho$ and U are freestream conditions used in the simulation.

These simulations were done for both one-temperature thermal model with
finite-rate chemistry (based on Gupta’s 5-species air model [11]) and two-
temperature thermal-non-equilibrium model with finite-rate chemistry (based on
Park’s 2-temperature 5-species air model [11]). The axisymmetric computational
grid was formed with 2400 cells.

Figure 6: Variation of C coefficient for a range of test conditions.

The results are shown in Fig. 6 and illustrate the variation of coefficient C
over a range of test conditions (i.e., enthalpies) and for two different
thermochemical models. It is demonstrated that the variation of the
coefficient C is most significant at lower enthalpies up to approximately 2.6
MJ/kg for the two-temperature thermal-non-equilibrium model, after which it
asymptotes at a value C $\approx$ 0.94. The data for one-temperature thermal
model is approximately 0.94. While there seems to be no reliable trend for C
vs. enthalpy, the approximate level of about 0.94 presents a useful starting
data point for researchers which can then be refined for any particular
operating point using a single blunt-body calculation. The commonly used value
C = 0.92 is not correct for all conditions, particularly at enthalpies higher
than approximately 2.6 MJ/kg, which makes up for the significant amount of
tests done in shock tunnels.

##  5 Conclusion

In this study, a range of CFD simulations were conducted to simulate flow in a
hypersonic nozzle. The nozzle-exit flow properties from these simulations
displayed departure from the equilibrium, highlighting the need to take the
influence of finite-rate thermochemistry into account. These results show that
the effects of thermochemical non-equilibrium are much stronger than expected
at a relatively low enthalpy condition of 2.67 MJ/kg. The values of constant C
in the reduced formulation of the Rayleigh-Pitot equation computed from our
simulations show departure from the commonly used value of 0.92. These
computed C values show a dependence on the non-equilibrium relaxation effects
such that C values, and thus computed Pitot pressures, decrease with the
increase in relaxation times, a result that is consistent with the findings of
Ref. [5, 6].

Our results show that the calculation of Pitot pressures using 0.92 value for
C lead to a maximum of 5% difference, a result that is not surprising since
values of Pitot pressure are driven by the momentum of the gas which is not
greatly influenced by the detailed distribution of energy within the gas
stream. However, our results also show that equilibrium calculations using C =
0.92 lead to much higher discrepancies in computed freestream properties
compared to the case than when the non-equilibrium nozzle simulations with
tuning the transition location are performed. The most significant effect is
on the static pressure and temperature with 11.9% and 14.4% difference,
respectively. This aligns with results from other reflected shock tunnel
facilities at similar enthalpies, for example [20]. Future experiments are
planned to confirm the degree of thermal non-equilibrium in the flow field,
both through optical techniques and measurement of the freestream static
pressure.

Additional calculations were done with a sweep of free-stream conditions,
varying only enthalpy. While the results did not show a reliable trend for C
vs. enthalpy, the approximate level of about 0.94 emerged from the data. This
value of C presents a useful starting data point for researchers and can be
further refined for any particular operating point using a single blunt-body
calculation. We recommend that suitable values for C should be estimated for
each operation point of a shock tunnel facility.

## Appendix

Across a normal shock wave, such as a bow shock at the nose of an aircraft,
maximum pressure coefficient is:

$C_{p,\,max}=\frac{\textit{p}_{0,2}-\textit{p}_{\infty}}{\frac{1}{2}\rho_{\infty}\textit{U}_{\infty}^{2}}$
(5)

Rewriting this in terms of freestream Mach number for supersonic conditions:

$C_{p,\,max}=\frac{\textit{p}_{0,2}-\textit{p}_{\infty}}{\frac{1}{2}\rho_{\infty}\textit{U}_{\infty}^{2}}=\frac{\textit{p}_{0,2}-\textit{p}_{\infty}}{\frac{\gamma}{2}\textit{p}_{\infty}M_{\infty}^{2}}=\frac{1}{{\frac{\gamma}{2}}M_{\infty}^{2}}\left[\frac{\left(\frac{\gamma+1}{2}M_{\infty}^{2}\right)^{\left(\frac{\gamma}{\gamma-1}\right)}}{\left(\frac{2\gamma}{\gamma+1}M_{\infty}^{2}-\frac{\gamma-1}{\gamma+1}\right)^{\frac{1}{\gamma-1}}}-1\right]$
(6)

For an “infinite” Mach number this becomes:

$\begin{split}\frac{C_{p,\,max}}{\lim\limits_{M_{\infty}\to\infty}}&=\frac{1}{{\frac{\gamma}{2}}\left(\frac{M_{\infty}}{\lim\limits_{M_{\infty}\to\infty}}\right)^{2}}\left[\frac{\left(\frac{\gamma+1}{2}\left(\frac{M_{\infty}}{\lim\limits_{M_{\infty}\to\infty}}\right)^{2}\right)^{\left(\frac{\gamma}{\gamma-1}\right)}}{\left(\frac{2\gamma}{\gamma+1}\left(\frac{M_{\infty}}{\lim\limits_{M_{\infty}\to\infty}}\right)^{2}-\frac{\gamma-1}{\gamma+1}\right)^{\frac{1}{\gamma-1}}}-1\right]\\\
&=\frac{1}{{\frac{\gamma}{2}}\left(\frac{M_{\infty}}{\lim\limits_{M_{\infty}\to\infty}}\right)^{2}}\left[\frac{\left(\frac{\gamma+1}{2}\left(\frac{M_{\infty}}{\lim\limits_{M_{\infty}\to\infty}}\right)^{2}\right)^{\left(\frac{\gamma}{\gamma-1}\right)}}{\left(\frac{2\gamma}{\gamma+1}\left(\frac{M_{\infty}}{\lim\limits_{M_{\infty}\to\infty}}\right)^{2}\right)^{\left(\frac{1}{\gamma-1}\right)}}\right]\\\
&=\frac{\frac{2}{\gamma}{\left(\frac{\gamma+1}{2}\right)}^{\left(\frac{\gamma}{\gamma-1}\right)}}{{\left(\frac{2\gamma}{\gamma+1}\right)}^{\left(\frac{1}{\gamma-1}\right)}}\left[\frac{\left(\frac{M_{\infty}}{\lim\limits_{M_{\infty}\to\infty}}\right)^{\left(\frac{2\gamma}{\gamma-1}\right)}}{\left(\left(\frac{M_{\infty}}{\lim\limits_{M_{\infty}\to\infty}}\right)^{2(\gamma-1)}\left(\frac{M_{\infty}}{\lim\limits_{M_{\infty}\to\infty}}\right)^{2}\right)^{\frac{1}{\gamma-1}}}\right]\end{split}$
(7)

Collecting exponents on
$\frac{M_{\infty}}{\lim\limits_{M_{\infty}\to\infty}}$:

$\begin{split}\frac{C_{p,\,max}}{\lim\limits_{M_{\infty}\to\infty}}&=\frac{\frac{2}{\gamma}{\left(\frac{\gamma+1}{2}\right)}^{\left(\frac{\gamma}{\gamma-1}\right)}}{{\left(\frac{2\gamma}{\gamma+1}\right)}^{\left(\frac{1}{\gamma-1}\right)}}\left[\frac{M_{\infty}}{\lim\limits_{M_{\infty}\to\infty}}^{\left(\frac{2\gamma}{\gamma-1}\right)-(2(\gamma-1)+2)\left(\frac{1}{\gamma-1}\right)}\right]\\\
&=\frac{\frac{2}{\gamma}{\left(\frac{\gamma+1}{2}\right)}^{\left(\frac{\gamma}{\gamma-1}\right)}}{{\left(\frac{2\gamma}{\gamma+1}\right)}^{\left(\frac{1}{\gamma-1}\right)}}\left[\left(\frac{M_{\infty}}{\lim\limits_{M_{\infty}\to\infty}}\right)^{0}\right]\\\
&=\frac{\frac{2}{\gamma}{\left(\frac{\gamma+1}{2}\right)}^{\left(\frac{\gamma}{\gamma-1}\right)}}{{\left(\frac{2\gamma}{\gamma+1}\right)}^{\left(\frac{1}{\gamma-1}\right)}}\\\
\end{split}$ (8)

Then doing some manipulations:

$\begin{split}\frac{C_{p,\,max}}{\lim\limits_{M_{\infty}\to\infty}}&=\frac{1}{\gamma}\cdot\frac{2}{2^{\left(\frac{\gamma}{\gamma-1}\right)}\cdot
2^{\left(\frac{1}{\gamma-1}\right)}}\cdot\frac{{(\gamma+1)^{\left(\frac{\gamma}{\gamma-1}\right)}}}{{\left(\frac{\gamma}{\gamma+1}\right)}^{\left(\frac{1}{\gamma-1}\right)}}\\\
&=\frac{1}{\gamma}\cdot\frac{1}{2^{\left(\frac{2}{\gamma-1}\right)}}\cdot\frac{(\gamma+1)^{\left(\frac{\gamma}{\gamma-1}\right)}}{{\left(\frac{\gamma}{\gamma+1}\right)}^{\left(\frac{1}{\gamma-1}\right)}}\\\
&=\frac{1}{\gamma^{\left(1+\frac{1}{\gamma-1}\right)}}\cdot\frac{1}{2^{\left(\frac{2}{\gamma-1}\right)}}(\gamma+1)^{\left({\frac{\gamma}{\gamma-1}}+\gamma-1\right)}\\\
&=\frac{\gamma^{\left(\frac{\gamma}{1-\gamma}\right)}}{2^{\left(\frac{2}{\gamma-1}\right)}}(\gamma+1)^{\left(\frac{\gamma+1}{\gamma-1}\right)}\end{split}$
(9)

If we substitute $\gamma$=1.4 in Eq.9 then $C_{p\,max}$ becomes:

$\begin{split}\frac{C_{p,\,max}}{\lim\limits_{M_{\infty}\to\infty}}=1.8394\end{split}$
(10)

## References

  * Morgan [2001] Morgan, R. G., _Free Piston-Driven Reflected Shock Tunnels_ , Academic Press, Burlington, VA, 2001, Chap. 4.2, pp. 587–601.  10.1016/B978-012086430-0/50013-0.
  * Ames research staff [1953] Ames research staff, “Equations, tables, and charts for compressible flow,” Technical report NACA RP 1135, NACA, 1953.
  * Anderson, J. D. [2006] Anderson, J. D., _Hypersonic and High Temperature Gas Dynamics_ , 2nd ed., American Institute of Aeronautics and Astronautics, 2006.  10.2514/4.861956.
  * Fuller [1961] Fuller, F. B., “Numerical solutions for supersonic flow of an ideal gas around blunt two-dimensional bodies,” Technical report NASA TN D791, NASA, 1961\.
  * Kantrowitz [1946] Kantrowitz, A., “Heat-capacity lag in gas dynamics,” _The Journal of Chemical Physics_ , Vol. 14, No. 3, 1946, pp. 150–164.  10.1063/1.1724115.
  * Huber and Kantrowitz [1947] Huber, P., and Kantrowitz, A., “Heat-Capacity Lag Measurements in Various Gases,” _The Journal of Chemical Physics_ , Vol. 15, No. 5, 1947, pp. 275–284.  10.1063/1.1746496.
  * Wyllie [2013] Wyllie, E., “Investigation of Hypersonic and High Temperature Gas Effects on Pitot Pressure,” Undergraduate thesis, The University of Queensland, 2013.
  * Candler et al. [2015] Candler, G. V., Johnson, H. B., Nompelis, I., Gidzak, V. M., Subbareddy, P. K., and Barnhardt, M., “Development of the US3D Code for Advanced Compressible and Reacting Flow Simulations,” _53rd AIAA Aerospace Sciences Meeting_ , 2015.  10.2514/6.2015-1893, AIAA Paper 2015–1893.
  * Gibbons et al. [2023] Gibbons, N. N., Damm, K. A., Jacobs, P. A., and Gollan, R. J., “Eilmer: an open-source multi-physics hypersonic flow solver,” _Computer Physics Communications_ , Vol. 282, 2023.  10.1016/j.cpc.2022.108551.
  * Gordon and McBride [1994] Gordon, S., and McBride, B., “Computer program for complex chemical equilibrium compositions and applications,” Technical report NASA-RP-1311, NASA, 1994.
  * Gupta et al. [1990] Gupta, R. N., Yos, J. M., Thompson, R. A., and Lee, K., “A review of reaction rates and thermodynamic and transport properties for an 11-species air model for chemical and thermal nonequilibrium calculations to 30000 K,” Technical report NASA-RP-1232, NASA, 1990.
  * Wilke [1950] Wilke, C. R., “A viscosity equation for gas mixtures,” _Pure and Applied Chemistry_ , Vol. 23, No. 1, 1950, pp. 73–98.
  * Chan et al. [2018] Chan, W. Y. K., Jacobs, P. A., Smart, M., Grieve, S., Craddock, C., and Doherty, L., “Aerodynamic Design of Nozzles with Uniform Outflow for Hypervelocity Ground-Test Facilities,” _Journal of Propulsion and Power_ , Vol. 34, No. 6, 2018, pp. 1467–1478.  10.2514/1.B36938.
  * Jacobs et al. [2020] Jacobs, P. A., Chan, W. Y. K., Sopek, T., Zander, F., Damm, K. A., Gibbons, N. N., and Gollan, R. J., “Calculation of Test Flow Conditions for Reflected-Shock Tunnels,” , 2020.  URL https://gdtk.uqcloud.net/pdfs/eilmer-talk-pj-2020-dec.pdf, Accessed 2023-03-30.
  * Wilcox [2006] Wilcox, D. C., _Turbulence Modelling for CFD_ , 3rd ed., DCW Industries, Inc., La Cañada, CA, 2006.
  * Nichols and Nelson [2004] Nichols, R. H., and Nelson, C., “Wall function boundary conditions including heat transfer and compressibility,” _AIAA Journal_ , Vol. 42, No. 6, 2004, pp. 1107–1114.  10.2514/1.3539.
  * Collen et al. [2021] Collen, P., Doherty, L. J., Subiah, S.-D., Sopek, T., Jahn, I., Gildfind, D., Penty Geraets, R., Gollan, R. J., Hambidge, C., Morgan, R. G., and McGilvray, M., “Development and Commissioning of the T6 Stalker Tunnel,” _Experiments in Fluids_ , Vol. 62, 2021.  10.1007/s00348-021-03298-1.
  * McGilvray et al. [2009] McGilvray, M., Jacobs, P. A., Morgan, R. G., Gollan, R. J., and Jacobs, C. M., “Helmholtz Resonance of Pitot Pressure Measurements in Impulsive Hypersonic Test Facilities,” _AIAA Journal_ , Vol. 47, No. 10, 2009, pp. 2430–2439.  10.2514/1.42543.
  * Jacobs et al. [2011] Jacobs, P. A., Gollan, R. J., Potter, D. F., Zander, F., Gildfind, D. E., Blyton, P., Chan, W. Y. K., and Doherty, L. J., “Estimation of high-enthalpy flow conditions for simple shock and expansion processes using the ESTCj program and library,” techreport, The University of Queensland, 2011\.
  * Hannemann et al. [2010] Hannemann, K., Martinez Schramm, J., Wagner, A., Karl, S., and Hannemann, V., “A Closely Coupled Experimental and Numerical Approach for Hypersonic and High Enthalpy Flow Investigations Utilising the HEG Shock Tunnel and the DLR TAU Code,” Tech. report RTO-EN-AVT-186, NATO, 2010.


# Collective Free Electron Excitations in Half-Space Configuration

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

In this research, using the suitably linearized kinetic-corrected Schrödinger-
Poisson system, we study the first-order density and potential perturbations
in a half-space electron gas of arbitrary degeneracy in the jellium model. By
solving the eigenvalue problem in coupled pseudoforce model with appropriate
boundary conditions, pure states of plasmonic excitations in the quasineutral
electron gas are determined from which mixed states of the system
characterizing the thermodynamic equilibrium is constructed by enumeration of
quasiparticle orbital. It is shown that the Friedel-like electron density
fringes appear inside the electron gas due to interference effect. It is also
remarked that screened electron packet appears in front of the physical
boundary which is different from continuous density spectrum. The
electrostatic potential coupled to the electron density beyond the physical
boundary becomes oscillatory for density-temperature regimes relevant to fully
degenerate electron gas. The potential minimum found in front of physical
boundary confirms Casimir-Polder-like force between parallel metallic plates
and may be responsible for agglomeration of metallic nanoparticles in quantum
contacts. The dependence of potential and electron density at the electron gas
boundary is investigated for a wide range of parametric density-temperature
regimes. It is also found that the kinetic correction to the model has a
profound effect on the particle-like low phase-speed matter-wave dispersion
branch of plasmon excitations, where, the increase in the density of electrons
leads to the increase in the effective electron mass, and consequently,
decrease in the electron mobility, while, increase in the temperature causes
an increase in the effective electron mobility.

###### pacs:

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

##  I Introduction

Plasmonics man1 ; maier  is a new interdisciplinary field of research with
important applications in nanoelectronics haug ; gardner , optoelectronics
wooten , and semiconductor integrated circuit industry markovich ; manvir .
Plasmons are elementary collective excitations of quantum plasmas which play
an inevitable role in physical phenomena krall ; stenf1 . These entities
provide ideal platform for fast THz device communications umm  and beyond
where the conventional wired communications fail to operate. Plasmonic energy
conversion is an alternative way of solar power extraction in future
photovoltaic and catalytic designs due to their much higher efficiency
compared to semiconductor technologies cesar ; jac . The state of the art
nanocatalitic engineering for energy extraction is based on conversion of
energy of collective surface electron plasma oscillations caused by localized
electromagnetic radiation into hot electron generation in plasmonic devices
which operates in ultraviolet-visible (UV-VIS) frequency range. However, the
conventional solar cell devices are based on low energy interband electron-
hole transitions in semiconductors which operate at visible-near infrared
(VIS-NIR) low frequency range. However, plasmonic devices operating in
ultraviolet spectrum in the wavelength range ($10-380$)nm requires special
technological considerations and nanofabrication costs at  due to low
contribution of plasmonic device operation range to the solar radiation
spectrum. Recent studies have revealed that nanostructured materials can be
used to directly convert collected electromagnetic energy into electricity by
local surface plasmon resonance (LSPR) in the hot electron generation
mechanism at2 . The generated hot electrons are then collected by putting the
metallic nanoparticles in contact with a good electron acceptor like $TiO_{2}$
in Schottky configuration tian . Studies also reveal that hot electron
extraction efficiency is strongly dependent on the size and geometry of
nanoparticles as well as their composition vin .

A wide range of quantum plasma regime extends from doped semiconductors hu1 ;
seeg  with moderately high electron number density at low temperature to high
density metallic compounds up to high-temperature and density so-called the
warm dense matter (WDM) glenzer2 ; koenig . Quantum effect arise naturally
when the interspecies distances are lowered beyond the thermal de Broglie’s
wavelength bonitz . In inertial confinement fusion experiment quantum effects
vary in strength due to change in equation of state (EoS) of matter during
shock compression causing increase in temperature and number density mithen .
The change of EoS has fundamental effect on the faith of compact stellar
object setting distinct limit on the mass of such solar entities
chandrasekhar1 . Another interesting feature of charged particles is their
complex collective electromagnetic interactions with their average
interparticle distances beyond the quantum limit. The later is due to the
nonlocal nature of interactions via the Bohm’s quantum potential. Such a
collective interactions lead to many interesting features of electron
interactions with ionic lattice in solids kit ; ash . The theoretical
development of quantum plasmas has a long history with pioneering
contributions over the past century bohm ; pines ; klimontovich . Because of
the dominant quantum effects caused by the EoS of a Fermi gas, semiclassical
plasma theories which incorporate the quantum statistical pressure effects
ignoring the electron diffraction can lead to some interesting collective
features in dense plasmas ichimaru1 .

There has been an increasing interest in the quantum plasma research over the
recent decade due to their better description of collective features in dense
plasmas se ; sten ; ses ; haas1 ; brod1 ; mark1 ; man3 ; mold2 ; sm ; fhaas .
Quantum kinetic man2 , quantum (magneto)hydrodynamic haasbook  and gradient
corrected theories bonitz0  has been used to investigate different aspects of
many-fermion interactions in plasmas. Application of quantum hydrodynamic has
led to discovery of many interesting collective properties of plasmas which
has not been captured by theories ignoring the quantum electron diffraction
feature. The quantum hydrodynamic formulations is cast into the Schrödinger-
Poisson model manfredi  using the Madelung transformation madelung . The leter
model has been used to study various linear and nonlinear features of plasmon
excitations in quantum plasmas mannew ; hurst . Recently, the pseudoforce
system derived from linearized Schrödinger-Poisson model, has been used to
study different features of pure states of elementary excitations in an
electron gas of arbitrary degenerate electron gas akbquant . The pseudoforce
model has a fundamental property of taking into account the dual lengthscale
nature of plasmons in the electron gas. The discrepancy between quantum
hydrodynamic and density functional theories in the linear limit akbhd ; mold
; michta  has led to some kinetic correction to the fluid theories in the
quantum potential term. In current analysis we take into account the kinetic
correction in the Schrödinger-Poisson model in order study the plasmonic
excitations in the half-space electron gas and investigate the effects of such
density-temperature dependent correction on the energy dispersion of plasmonic
excitations.

##  II The Mathematical Model

The closed set of quantum hydrodynamic equations including the kinetic
correction to the electron quantum diffraction effect is as follows

$\displaystyle\frac{{\partial n}}{{\partial t}}+\frac{{\partial nu}}{{\partial
x}}=0,$ (1a) $\displaystyle\frac{{\partial u}}{{\partial t}}+u\frac{{\partial
u}}{{\partial x}}=\frac{e}{m}\frac{{\partial\phi}}{{\partial
x}}-\frac{1}{m}\frac{{\partial\mu}}{{\partial
x}}+{\rm{}}\frac{{{\zeta\hbar^{2}}}}{{2{m^{2}}}}\frac{\partial}{{\partial
x}}\left({\frac{1}{{\sqrt{n}}}\frac{{{\partial^{2}}\sqrt{n}}}{{\partial{x^{2}}}}}\right),$
(1b) $\displaystyle\frac{{{\partial^{2}}\phi}}{{\partial{x^{2}}}}=4\pi
e\left({n-{n_{0}}}\right),$ (1c)

in which $n$, $u$, $\phi$ and $\mu$ denote the electron number density, fluid
velocity, electrostatic potential and chemical potential. Here $\zeta$
represents the temperature and density dependent correction factor in order to
account for the correct contribution to the low-frequency and long-wavelength
plasmon excitationshaas2016 .

$\zeta=\frac{{{\rm{L}}{{\rm{i}}_{3/2}}\left[{-\exp(\beta{\mu_{0}})}\right]{\rm{L}}{{\rm{i}}_{-1/2}}\left[{-\exp(\beta{\mu_{0}})}\right]}}{3{{\rm{L}}{{\rm{i}}_{1/2}}{{\left[{-\exp(\beta{\mu_{0}})}\right]}^{2}}}},$
(2)

in which $\beta=1/k_{B}T$ and $\mu_{0}$ is the equilibrium value of the
chemical potential. The function ${\rm{Li}}_{k}$ is the polylog function
defined below

${\rm{L}}{{\rm{i}}_{k}}[-\exp(\eta)]=-\frac{{{F_{k-1}}(\eta)}}{{\Gamma(k)}},$
(3)

where $\Gamma$ is the gamma function and the polylog function definition in
terms of the Fermi function is

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

The one-dimensional Schrödinger-Poisson model manfredi  is mathematically
equivalent to the quantum hydrodynamic equations (1).

$\displaystyle i\sqrt{\zeta}\hbar\frac{{\partial\cal N}}{{\partial
t}}=-\zeta\frac{{{\hbar^{2}}}}{{2m}}\frac{{\partial^{2}{{\cal
N}}}}{{\partial{x^{2}}}}-e\phi{\cal N}+\mu(n,T){\cal N},$ (5a)
$\displaystyle\frac{{\partial^{2}{\phi}}}{{\partial{x^{2}}}}=4\pi e(|{\cal
N}|^{2}-n_{0}),$ (5b)

where ${\cal N}=\sqrt{n(x,t)}\exp[iS(x,t)/\hbar\sqrt{\zeta}]$ and
$u(x,t)=(1/m)\partial S(x,t)/\partial x$ are the Madelung transformations
madelung ; bohm0 . Also, $n_{0}$ is the unperturbed electron density which is
equal to that of the neutralizing positive background. It is straightforward
to show the equivalence of (5) and (1) by using the ansatz transformations.
Using ${\cal N}=\sqrt{n(x,t)}\exp[iS(x,t)/\hbar\sqrt{\zeta}]$ in 5(a) and
separating real and imaginary parts, we have

$\displaystyle m\frac{{\partial n(x,t)}}{{\partial t}}+\frac{{\partial
n(x,t)}}{{\partial x}}\frac{{\partial S(x,t)}}{{\partial
x}}+n(x,t)\frac{{{\partial^{2}}S(x,t)}}{{\partial{x^{2}}}}=0,$ (6a)
$\displaystyle\frac{{{\partial^{2}}S(x,t)}}{{\partial t\partial
x}}+\frac{1}{m}\frac{{\partial S(x,t)}}{{\partial
x}}\frac{{{\partial^{2}}S(x,t)}}{{\partial{x^{2}}}}=\frac{{e\partial\phi(x,t)}}{{\partial
x}}-\frac{{\partial\mu(x)}}{{\partial x}}+\frac{{\partial B(x,t)}}{{\partial
x}},$ (6b) $\displaystyle
B(x,t)=\frac{{\zeta{\hbar^{2}}}}{{8m{n^{2}}(x,t)}}\left\\{{2n(x,t)\frac{{{\partial^{2}}n(x,t)}}{{\partial{x^{2}}}}-{{\left[{\frac{{\partial
n(x,t)}}{{\partial x}}}\right]}^{2}}}\right\\}.$ (6c)

Now using $u(x,t)=(1/m)\partial S(x,t)/\partial x$, with little mathematical
work, one arrives at the continuity and momentum equations in (1). The
chemical potential, $\mu(n,T)$ of the arbitrary degenerate electron gas is
defined through the generalized parametric equation of state (EoS). The EoS
for isothermal compression of electron gas with arbitrary degeneracy is then
given in terms of the Fermi integrals

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

Note that in current analysis we use the Thomas-Fermi assumption where the
temperature is fixed in accordance with the Fermi-Dirac single-particle
distribution dya1 ; dya2 .

The main purpose of current research is to use the model for metallic
compounds with fully degenerate electrons and strongly coupled ions in static
background lattice, i.e. jellium model. The degeneracy parameter is
$\eta=T/T_{F}$ in which $T$ and $T_{F}$ are respectively the temperatures for
electron gas and Fermi energy level. Therefore, the parametric regime
$\eta<1$($\eta>1$) denotes the degenerate(nondegenrate) electron gas regime.
For fully degenerate metals or metallic nanoparaticles, we have the extreme
case of $\eta\ll 1$, or equivalently, $T_{F}\gg T$. For elemental metals at
room temperature, $E_{F}$ amounts to about $10^{4}$K or slightly higher,
depending on their electron concentrations. In this limit the chemical
potential equals the Fermi energy and all electrons reside below the well
defined Fermi level. On the other hand, ion dynamic properties are
characterized via the ion coupling parameter, $\Gamma=Ze^{2}/dk_{B}T_{i}$, in
which $d$ is the average inter-ion distances and $T_{i}$ is the ion
temperature which is much lower compared to that of electron gas, $T_{i}\ll
T_{e}$, due to the large mass ratio. A good criterion for the weak and strong
coupling cases are $\Gamma\ll 1$ and $\Gamma\gg 1$, respectively ichimaru0 .
Therefore, metallic compounds with $r_{s}=2-6$, where $r_{s}=d/r_{B}$ being
the Bruckener parameter and $r_{B}$ the Bohr radius, are categorized as
strongly coupled material with ions fixed in crystal lattice. For simplicity,
in current model we ignore the effect of phonon excitations on plasmonic
properties of the electron gas. Moreover, the electron fluid with $r_{s}\ll
10^{-2}$ is regarded as nonrelativistic electron gas, while semiconductors
have the parameter value of $r_{s}>25$. Because for semiconductors the Fermi
temperature is very close to that of the electrons, these material are
regarded as partially degenerate electron gas. The electron degeneracy starts
at electron number density of $n\simeq 10^{18}$cm-3 where the free electron
wavefunctions start to overlap. However, the electron density of doped
semiconductors can be much lower than this value. In fully degenerate
elemental metals the electron concentration is typically in the range
($10^{21}-10^{23}$)cm-3 characterizing them as fully degenerate quasineutral
electron gas.

Here we would like to give a brief comment on applicability of Schrödinger-
Poisson model to study fermion quasiparticle systems, which was first
developed by Manfredi and Haas from quantum hydrodynamic theory and Wigner-
Poisson system manfredi . The quantum hydrodynamic equations are obtained from
moments of the Wigner distribution function under liming assumptions. They
used a first-order plane-wave approximation for N-body wavefunctions with
amplitudes independent of quasiparticle orbital. Therefore, the only orbital
dependent part comes from the phase function associated with each orbital.
However, such assumption may put limitations on the applicability of the model
to systems in thermal contact with surrounding environment, since, the density
matrix which is used to study the temporal evolution becomes independent of
plasmon orbital and only depends on spacial density variations. Bonitz et. al.
bonitz0  have given an extensive discussion on different formulations and
possible improvements of the quantum hydrodynamic model. It has been noted
that, due to the mentioned limitations, quantum hydrodynamic model gives an
expression for plasmon dispersion which differs from that of the random phase
approximation (RPA) in the low-frequency static limit michta ; akbhd ; mold .
Therefore, a kinetic correction factor, $\zeta$, is applied to the quantum
force term in order for the hydrodynamic theory to be consistent with RPA. In
the proceeding analysis, we use the kinetic-corrected Schrödinger-Poisson
model in order to obtain statefunctions of pure quasiparticle orbital. The
energy eigenvalues of these orbital are used to construct appropriate mixed
states using quasiparticle orbital enumeration and determine thermodynamic
quantities, exactly in a similar manner as is done in the free electron theory
of solids ash .

##  III Pure Quasiparticle States

In order to study the system (5) in the linear limit, for simplicity we ignore
minor interaction terms such as the electron exchange and correlations.
However, we will account to antisymmetric wavefunctions later on in this
analysis. Appropriately normalized time-independent half-space jellium model
of quasineutral electron gas of arbitrary degeneracy, constitutes of a coupled
system of second-order differential equations, the so-called pseudoforce
system

$\displaystyle\zeta\frac{{d^{2}{\Psi(x)}}}{{d{x^{2}}}}+\Phi(x)+E\Psi(x)+E=0,$
(8a) $\displaystyle\frac{{d^{2}{\Phi(x)}}}{{d{x^{2}}}}-{\Psi}(x)=0,$ (8b)

where we have used a variable separation condition, ${\cal
N}(x,t)=\psi(x)\psi(t)$, for quasi-stationary collective excitations and
applied the normalization scheme, $\Psi(x)=\psi(x)/\sqrt{n_{0}}$ and
$\Phi(x)=e\phi(x)/E_{p}$ with $E_{p}=\hbar\sqrt{4\pi e^{2}n_{0}/m}$ being the
plasmon energy of the system. Also, $E=(\epsilon-\mu_{0})/E_{p}$ with
$\epsilon$ being the energy eigenvalue of the collective quasiparticle
(plasmons) as defined through
$\epsilon\psi(t)=\hbar\omega\psi(t)=i\sqrt{\zeta}\hbar d\psi(t)/dt$ with
$\omega$ being the eigenfrequency of quasiparticle excitations. Note that the
space and time variables are also normalized, respectively, to the plasmon
length $1/k_{p}$ with $k_{p}=\sqrt{2mE_{p}}/\hbar$ being the plasmon
wavenumber and $\hbar/E_{p}$. Also, the Thomas-Fermi assumption is used for
the chemical potential in the fully degenerate limit, where, the chemical
potential variations are compensated by the local electrostatic potential of
the electron gas. The solution to psudoforce model, (8), may be written as

$\left[{\begin{array}[]{*{20}{c}}{\Phi(x)}\\\
{\Psi(x)}\end{array}}\right]=\frac{{A\zeta}}{\alpha}\left\\{{\begin{array}[]{*{20}{c}}{{\Psi_{0}}+k_{2}^{2}\left({{\Phi_{0}}+E}\right)}&{-\left[{{\Psi_{0}}+k_{1}^{2}\left({{\Phi_{0}}+E}\right)}\right]}\\\
{-\left[{\left({{\Phi_{0}}+E}\right)/\zeta+k_{1}^{2}{\Psi_{0}}}\right]}&{\left({{\Phi_{0}}+E}\right)/\zeta+k_{2}^{2}{\Psi_{0}}}\end{array}}\right\\}\left[{\begin{array}[]{*{20}{c}}{\cos({k_{1}}x)}\\\
{\cos({k_{2}}x)}\end{array}}\right],$ (9)

in which $\Phi_{0}$ and $\Psi_{0}$ characterize the initial functional values
for the homogenous system and $A$ is a normalizing factor, where, we assumed
$d\Psi(x)/dx|_{x=0}=d\Phi(x)/dx|_{x=0}=0$, for simplicity. The characteristic
excitation wavenumbers $k_{1}$ and $k_{2}$ are given as

${k_{1}}=\sqrt{\frac{{E-\alpha}}{{2\zeta}}},\hskip
8.53581pt{k_{2}}=\sqrt{\frac{{E+\alpha}}{{2\zeta}}},\hskip
8.53581pt\alpha=\sqrt{{E^{2}}-4\zeta}.$ (10)

Note the complementarity-like relation between the wave- and particle-like
excitation wavenumbers, $k_{1}k_{2}=1$. The solution (9) leads to the matter-
wave energy dispersion $E=(1+\zeta k^{4})/k^{2}$ in which $E$ and $k$ are
normalized to $E_{p}$ and plasmon wavenumber $k_{p}$, respectively. The
solution (9) has the same characteristic as the proposed de Broglie’s double
solution in the pilot-wave theory bohm0 ; nikolic , where a single electron is
guided through by collective electrostatic interactions among other electrons.
The kinetic correction parameter depends on normalized density and temperature
parameters

$\zeta(\sigma,\theta)=\frac{{{\rm{L}}{{\rm{i}}_{3/2}}\left[{-\exp(\sigma/\theta)}\right]{\rm{L}}{{\rm{i}}_{-1/2}}\left[{-\exp(\sigma/\theta)}\right]}}{3{{\rm{L}}{{\rm{i}}_{1/2}}{{\left[{-\exp(\sigma/\theta)}\right]}^{2}}}},$
(11)

in which $\sigma=\mu_{0}/E_{p}$, $\theta=T/T_{p}$ are fractional chemical
potential and electron temperature parameters with the plasmon temperature
defined as, $T_{p}=E_{p}/k_{B}$.

##  IV Damped Quasiparticle States

The quasipaticle excitations in (8) is now generalized to include the spacial
damping or screening effect at the jellium boundary, $x=0$. To this end, we
consider the damped pseudoforce system as follows akbnew .

$\displaystyle\zeta\frac{{{d^{2}}\Psi(x)}}{{d{x^{2}}}}+2\zeta\xi\frac{{d\Psi(x)}}{{dx}}+\Phi(x)+{E}\Psi(x)+E=0,$
(12a)
$\displaystyle\frac{{{d^{2}}\Phi(x)}}{{d{x^{2}}}}+2\xi\frac{{d\Phi(x)}}{{dx}}-\Psi(x)=0,$
(12b)

where $\xi=k_{sc}/k_{p}$ is the normalized one-dimensional screening parameter
depending on the fractional parameters via,
$\xi^{2}=({E_{p}}/2{n_{0}})\partial
n/\partial\mu=(1/2\theta){\rm{L}}{{\rm{i}}_{1/2}}\left[{-\exp({\sigma}/\theta)}\right]/{\rm{L}}{{\rm{i}}_{3/2}}\left[{-\exp({\sigma}/\theta)}\right]$
akbnew . The general solution to damped pseudoforce system (12) is

$\displaystyle\begin{array}[]{*{20}{l}}{{\Phi_{d}}(x)=\frac{{A\zeta{{\rm{e}}^{-\xi
x}}}}{\alpha}\left\\{{\begin{array}[]{*{20}{l}}{\left[{k_{2}^{2}\left({{\Phi_{0}}-E}\right)+{\Psi_{0}}}\right]\left[{\cos({\beta_{1}}x)+\frac{\xi}{{{\beta_{1}}}}\sin({\beta_{1}}x)}\right]-}\\\
{\left[{k_{1}^{2}\left({{\Phi_{0}}-E}\right)+{\Psi_{0}}}\right]\left[{\cos({\beta_{2}}x)+\frac{\xi}{{{\beta_{2}}}}\sin({\beta_{2}}x)}\right]}\end{array}}\right\\},}\\\
{{\Psi_{d}}(x)=\frac{{A\zeta{{\rm{e}}^{-\xi
x}}}}{\alpha}\left\\{{\begin{array}[]{*{20}{l}}{\left[{\left({{\Phi_{0}}-E}\right)/\zeta+k_{2}^{2}{\Psi_{0}}}\right]\left[{\cos({\beta_{2}}x)+\frac{\xi}{{{\beta_{2}}}}\sin({\beta_{2}}x)}\right]-}\\\
{\left[{\left({{\Phi_{0}}-E}\right)/\zeta+k_{1}^{2}{\Psi_{0}}}\right]\left[{\cos({\beta_{1}}x)+\frac{\xi}{{{\beta_{1}}}}\sin({\beta_{1}}x)}\right]}\end{array}}\right\\},}\end{array}$
(13g)

where $\beta_{1}=\sqrt{k_{1}^{2}-\xi^{2}}$ and
$\beta_{2}=\sqrt{k_{2}^{2}-\xi^{2}}$. The boundary values $\Phi_{0}$ and
$\Psi_{0}$ are defined using the half-space electron density distribution
condition. For the spatial region, $x\leq 0$, we use the solution (9), while,
for $x>0$ the solution (13) is used. Let us assume $\Phi(x=0)=0$ at the
boundary which leads to $\Phi_{0}=-E$. On the other hand, we assume that the
electron gas is extended to $-L$, thus, leading to equation

$\int\limits_{-L}^{0}{A|\Psi(x){|^{2}}dx=L},$ (14)

in which $A$ is a dimensionless normalization factor for both $\Psi(x)$ and
$\Phi(x)$. The condition (14) gives rise to the following expression in terms
of $A$ and $\Psi_{0}$

$\displaystyle\frac{{4k_{1}^{2}k_{2}^{3}\cos\left({{k_{1}}L}\right)\sin\left({{k_{2}}L}\right)}}{{L\left({k_{1}^{2}-k_{2}^{2}}\right)\left({k_{1}^{4}+k_{2}^{4}}\right)}}-\frac{{4k_{1}^{3}k_{2}^{2}\cos\left({{k_{2}}L}\right)\sin\left({{k_{1}}L}\right)}}{{L\left({k_{1}^{2}-k_{2}^{2}}\right)\left({k_{1}^{4}+k_{2}^{4}}\right)}}+\frac{{\left({k_{1}^{3}+k_{2}^{3}}\right)\sin\left({2{k_{1}}L}\right)}}{{2L\left({k_{1}^{4}+k_{2}^{4}}\right)}}=\frac{{2{\alpha^{2}}{{\left({k_{1}^{4}+k_{2}^{4}}\right)}^{-1}}}}{{{A^{2}}{\zeta^{2}}\Psi_{0}^{2}}}.$
(15a)

In the limit large system size as compared to plasmon length, $L\gg 1$, we
arrive at a simplified expression

$A(E,\sigma,\theta)=\frac{{\alpha\sqrt{2}}}{{\zeta(\sigma,\theta){\Psi_{0}}\sqrt{k_{1}^{4}+k_{2}^{4}}}}.$
(16)

We have taken $L=100$ throughout this analysis. It is found that the solution
for the undamped region, $x<0$, reduces to that of the electron gas confined
to infinite potential akbquant

$\left[{\begin{array}[]{*{20}{c}}{{\Phi}(x)}\\\
{{\Psi}(x)}\end{array}}\right]=\sqrt{\frac{2}{{k_{1}^{4}+k_{2}^{4}}}}\left[{\begin{array}[]{*{20}{c}}1&{-1}\\\
{-k_{1}^{2}}&{k_{2}^{2}}\end{array}}\right]\left[{\begin{array}[]{*{20}{c}}{\cos({k_{1}}x)}\\\
{\cos({k_{2}}x)}\end{array}}\right].$ (17)

The energy eigenvalues of the pure quasiparticle quantum states are quantized
with the condition, $k_{2}\pm k_{1}=2\pi n$ (the quantum number $n$ not to be
confused with electron number density) leading to the energy eigenvalues

${E_{n}}=\frac{{1+\zeta{{\left({n\pi+\sqrt{1+{n^{2}}{\pi^{2}}}}\right)}^{4}}}}{{{{\left({n\pi+\sqrt{1+{n^{2}}{\pi^{2}}}}\right)}^{2}}}},$
(18)

with the $n=0$ ground state level corresponding to $k_{1}=k_{2}$ is the
quantum beating level which coincides with the quasiparticle conduction level
minimum at $E_{0}=1+\zeta$.

In order to develop a theory for mixed quasiparticle quantum states we follow
exactly the free electron analogy by numeration of the quasiparticle orbitals
considering that the pure state wavefunctions are complete, orthonormal and
antisymmetric (obtained via N-dimensional Slater determinant) due to the pauli
exclusion acting on electrons. However, an important difference between
electrons and collective fermion quasiparticles is that the later species do
not have spin degeneracy like electrons, hence, each quasiparticle occupies a
distinct energy level. This is implied by the fact that each quasiparticle
orbital corresponds to a free electron (spin-up and spin-down) level akbquant
, so that, that metallic density is characterized by the resonant condition,
$E_{p}=2E_{F}$ akbground . In other words, the change in single-electron spin
in antisymmetric N-electron wavefunction modifies the corresponding
quasiparticle orbital to a new level.

##  V Mixed Quasiparticle States

Once the pure quasiparticle states are known, it is an easy task to compose
appropriate mixed states for equilibrium condition. These states are given as

$\displaystyle\Phi(x)=\sum\limits_{j}{f({E_{j}}){\Phi_{j}}(x)},\hskip
8.53581pt\Psi(x)=\sum\limits_{j}{f({E_{j}}){\Psi_{j}}(x)},$ (19a)
$\displaystyle n(x)=\sum\limits_{j}{{n_{j}}},\hskip
8.53581pt{n_{j}}={g_{j}}|{\Psi_{j}}{|^{2}},$ (19b)

where $f(E_{j})$ is the probability of occupying the $j$-th orbital by a
quasiparticle and $g_{i}=f(E_{j})D(E_{j})/\sum\limits_{j}f(E_{j})D(E_{j})$
($\sum_{j}g_{j}=1$), with $D(E_{j})$ being the density of states at given
energy eigenvalue, $E_{j}$, is the occupation number at each orbital. Note
that $n_{j}$ gives the number density of quasiparticles at given orbital which
is exactly equal to the number of electrons. Note also that the mixed state-
vector, $\\{\Phi(x),\Psi(x)\\}$, satisfies the corresponding pseudoforce
systems (8) and (12). This can be shown by appropriate weighted summation over
orbital energy eigenvalues.

Let us now study the system by enumeration of quasiparticle modes. For $N$
electron orbital there corresponds $N$ quuasiparticle orbital. The
quasiparticle wavenumber modes are given by
$k=(N\pi+\sqrt{1+{N^{2}}{\pi^{2}}})/L$, corresponding to maximum orbital
energy of

${E_{m}}=\frac{{{L^{4}}+\zeta{{\left({N\pi+\sqrt{1+{N^{2}}{\pi^{2}}}}\right)}^{4}}}}{{{L^{2}}{{\left({N\pi+\sqrt{1+{N^{2}}{\pi^{2}}}}\right)}^{2}}}}.$
(20)

The quasiparticle DoS is $dN/dE$ which follows

$D(E)=\frac{{\left({{E^{2}}-4\zeta}\right)\left({{L^{4}}+\zeta}\right)+E\left({{L^{4}}-\zeta}\right)\sqrt{{E^{2}}-4\zeta}}}{{4\pi
L\sqrt{2}\left({{E^{2}}-4\zeta}\right)\sqrt{\zeta\left({{L^{4}}-\zeta}\right)\sqrt{{E^{2}}-4\zeta}-4{L^{2}}{\zeta^{2}}+E\left({{L^{4}}+\zeta}\right)}}}.$
(21)

Therefore, we have $g_{j}=f(E_{j})D(E_{j})/\sum\limits_{j}f(E_{j})D(E_{j})$
where $f(E_{j})=1/[1+\exp(E_{j}/\theta)]$ is the Fermi occupation number which
we assume also applies to quasiparticles for reasons mentioned earlier.
However, in the continuum limit for a large spatial system size compared to
the plasmon length, $L\gg 1$, the electron number density distribution is
given by

$n(x)=\int\limits_{2\sqrt{\zeta}}^{\infty}{D(E)f(E){{\left|{\Psi_{E}(x)}\right|}^{2}}dE},$
(22)

and the electrostatic energy distribution can be obtained also as follows

$\Phi(x)=\int\limits_{2\sqrt{\zeta}}^{\infty}f(E)\Phi_{E}(x)dE.$ (23)

##  VI Numerical Analysis and Discussion

Figure 1: 1(a) The variation of normalized screening (damping) parameter as a
function of normalized chemical potential, $\sigma$ for different values of
normalized electron temperature, $\theta$. 1(b) The variation of kinetic
correction as a function of normalized chemical potential, $\sigma$ for
different values of normalized electron temperature, $\theta$. 1(c) Variation
in normalized temperature $\theta$ in terms of electron number density in
logarithmic scale. 1(d) The plasmon energy variation in terms of the electron
number density in logarithmic scale. The increase in the thickness of curves
in each plot is meant to represent an increase in the varied parameters above
each panel.

Figure 1 shows variations of the characteristic plasmon parameters for the
electron gas of arbitrary degeneracy. Figure 1(a) depicts the variations of
normalized one-dimensional screening parameter, $\xi(\mu,\theta)$ in terms of
the normalized chemical potential, $\mu$, for various values of the normalized
plasmon temperature, $\theta$. It is remarked that with increase in the
chemical potential of the gas the screening parameter decreases substantially.
It is also seen that as the normalized temperature increases, variations in
$\xi$ over the chemical potential range decreases. This plot also shows that
for the case of complete degeneracy the dependence of the normalized screening
parameter to $\theta$ becomes insignificant. Figure 1(b) shows the dependence
of kinetic correction parameter on the normalized chemical potential,
$\sigma$, for different values of the normalized plasmon temperature,
$\theta$. It is remarked that the correction decreases with decrease in
$\sigma$. However, for given value of $\sigma>0$ the correction higher for
larger value of $\theta$. The value of this parameter is always in the range
$1/9<\zeta<0.265$. Figure 1(c) depicts the variations in the fractional
electron temperature, $\theta$, as a function of the electron number density
for various absolute electron temperatures. It is seen that for a given
electron temperature the normalized electron temperature decreases
significantly with increase in the electron number density. The variation over
density becomes even more significant when the electron temperature increases.
It is remarked that in the complete degeneracy limit dependence of this
parameter to electron temperature becomes insignificant. Moreover, Fig. 1(d)
shows the variation in plasmon energy in terms of electron number density in a
logarithmic scale. It is remarked that the increase in plasmon energy becomes
sharp as the degeneracy limit, ($n_{0}\simeq 10^{18}$cm${}^{-3})$, sets in.
For typical metals the plasmon energy amounts a few electron Volts. For
instance, for cesium with the plasmon energy as low as $E_{p}\simeq 2.9$eV it
can reach much higher value for aluminium with $E_{p}\simeq 15$eV. However, in
order to compare the screening effect in these the difference in their
chemical potential which is their Fermi energy at zero-temperature limit
($E^{Cs}_{F}=1.59$eV, $E^{Al}_{F}=11.7$eV) must be taken into account.

Figure 2: Profiles of pure quasiparticle state statefunctions of half-space
electron gas of arbitrary degeneracy for given energy eigenvalue.

Figure 2 shows the variations in pure quasipaticle state functions, $\Psi(x)$,
and electrostatic potential energy, $\Phi(x)$, for given values of normalized
electron density and temperature and different energy eigenvalues. It is
clearly evident that values and derivatives of state functions match at the
half-space boundary. As remarked in Figs. 2(a) and 2(b), for $x<0$ there are
dual-tone oscillations in quasiparticle probability, $\Psi(x)$ due to dual
nature of these entities. Meanwhile, in the screened region, $x>0$, there is
an oscillatory exponential decrease in amplitude of the quasiparticle state,
$\Psi(x)$. In Fig. 2(b) the increase in energy eigenvalue leads to a relative
decrease in the wavefunction amplitude. Local variations in the electrostatic
energy of the system is shown in Figs. 2(c) and 2(d). It seems that the
variations in the region, $x<0$, for $\Phi(x)$ is also dual-tone as compared
to that of $\Psi(x)$. This is due to the fact that plasmon excitations are of
dual-wavenumber character with the larger wavenumber corresponding to the
single particle oscillations and the smaller one to the collective oscillation
in the electron gas. It is therefore seen that at high values of quasiparticle
energy the modulated fine oscillations in electrostatic energy profile
correspond to the single particle excitationsand the larger oscillatory
pattern is due to the collective plasmonic oscillations in the electron gas.
On the other hand, in the region, $x>0$, the electrostatic energy decay has
nearly a monotonic exponential decay without oscillations. Figure 2(d) reveals
that unlike for the case of number density, increase of energy eigenvalue also
significantly lowers the amplitude of electrostatic energy profile. However,
in the region, $x>0$, the decay rate for lower energy eigenvalue in Fig. 2(c)
is slightly faster as compared to that of the high energy profile, shown in
Fig. 2(d).

Figure 3: 3(a) and 3(c) show the normalized perturbed electron number-density
at a thermal equilibrium condition for given parameter values of $\theta$ and
$\mu$. 3(b) and 3(d) show normalized perturbed electrostatic potential energy
at a thermal equilibrium condition for given parameter values of $\theta$ and
$\mu$.

Figure 3 depicts the spatial distribution of normalized perturbed electron
number-density and electrostatic energy for given normalized chemical
potential, $\sigma$, and temperature, $\theta$. Figure 3(a) shows the electron
number density profile for given values of $\mu=0.3$ and $\theta=0.2$. The
Fig. 3(a) illustrates some interesting features of half-space plasmon
excitations at equilibrium temperature. It is seen that a density hump forms
just before the jellium boundary, meanwhile, an electron depleted region
beyond the jellium boundary appears just before a small electron density hump
occurs. Another quite unique feature is the well-defined periodic density
oscillations extended to the electron gas region ($x<0$) which is due to
resonant interaction between single electron and collective excitations. Such
feature has been shown to be a main characteristic feature of an electron gas
confined to an infinite potential well akbquant . The electron packet forming
in front of the jellium boundary differs significantly form previous reported
spill-out effect ciraci ; mario2 . Figure 3(b) depicts the electrostatic
energy profile corresponding to the parameter values used in Fig. 3(a). It is
seen that the periodic structure is also present in this case. Close to the
jellium boundary the amplitude of electrostatic energy is significantly
increased, while, beyond the boundary the potential energy exponentially
decreases towards infinity. Figures 3(c) and 3(d) depict the equilibrium mixed
statefunction profiles for higher value of fractional chemical potential. It
is evident from Fig. 3(c) that the density hump before physical boundary has
become sharper. Note the formation of a surface dipole before the jellium
boundary due to electron depletion befor the physical edge. It is furthe
remarked that density of the electron packet beyond the physical edge
increases significantly due to increase of the normalized chemical potential.
Quite similar trend appears for the electrostatic energy profile as the value
of $\sigma$ increases.

Figure 4: 4(a) and 4(c) show the variations of normalized equilibrium electron
number density and 4(b) and 4(d) show the normalized electrostatic potential
energy beyond the jellium boundary for different values of the normalized
electron temperature and chemical potential. Increase in thickness of curves
in each plot represents an increase in the varied parameter above each panel.

Figure 4 depicts the variations of number density and electrostatic energy
beyond the physical boundary, $x>0$. Figure 4(a) shows that the electron
density extends to the screened area beyond the jellium boundary. Also, a
there is density peak around $x=2$ which its amplitude increases with increase
in the normalized electron temperature with the maximum value unchanged.
Figure 4(b) shows the corresponding electrostatic energy indicating its
maximum increase with the increase in the temperature. Figures 4(c) and 4(d)
show the electron density and electrostatic energy profiles for different
values of normalized chemical potential and fixed temperature. In Fig. 4(c) it
is observed that the increase of chemical potential leads to sharp increase in
electron density peak with its maximum value moved towards boundary. It is
also remarked that Friedel-like oscillation appear for much higher chemical
potential value characteristics of metallic density regime. The corresponding
electrostatic energy profile is shown in Fig. 4(d) which reveals that increase
in chemical potential leads to formation of oscillatory energy minimum in
metallic chemical potential and temperature regime. The appearance of
potential minimum beyond the physical surface of perfectly conducting metals
is fully consistent with well-known Casimir-Polder effect niko1 ; niko2 ;
bressi . It can also describe the nanometallic particle agglomeration in
quantum contact stefan .

Figure 5: 5(a) Variations in electrostatic energy profile in metallic chemical
potential regime with the change in the normalized chemical potential of the
electron gas. 5(b) Variations in electrostatic energy profile in metallic
chemical potential regime with the change in the normalized electron
temperature. 5(c) Variation of the electron number density with normalized
chemical potential for different values of electron number density. 5(d)
Variation of the normalized electron temperature with normalized chemical
potential for different values of normalized electron temperature. The
increase in the thickness of curves in each plot is meant to represent an
increase in the varied parameter above each panel.

Figure 5 shows the variations in electron density and electrostatic potential
profile with change in normalized temperature and chemical potential of the
electron gas in dense metallic regime. Figure 5(a) reveals that increase of
the chemical potential for fixed value of the electron temperature leads to
sharp increase of the Friedel oscillation amplitude. Giant Friedel
oscillations has been reported to exist in metallic surfaces such as beryllium
spr . On the other hand, Fig. 5(b) indicates that the increase of electron gas
temperature in fixed chemical potential leads to increase in depth of the
attractive potential energy, moving the potential minimum closer to the
boundary. Based on our model, beryllium with a Fermi energy of $E_{F}\simeq
14.3$eV, as one of densest metallic elements, should exhibit a significant
Casimir-Polder and nanoparticle agglomeration effects. Figures 5(c) and 5(d)
show the variations of $\theta$ and $\mu$ with respect to the electron number
density of electron gas. In terms of normalized parameters $\theta$ and
$\sigma$ the electron number density can be written as
${n_{0}}(\sigma,\theta)={n_{p}}{\theta^{6}}{\rm{L}}{{\rm{i}}_{3/2}}{\left[{-\exp(2{\sigma}/\theta)}\right]^{4}}$
with $n_{p}=16e^{6}m_{e}^{3}/(\pi^{3}\hbar^{6})\simeq 5.66\times 10^{19}$cm-3
being a characteristic number-density. According to Fig. 5(c), the number
density in $\sigma>1$ regime becomes almost independent of the value of
$\theta$. This is reason for why the electron temperature plays the
insignificant role in fully degenerate electron gas. It is therefore concluded
that the oscillatory attractive potential energy may be present for electron
densities of strongly doped N-type semiconductors. It is further confirmed
from by Fig. 5(d) that for higher values of $\theta$, such as in the warm
dense matter regime, the Friedel effect should also be present.

Figure 6: 6(a) Variation in matter-wave dispersion of plasmon quasiparticle
with changes in normalized chemical potential due to kinetic correction. 6(b)
Variation in matter-wave dispersion of plasmon quasiparticle with changes in
the normalized electron temperature due to kinetic correction. 6(c) Variation
in the quasiparticle desity of states (DoS) with chemical potential and
temperature. 6(d) Variation in the quasiparticle desity of states (DoS) with
quasiparticle energy. The increase in the thickness of curves in each plot is
meant to represent an increase in the varied parameter above each panel.

Figure 6(a) shows variation in the energy dispersion of quasiparticle
excitations with changes in chemical potential for fixed electron temperature.
The plasmon quasiparticle conduction band minimum is at $E_{m}=2\sqrt{\zeta}$
and $k_{m}=\zeta^{-1/4}$. Therefore, it is seen that the variations come from
the kinetic correction parameter, $\zeta(\sigma,\theta)$. It is revealed that
the kinetic correction mostly affect the low phase speed particle-like branch
of the dispersion. With increase in the chemical potential this branch moves
to energy values. On the other hand, for fixed chemical potential, increase in
the normalized electron temperature moves the low phase-speed branch to higher
energy values. The curvature of the quasiparticle conduction band is related
to the electron effective mass and consequently their mobility in collective
motion. Figure 6(c) shows the variations in quasiparticle DoS with chemical
potential and temperature of electron gas for given energy eigenvalue. It is
remarked that the largest variation of DoS with the temeparture corresponds to
the metallic regime $\sigma\simeq 1$. As the chemical potential increases, the
DoS also increase for fixed temperature and reaches a maximum value for higher
$\sigma$ values. However, increase of the temperature always lowers the
quasiparticle DoS. Figure 6(d) shows the variations in DoS with energy
eigenvalue. It is remarked that with increase in $E$ the DoS sharply
decreases. This is contrary to the case of free electrons in which we have the
relation $D_{f}(E)\propto\sqrt{E}$. It is also shown that increase in chemical
potential significantly enhances the DoS at lower quasiparticle energies.

##  VII Conclusion

In this research we used the pseudoforce model obtained from the kinetic-
corrected quantum hydrodynamic theory in order to investigate the half-space
fermion quasiparticle excitations in pure states. By using the appropriate
boundary conditions and separation of variables the quantized energy
eigenvalues of the pure states are obtained and appropriate mixed states were
developed for thermal equilibrium state. The thermodynamic parameters such as
the local electron density and electrostatic potential energy are calculated
from the mixed quasiarticle states based on the density of states (DoS) and
quasiparticle occupation probabilities. It was found that electron density
packet develops in front of the free electron boundary profile of which is
strongly dependent of the electron chemical potential and temperature. It was
also revealed that for electron density and temperature range relevant to
metallic regime Friedel like oscillations take place in the surface density
distribution. Furthermore, attractive potential energy minimum appears in
front of metallic surfaces which is consistent with the Casimir-Polder effect.
The effect of quantum kinetic correction on the energy dispersion and DoS of
fermion quasiparticle was also investigated in detail.

## References

  * (1) G. Manfredi, Phys. Plasmas 25, 031701(2018); https://doi.org/10.1063/1.5026653
  * (2) S. A. Maier, Plasmonics: Fundamentals and Aplications, Springer Science Business Media LLC (2007).
  * (3) H. Haug and S. W. Koch, ”Quantum theory of the optical and electronic properties of semiconductors”, World Scientific, 2004,
  * (4) C. Gardner, SIAM, J. Appl. Math. 54 409(1994).
  * (5) F. Wooten, Optical Properties of Solids, Academic (1972).
  * (6) P.A. Markovich, C.A. Ringhofer, C. Schmeister, Semiconductor Equations, Springer, Berlin, 1990.
  * (7) Manvir S. Kushwaha, Surf. Sci. Rep. 41, 1 (2001); doi.org/10.1016/S0167-5729(00)00007-8
  * (8) N. A. Krall and A. W. Trivelpeice, ”Principles of Plasma Physics”, (San francisco Press, San francisco 1986).
  * (9) Stenflo, L. 1976 Phys. Scripta 14, 320.
  * (10) S. Ummethala, T. Harter, K. Koehnle, Z. Li, S. Muehlbrandt, Y. Kutuvantavida, J. Kemal, P. Marin-Palomo, J. Schaefer, A. Tessmann, et al., Nat. Photonics 13, 519(2019); doi.org/10.1038/s41566-019-0475-6
  * (11) C. Calvero, Nat. Photonic 8, 95(2014).
  * (12) Jacob B. Khurgin, Nanophotonics 9(2), 453(2020).
  * (13) H. A. Atwater, Sci. Am. 296, 56(2007).
  * (14) H. A. Atwater and A. Polman, Nat. Mater. 9, 205(2010).
  * (15) Y. Tian and T. Tatsuma, J. Am. Chem. Soc. 127, 7632(2005).
  * (16) Vincenzo Amendola and Moreno Meneghetti, J. Phys. Chem. C 113, 4277(2009).
  * (17) C. Hu, Modern Semiconductor Devices for Integrated Circuits (Prentice Hall, Upper Saddle River, New Jersey, 2010) 1st ed.
  * (18) K. Seeger, Semiconductor Physics (Springer, Berlin, 2004) 9th ed.
  * (19) S. Glenzer, O. L. Landen, P. Neumayer, R. W. Lee, K. Widmann, S. W. Pollaine, R. J. Wallace, G. Gregori, A. Höll, T. Bornath, R. Thiele, V. Schwarz, W.-D. Kräft, and R. Redmer, Phys. Rev. Lett. 98, 065002 (2007).
  * (20) M. Koenig, A. Benuzzi-Mounaix, A. Ravasio, T. Vinci1, N. Ozaki, S. Lepape, D. Batani, G. Huser, T. Hall, D. Hicks, A. MacKinnon, P. Patel, H. S. Park, T. Boehly, M. Borghesi, S. Kar and L. Romagnani, Plasma Phys. Control. Fusion 47, B441 (2005).
  * (21) M. Bonitz, D. Semkat, A. Filinov, V. Golubnychyi, D. Kremp, D. O. Gericke, M. S. Murillo, V. Filinov, V. Fortov, W. Hoyer, J. Phys. A, 36 5921(2003).
  * (22) James P. Mithen, J$\acute{e}$r$\hat{o}$me Daligault, and Gianluca Gregori, Phys. Rev. E 82 015401(R) (2012).
  * (23) S. Chandrasekhar, ”An Introduction to the Study of Stellar Structure”, Chicago, Ill. (The University of Chicago press), (1939), p.392.
  * (24) C. Kittel, Introduction to Solid State Physics, (John Wiely and Sons, New York, 1996), 7th ed.
  * (25) N. W. Ashcroft and N. D. Mermin, Solid state physics (Saunders College Publishing, Orlando, 1976).
  * (26) D. Bohm and D. Pines, Phys. Rev. 92 609(1953).
  * (27) D. Pines, Phys. Rev. 92 609(1953).
  * (28) Y. Klimontovich and V. P. Silin, in Plasma Physics, edited by J. E. Drummond (McGraw-Hill, New York, 1961).
  * (29) S. Ichimaru, Rev. Mod. Phys. 54, 1017 (1982).
  * (30) P. K. Shukla and B. Eliasson, Phys. Rev. Lett. 99, 096401(2007).
  * (31) L Stenflo Phys. Scr. T50 15(1994).
  * (32) P. K. Shukla, B. Eliasson, and L. Stenflo Phys. Rev. E 86, 016403(2012).
  * (33) F. Haas, Quantum Plasmas: An Hydrodynamic Approach (Springer, New York, 2011).
  * (34) G. Brodin and M. Marklund, New J. Phys. 9, 277(2007).
  * (35) M. Marklund and G. Brodin, Phys. Rev. Lett. 98, 025001(2007).
  * (36) N. Crouseilles, P. A. Hervieux, and G. Manfredi, Phys. Rev. B 78, 155412 (2008).
  * (37) Z. Moldabekov, Tim Schoof, Patrick Ludwig, Michael Bonitz, and Tlekkabul Ramazanov, Phys. Plasmas, 22, 102104(2015); doi.org/10.1063/1.4932051
  * (38) L. Stanton and M. S. Murillo, Phys. Rev. E 91, 033104(2015).
  * (39) F. Haas, G. Manfredi, P. K. Shukla, and P.-A. Hervieux, Phys. Rev. B, 80, 073301 (2009).
  * (40) G. Manfredi, “How to model quantum plasmas,” Fields Inst. Commun. 46, 263–287 (2005); in Proceedings of the Workshop on Kinetic Theory (The Fields Institute, Toronto, Canada 2004): http://arxiv.org/abs/quant–ph/0505004.
  * (41) F. Haas, Quantum Plasmas: An Hydrodynamic Approach (Springer, New York, 2011).
  * (42) M. Bonitz, Zh. A. Moldabekov, and T. S. Ramazanov, Phys. Plasmas, 26, 090601 (2019); doi: 10.1063/1.5097885
  * (43) G. Manfredi and F. Haas, Phys. Rev. B 64, 075316 (2001);
  * (44) E. Madelung, Quantentheorie in hydrodynamischer Form. Zeitschr. f. Phys. 40, 322–326 (1926).
  * (45) G. Manfredi, 86(2), 825860201(2020); doi:10.1017/S0022377820000240
  * (46) J. Hurst, K. L. Simon, P. A. Hervieux, G. Manfredi and F. Haas, Phys. Rev. B 93, 205402(2016).
  * (47) M. Akbari-Moghanjoughi, Phys. Plasmas, 26, 012104 (2019); doi.org/10.1063/1.5078740
  * (48) M. Akbari-Moghanjoughi, Phys. Plasmas 22, 022103 (2015); ibid. 22, 039904 (E) (2015).
  * (49) Zh. Moldabekov, T. Schoof, P. Ludwig, M. Bonitz, and T. Ramazanov, Phys. Plasmas 22, 102104 (2015).
  * (50) D. Michta, F. Graziani, and M. Bonitz, Contrib. Plasma Phys. 55, 437 (2015).
  * (51) D. Bohm, Phys. Rev. 85, 166 (1952).
  * (52) S. A. Dyachkov, P. R. Levashov, Phys. Plasmas 21 052702 (2014).
  * (53) S. A. Dyachkov, P. R. Levashov, D. V. Minakov, Phys. Plasmas 23 112705 (2016).
  * (54) Setsuo ICHIMARU, Hiroshi IYETOMI and Shigenori TANAKA, Phys. Rep. 205, 91(1987).
  * (55) M. Bonitz, Zh. A. Moldabekov, and T. S. Ramazanov, Phys. Plasmas 26, 090601 (2019); doi.org/10.1063/1.5097885
  * (56) H. Nikolic, Found. Phys. Lett. 17 363(2004).
  * (57) Cristian Ciracì, Radoslaw Jurga, Muhammad Khalid and Fabio Della Sala, Nanophotonics 8(10), 821(2019); doi.org/10.1515/nanoph-2019-0199
  * (58) Mario Zapata Herrera, Javier Aizpurua, Andrey K. Kazansky, and Andrei G. Borisov, Langmuir 32, 2829(2016); doi:10.1021/acs.langmuir.6b00112
  * (59) Hrvoje Nikolic, Annals of Physics 383, 181(2017).
  * (60) Hrvoje Nikolic, Physics Letters B 761, 197(2016).
  * (61) G. Bressi, G. Carugno, R. Onofrio, and G. Ruoso, Phys. Rev. Lett. 88, 041804(2002).
  * (62) Aleksandar D. Rakic, Aleksandra B. Djuris, Jovan M. Elazar, and Marian L. Majewski, Apll. Optics, 37(22) 5271(1998).
  * (63) Y. Ohko, Yoshihisa Ohko, Tetsu Tatsuma, Tsuyoshi Fujii, Kenji Naoi, Nat. Mater. 2, 29(2003).
  * (64) C. Sonnichsen, T. Franzl, T. Wilk, G. von Plessen, J. Feldmann, O. Wilson, and P. Mulvaney, Phys. Rev. Lett. 88, 077402(2002).
  * (65) M. Akbari-Moghanjoughi, arXiv:2009.08098 [physics.plasm-ph] (under review).
  * (66) Tigran V. Shahbazyan, Phys. Rev. B 94, 235431(2016).
  * (67) Bengt Eliasson, and Chuan Sheng Liu, Phys. Plasmas 25, 012105(2018).
  * (68) F. Haas and S. Mahmood, Phys. Rev. E 92, 053112 (2015).
  * (69) M. Akbari-Moghanjoughi, Phys. Plasmas, 26, 112102 (2019); doi.org/10.1063/1.5123621
  * (70) M. Akbari-Moghanjoughi, Alireza Abdikian, and Arash Phirouznia1, Phys. Plasmas 27, 042107 (2020); doi.org/10.1063/5.0004857
  * (71) Stefan Christian Endres, Lucio Colombi Ciacchi and Lutz Madler, J. Aerosol Science 153, (105719)2021; doi.org/10.1016/j.jaerosci.2020.105719
  * (72) P. T. Sprunger, L. Petersen, E. W. Plummer, E. Lægsgaard, F. Besenbache, Science 275 1764(1997); doi:10.1126/science.275.5307.1764