# Uncertainties in the Transport Properties of Helium Gas at Cryogenic
Temperatures Determined Using Molecular Dynamics Simulation

Kaanapuli Ramkumar  Swati Swagatika Mishra  Sudeep Bhattacharjee

###### Abstract

The transport properties of 4He are determined in the temperature range of 10
K to 150 K and pressure range of 0.101 MPa (1 atm) to 0.210 MPa, in classical
and quantum frameworks, employing numerical and molecular dynamics simulation
(MDS) techniques. Although classical MDS is efficient and comprehensive in
understanding the transport properties, the increasing effect of quantum
diffraction in the temperature range of 10 - 40 K is not accounted for in the
MDS. The statistical error associated with the Green Kubo (GK) approach also
adds to the uncertainty in the simulated transport coefficients. The article
presents a detailed analysis of the uncertainties in the classical MDS to
study the transport phenomena of 4He at low temperatures. In particular, it
addresses two questions, how the quantum effects manifest at different
temperature and pressure regimes and how the statistical error in the
autocorrelation functions of the GK approach can be minimised.

###### keywords:

Molecular dynamics simulation , Transport properties , Green-Kubo method ,
Quantum effects , Uncertainties

[inst]organization=Department of Physics, Indian Institute of Technology
Kanpur,city=Kanpur, postcode=208016, state=Uttar Pradesh, country=India

##  1 Introduction

knowledge of transport properties of simple gases to complex fluidic systems
is important in various interdisciplinary fields such as fuel combustion
technologies [1], nuclear reactor engineering [2], extraterrestrial
explorations [3] and nuclear fusion [4] including plasmas technologies [5].
Transport coefficients such as diffusion, thermal conductivity, and shear
viscosity that quantify the transfer of energy, momentum, and mass in a fluid
act as a bridge connecting the microscopic molecules in a fluid with the
macroscopic system as a whole. The interaction and the properties of the
individual molecules in a fluid determine the transport coefficients, which
can be measured at the macroscopic level. The design of thermal protection
systems for planetary exploration probes [3, 6], modeling of atmospheres of
other planets [7, 5], designing inertial confinement fusion reactors [4, 5]
and the modeling of plasma processes used in gasification of biomass [8]
require extensive data on the transport properties of various species.
Furthermore, the transport properties of neutral gases become very important
for the understanding of the dynamics of charged species in a weakly ionized
plasma medium, especially at cryogenic temperatures (4 K to 10 K) and at
atmospheric pressure [9, 10].

In recent years, transport properties of different systems have been studied
in detail employing molecular dynamics simulations (MDS). The advantage of
using simulation is that it has flexibility and allows one to manipulate any
parameter of a system, thereby permitting one to understand the effects of
different parameters in greater detail without the limitations of conventional
experimental investigations. With the help of MDS, Schaefer, and Kohns studied
the ion clustering in concentrated electrolytes [11] and Takahashi et al.
studied the quantum effect of liquid hydrogen on bubble nucleation [12].
Yiannourakou et al. performed an MDS study to optimize the design and
performance of acid gas removal units [13]. González-Salgado et al. studied
the anomalous behavior of water at 4oC with the help of MDS [14]. Cunha and
Robbins studied the pressure-viscosity relation of 2,2,4-trimethylhexane using
MDS [15]. Recently, the temperature guided behavioral transitions in confined
helium at cryogenic temperatures (30 K to 150 K) were looked at by our group
[16]. In the current work, we extend the investigation to even lower
temperatures where the quantum effects would be important.

All these studies are a testimony to the importance of MDS in understanding
physical phenomena. The most prominently used approach to calculating
transport properties using MDS is the Green Kubo (GK) approach. However, the
GK approach is known to have a significant statistical error associated with
it [17, 18]. Furthermore, at low temperatures and high pressure regimes, the
gas molecules tend to show their quantum nature and the classical MDS are
unable to provide accurate results, leading to the divergence in the
classically calculated transport coefficients from the experimental values.
Therefore, the quantum nature of the molecules is also a source of error while
using the classical MDS [19, 20]. For the efficient usage of any simulation
result in practical applications, one has to be aware of the different
uncertainties involved with the simulation. The purpose of the present article
is to understand the uncertainties associated with the calculation of
transport properties in various temperature and pressure regimes. The primary
focus of the article is on the uncertainties associated with the GK approach
and the manifestation of quantum effects at low temperature and high pressure
regimes. For this purpose, the Lennard-Jones (LJ) potential is used to
describe the intermolecular interactions and study the transport coefficients
such as viscosity ($\eta$), thermal conductivity ($\kappa$), and diffusivity
($\mathcal{D}$) of 4He gas.

Figure 1: Calculation of transport coefficients within different frameworks.

This article is divided into five sections. Section II briefly describes the
theoretical concepts of calculating transport coefficients and quantum effects
manifesting at low temperatures. Section III contains the details of the MDS
simulations and numerical quantum calculations. Section IV discusses the
results of the simulation and numerical calculations, and finally, the article
is concluded in Section V.

##  2 Theoretical Background

###  2.1 Transport Properties

The article studies transport properties of 4He gas in a temperature regime
where the quantum effects starts to emerge (around 30 K). Hence, it is
necessary to study them in both the classical and the quantum formalisms, as
the investigated temperature range in this work is from 10 K to 150 K. Within
the classical formalism, the transport properties could be calculated in two
ways: (i) by simulation methodologies such as MDS, using the Green-Kubo method
[21, 22] and (ii) through numerical methods employing Chapman-Cowling
approximations [23]. For simple systems such as noble gases where binary
interactions are predominant, one could use numerical methods, which involves
calculation of scattering angles in binary collisions to be used in Chapman-
Enskog theory and Chapman-Cowling approximations [23]. In the quantum
formalism, we use the same numerical method used in the classical method, but
the Schrödinger’s equation is solved to calculate the phase shifts of
scattering processes, which play the quantum analogous role of scattering
angle computed in the classical formalism. The flow chart in Fig. 1 depicts
how these different methodologies are employed in the calculation of the
transport coefficients and the following subsections briefly summarises them.

####  2.1.1 Green Kubo Theory

The GK relations [21, 22], uses the autocorrelation functions (ACF) of the
flux corresponding to the transport properties to evaluate the transport
coefficients. The viscosity ($\eta$) is calculated from the ACF of the off-
diagonal pressure tensor ($P_{ij}(t)$), given by [24],

$\eta=\frac{V}{k_{B}T}\int_{0}^{\infty}\langle
P_{ij}(t_{0})P_{ij}(t+t_{0})\rangle_{t_{0}}dt,$ (1)

and the thermal conductivity ($\kappa$) from the heat flux ($J(t)$) ACF, given
by [25],

$\kappa=\frac{V}{3k_{B}T}\int_{0}^{\infty}\langle J(t_{0})\cdot
J(t+t_{0})\rangle_{t_{0}}dt,$ (2)

and the diffusion coefficient ($\mathcal{D}$) from the velocity ($v(t)$) ACF,
given by [26],

$\mathcal{D}=\frac{1}{3}\int_{0}^{\infty}\langle v(t_{0})\cdot
v(t+t_{0})\rangle_{t_{0}}dt,$ (3)

where $t_{0}$ and $t$ are any arbitrary initial and instantaneous times,
respectively.

####  2.1.2 Chapman-Cowling Approximation

To obtain the transport properties numerically from the Chapman-Cowling
approximation, the scattering angles for binary collisions between the gas
molecules are used to evaluate the collision cross-sections given by [23]

$Q^{(r)}=2\pi\int_{0}^{\infty}\left(1-cos^{r}\zeta\right)\ b\ db,$ (4)

where $\zeta$ is the scattering angle and $b$ is the scattering parameter and
$r$ takes the values 1,2. The collision cross-section thus evaluated is used
to calculate the collision integrals ($\Omega$), given by [23],

$\Omega^{(r,s)}=\left(\frac{k_{B}T}{\pi
m}\right)^{1/2}\int_{0}^{\infty}Q^{(r)}e^{-g^{2}}g^{2s+3}dg,$ (5)

where $s$ takes the values 1,2 and $g$ is the non-dimensional relative
velocity of the scattering particles, given by $g=\hbar\ k/\sqrt{2k_{B}T\mu}$,
where $\mu$ is the reduced mass and $k$ is the wave number of the particles.
The Boltzmann’s transport equations are solved using the Chapman-Enskog theory
[27], and the transport coefficients are obtained using the Chapman-Cowling
approximation, to the first order, as [27],

$\displaystyle\eta=\frac{5k_{B}T}{8\Omega^{(2,2)}},$ (6a)
$\displaystyle\kappa=\frac{5k_{B}T}{8\Omega^{(2,2)}}\frac{15k_{B}}{4m},$ (6b)
$\displaystyle\mathcal{D}=\frac{3(k_{B}T)^{2}}{8Pm\Omega^{(1,1)}}.$ (6c)

Here $m$ is the mass of the interacting particles and $P$ corresponds to the
pressure of the system, in SI units, at which the transport coefficients are
calculated.

In quantum mechanics, due to the uncertainty principle, the relation between
the scattering angle ($\chi$) and the impact parameter ($b$) cannot be defined
to calculate the collision cross-sections (Eq. 4). Hence to do a quantum
mechanical calculation, instead of scattering angles, the phase shifts
($\delta_{l}$) in the asymptotic region of scattering state wave functions are
used to obtain the quantum mechanical collision cross-sections, given by [28,
23, 29],

$Q^{(1)}=\frac{4\pi}{k^{2}}\sum_{l=0}^{\infty}(l+1)\>sin^{2}(\delta_{l+1}-\delta_{l}),$
(7)

and

$Q^{(2)}=\frac{2\pi}{k^{2}}\sum_{l=0}^{\infty}\frac{(l+1)(l+2)}{l+3/2}\>sin^{2}(\delta_{l+2}-\delta_{l}),$
(8)

which can be used in the collision integral equation (Eq. 5) to evaluate the
transport coefficients from the Chapman-Cowling expressions (Eqs. 6). Here,
$k$ is the wave number, and $l$ is the angular momentum quantum number of the
state considered. To include the effects of quantum statistics, the summations
in the above formulas are carried out only over the even $l$ values for the
Bose-Einstein case and over the odd $l$ values for the Fermi-Dirac case due to
symmetry and in the case of collision between indistinguishable particles, the
formulas are multiplied by a factor of 2 [23, 28]. A detailed description of
the quantum effects is provided in the next subsection.

###  2.2 Quantum Effects

Due to the wave-particle duality of quantum particles, at low temperatures,
when the thermal de Broglie wavelength ($\lambda_{th}$), given by the
expression
$\lambda_{th}=\sqrt{2\pi\hbar^{2}/(k_{B}mT)}$, where $k_{B}$ is the
Boltzmann’s constant, $m$ is the mass of the molecule, $T$ is the temperature,
and $\hbar$ is the reduced Planck’s constant, is comparable to the particle’s
dimensions, diffraction effects become dominant [29]. Furthermore, at
cryogenic temperatures, the particles with spin half follow Fermi-Dirac
statistics, and particles with integer spin follow Bose-Einstein statistics
and not the classical Boltzmann distribution. When $\lambda_{th}$ is
comparable to the inter-particle distances, the effects due to the quantum
statistics (symmetry effects) become dominant [29].

For 4He atoms, the Van der Waals radius ($r_{v}$) is 1.4 Å. This corresponds
to $\lambda_{th}$ of 4He at $\sim$ 45 K. Fig. 2 shows how the mean inter-
particle distances ($d$) in 4He gas, obtained from the number density ($n$)
data at atmospheric pressure [30], from the relation $d=n^{-1/3}$, approaches
$\lambda_{th}$ at a much lower temperature, less than 4.222 K, which is the
boiling point of 4He gas at atmospheric pressure [30]. Hence, while
considering quantum effects, quantum statistics effects may be neglected for
the temperature range of 10 - 150 K.

Figure 2: Comparing Thermal de Broglie wavelength ($\lambda_{th}$) (solid
line) with mean inter-particle distance ($d$) (dotted line). The Van der Waals
radius ($r_{v}$) is marked by the dashed line for reference.

It should be noted that under the quantum framework, the calculation of the
diffusion coefficient needs special considerations. Classically, one can track
the motion of individual particles to find the self-diffusion coefficient but
quantum mechanically, when identical particles interact, they lose their
identity after the interaction (indistinguishability) [29, 27]. This causes
the quantum mechanically obtained values for self-diffusion scattering cross-
section to be almost twice that of classically obtained cross-section and
surprisingly, the classical values are closer to the experimental values than
the quantum mechanical values, even at the cryogenic temperatures. This is
because experimentally it becomes impossible to track the motion of individual
identical particles and hence, some property of the particles, which do not
affect the collision process, is used to differentiate the particles in the
system [31]. For instance, in the case of 3He self-diffusion experiments, the
molecules are categorised based on their spin state (+½ and -½) [32] as spin
does not affect collision interactions. Similarly, different isotopes with
mass ratios very close to unity are also used in calculating self-diffusion
coefficients experimentally [33]. By doing this, the interacting particles are
made distinguishable and the experimentally measured self-diffusion
coefficients are not that of true self-diffusion coefficients of identical
particles. Theoretically, one can calculate these values by assuming a
distinguishable scattering process where the self-diffusion coefficient is
calculated as a limiting case of binary diffusion [29].

##  3 Computational Details

###  3.1 Interaction Potential

In all the calculations done in this article, LJ potential is considered to
depict the interaction between the 4He atoms [34], which is given by

$V_{LJ}(r)=4\epsilon\left[\left(\frac{r}{\sigma}\right)^{-12}-\left(\frac{r}{\sigma}\right)^{-6}\right],$
(9)

with the energy parameter $\epsilon=$9.44$$ meV and the length parameter
$\sigma=$2.64$$ $\mathrm{\SIUnitSymbolAngstrom}$ [35].

###  3.2 Classical Molecular Dynamics Simulation

The Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [36] is
used to perform the MDS. While the GK relations (Eq. 1,2,3) are used to
calculate the transport properties, the diffusion coefficients are also
calculated from the asymptotic time slope of mean square displacement (MSD) of
the particles in the system, using the Einstein’s relation given by,
$\mathcal{D}=\lim_{t\to\infty}\frac{1}{6t}MSD$ [26]. This has been done to
show that the uncertainties due to quantum effects are innate to classical
formalism and does not depend on the exact method used to calculate the
transport coefficient values.

####  3.2.1 Simulation Parameters

The simulation system size was set at
$660$\mathrm{\SIUnitSymbolAngstrom}$\times
660$\mathrm{\SIUnitSymbolAngstrom}$\times 660$\mathrm{\SIUnitSymbolAngstrom}$$
and is kept constant throughout all the simulations performed in this work.
The LJ potential used is truncated at a cut-off distance of 10
$\mathrm{\SIUnitSymbolAngstrom}$ for all simulations and this value is based
on the prescribed cutoff length for LJ interactions, which is
$r_{c}>2.5\sigma$ [34]. Beyond this cut-off length, the interaction potential
is small enough (< 0.02 times $\epsilon$) and can be neglected. The time
interval between the two simulation steps is 1 fs (step size) and the initial
system equilibration is performed for 5 ns as an NPT ensemble and the
thermalization of the system is performed under NVE ensemble setting for 1 ns.
The periodic boundary conditions have been applied to the simulation system.

A number of simulations were performed to calculate the transport properties
at various temperatures in the range of 10 K to 150 K for 0.101 MPa pressure
(1 atm) and various other pressures in the range of 0.14 MPa to 0.21 MPa for
temperatures 30 K and 100 K. Tables 1 and 2 show the number of atoms simulated
in each simulation.

Table 1: Number of atoms present in the system for simulation at various temperatures at 0.101 MPa. Temperature (K) |  No. of atoms
---|---
10 |  211090
30 |  70364
40 |  52773
50 |  42218
70 |  30156
100 |  21109
150 |  14073
Table 2: Number of atoms present in the system for simulations at various pressures at 30 K and 100 K. Pressure (MPa) | No. of atoms
---|---
|  30 K |  100 K
0.14 |  97221 |  29166
0.16 |  111110 |  33333
0.18 |  124998 |  37499
0.21 |  145831 |  43749

####  3.2.2 Green Kubo ACF Integration Details

The primary source of uncertainty in the GK approach is the time integral over
the ACF in Eqs. 1,2 and 3 [17]. Considering the ACF for longer time periods
above the correlation length leads to increased error due to the accumulation
of fluctuations in the tail of the ACF. Prematurely truncating the ACF before
it sufficiently decays also leads to error in the calculated transport
properties. Various methods based on exponential fitting or appropriately
truncating the ACF have been proposed to minimise this integration error [24,
37]. In this article, to minimise the error due to random fluctuations in the
ACF, the ACF is calculated as a cumulative average of a few thousand
simulation runs. Also, for simple fluids in the gaseous state, the decay in
the ACF is exponential in nature [18]. Hence, the ACF is fit to an exponential
function, $a\ e^{-t/\tau}$, where $\tau$ gives the correlation length. This
removes the dependence of the transport properties on the ACF integration
cutoff and any residual white noise, after the cumulative averaging, on the
ACF tail.

###  3.3 Numerical Quantum Calculation

One of the disadvantages of the purely classical MDS is that at low
temperatures, the QM effects are not taken into account. In the case of 4He
gas, the QM effects start to dominate at relatively higher temperatures (at
around 40 K) than for any other gases. Using the quantum mechanical theory,
the transport properties can be obtained by solving for the phase shifts in
the asymptotic region of the scattering state wave function, which then can be
used to calculate the collision cross-section in the Chapman-Cowling relations
[27].

To obtain the phase shifts, the Schrödinger’s equation, given by

$\frac{\partial^{2}\mathcal{U}(r)}{\partial r^{2}}=f(r)\ \mathcal{U}(r),$ (10)

is solved numerically using the Numerov’s integration method [38], where the
wave function is discretised as,

$\mathcal{U}_{i+1}=\frac{\mathcal{U}_{i-1}\ (12-d^{2}f_{i-1})-2\
\mathcal{U}_{i}\ (5d^{2}f_{i}+12)}{d^{2}f_{i+1}-12}.$ (11)

Here, $d$ is the step size considered for discretizing the wave function, $l$
is the angular momentum state considered, $\mu$ is the reduced mass of two
colliding 4He atoms and for the case of LJ potential, the function $f(r)$ is
given by,

$f(r)=\frac{2\mu}{\hbar^{2}}\left(l(l+1)\frac{\hbar^{2}}{2\mu
r^{2}}+V_{LJ}(r)-E\right),$ (12)

where $E$ is the scattering state energy eigenvalue. In our calculations, the
wave functions are discretized with a step size of 0.025 a.u. and are
evaluated up to a distance of 1500 a.u.

The analytical expression for the asymptotic scattering state wave function,
for a given angular momentum state $l$, is given by [39],

$\mathcal{U}(r)=A\>kr\ (j_{l}(kr)\ \cos\ \delta_{l}-n_{l}(kr)\ \sin\
\delta_{l}),$ (13)

where $A$ is a normalization factor, $\delta_{l}$ is the phase shift
corresponding to a given angular momentum state $l$, $j_{l}$ and $n_{l}$ are
the spherical Bessel and Neumann functions, and $\mathcal{U}(r)=R(r)\ r$, with
$R(r)$ being the radial part of the scattering wave function. The phase shift
($\delta_{l})$ can be evaluated using,

$\Gamma=\frac{r_{2}\
\mathcal{U}(r_{1})}{r_{1}\mathcal{U}(r_{2})}=\frac{j_{l}(kr_{1})\ \cos\
\delta_{l}-n_{l}(kr_{1})\ \sin\ \delta_{l}}{j_{l}(kr_{2})\ \cos\
\delta_{l}-n_{l}(kr_{2})\ \sin\ \delta_{l}},$ (14)

and

$\delta_{l}=\frac{j_{l}(kr_{1})-\Gamma\ j_{l}(kr_{2})}{n_{l}(kr_{1})-\Gamma\
n_{l}(kr_{2})}.$ (15)

To obtain the quantum scattering cross-section, the phase shift values must be
calculated for all even angular momentum states at different energies from 0
to $\infty$, as the quantum cross-sections have summation up to infinite
angular momentum state (Eqs. 7 and 8) and the collision integral has integral
over wave number from zero to infinity (Eq. 5). This is approximated by
carrying out the summation over angular momentum states in Eqs. 7 and 8, and
integration over the wave vector, in Eq. 5, until the change in the cross-
section and collision integral is not more than a factor of $10^{-8}$ after
the inclusion of the phase shifts for higher angular momentum states and
energy values. The integration in collision integral (Eq. 5) is performed
using the trapezoidal method as it varies smoothly with wave number. For the
temperature region of interest in this article (10 K to 150 K), a maximum
angular momentum state of $l=22$ and a maximum wave vector of 9 a.u. is
considered. Here, it is assumed that only binary interactions take place
during a collision, i.e. only two particle collisions are considered [23, 27].

##  4 Results

The values of transport coefficients, namely viscosity ($\eta$), thermal
conductivity ($\kappa$), and diffusivity ($\mathcal{D}$) calculated using
different methods in the temperature range of 10 K to 150 K and pressure range
of 0.101 MPa to 0.210 MPa are provided in Tables 3, 4 and 5. The statistical
error associated with the GK approach in calculating viscosity and thermal
conductivity is mentioned along with their values. The diffusion coefficient,
which was obtained from the asymptotic slope of mean square distance, has
negligible uncertainty in its slope and hence is not mentioned. The quantum
mechanical calculations have no statistical error associated with them and can
be obtained to any degree of accuracy.

###  4.1 Fluctuations in ACF and statistical error associated with Green Kubo
approach

As mentioned in the previous sections, to minimize the error due to
fluctuations in the ACF, it is calculated as a cumulative average of the
previous runs. The random fluctuations in the ACFs are white noises which
cancel out when averaged over a large data sample. Here a single run refers to
the simulation duration over which one ACF is calculated. The duration of this
changes for each simulation depending on the temperature and pressure
conditions as the time the ACF takes to decay varies with these parameters.
Fig. 3 shows how the heat flux ACF converges to an exponentially decaying
function when cumulatively averaged over several runs. The inset in Fig. 3
shows the thermal conductivity calculated at the end of each run using direct
integration of the cumulatively averaged ACF.

Figure 3: Heat flux autocorrelation function (ACF) calculated in a cumulative
manner over several runs at 30 K and 0.101 MPa. The inset shows the thermal
conductivity value calculated at the end of each run.  Figure 4: Viscosity
calculated using direct integration of the autocorrelation and exponentially
fitting autocorrelation function (ACF) at 30 K and 0.210 MPa. The inset shows
the frequency distribution curve of values of viscosity calculated after 1500
runs.  Figure 5: Percentage deviation in the calculated transport properties,
from the quantum mechanically calculated values, for different pressures at,
(a) 30 K and (b) 100 K.

In principle, an average over infinite data sets is required to remove white
noise completely from the system. Hence to minimise any residual fluctuations
in the ACF and errors due to cut-offs in the ACF integration, cumulatively
averaged ACF is fit to an exponentially decaying function ($a\ e^{-t/\tau}$)
and the integral over ACF is analytically given by $a\tau$. Fig. 4, shows the
viscosity calculated at 30 K and 0.21 MPa by direct integration (red line) of
the cumulatively averaged ACF and by analytical integration of exponentially
fit cumulatively averaged ACF. It can be observed that a better convergence is
obtained in the calculated viscosity values. In the inset of Fig. 4, the
number of times a particular viscosity value appears when the simulation is
run several times ($\sim 1500\ runs$) is given in the form of a frequency
distribution plot. The red curve shows the frequency distribution for the
values by direct integration (trapezoidal method) of ACF and the blue curve
shows the distribution for the values obtained by integration of the
exponential fit. The width of the frequency distribution plot obtained for the
exponentially fit ACF is lower than that obtained by direct integration. This
shows a reduction in the standard deviation (statistical error) of the
obtained viscosity values. It can also be observed that the two curves peak at
different values. This shows that the exponential fit also accounts for the
data lost from the tail of the ACF if an arbitrary cutoff is used.

###  4.2 Uncertainty due to quantum effects

Fig. 5 (a) and (b) show the percentage deviation in the transport
coefficients, obtained through (I) MDS and (II) from the quantum mechanically
calculated values at different temperatures and pressures. Fig 6 shows how the
deviation in the transport coefficients increase with a decrease in
temperature. Below 40 K, at 0.101 MPa, the quantum mechanical effects are
considerable, causing deviation in the classically simulated values. In Fig. 5
(a), at a temperature of 30 K, the deviations increase with the increase in
pressure. This can be attributed to the fact that with an increase in
pressure, the mean interparticle distance decreases and hence approaches the
$\lambda_{th}$ at higher temperatures. In Fig. 5 (b), the deviations in the
transport coefficients remain roughly the same with pressure. This shows that
the pressure range considered is not high enough to make the quantum effects
prominent at 100 K. It is also evident that the uncertainties due to quantum
effects are present independent of the calculation procedure, GK or MSD, used
within the classical MDS.

Figure 6: Percentage deviation in the calculated transport properties, from
the quantum mechanically calculated values, at different temperatures at 0.101
MPa pressure.

##  5 Discussion

In this article, the transport properties such as viscosity, diffusion
coefficient, and thermal conductivity of helium gas are calculated in the
temperature range of 10 K to 150 K and the pressure range of 0.101 MPa to
0.210 MPa. Temperatures below the critical temperature (5.2 K), where the
effects due to quantum statistics are considerable, and pressures above the
critical pressure (0.22 MPa), where the 4He is in the supercritical fluid
state are not considered in this study.

The most prominent sources of uncertainty in the transport coefficient
calculations at cryogenic temperatures using MDS are the quantum effects and
the statistical errors associated with the GK approach. From the simulations
performed, it is observed that the uncertainties due to quantum effects
increase with a decrease in temperature and an increase in pressure due to the
mean interparticle distance approaching the thermal de Broglie wavelength.

The uncertainties from the GK approach can be minimised to some extent using
different integration methods for the ACF but the uncertainties due to the
quantum effects can’t be removed as long as one confines themselves to the
classical MDS. One can perform quantum MDS, which includes the quantum
effects, but quantum simulations are orders of magnitude more computationally
intensive than classical simulations. Another way around this problem is to
use existing quantum correction procedures for classical MDS, but they are
limited in their applicability. The well-known Feynman-Hibbs (FH) variational
approach, which accounts for the positional uncertainty of quantum particles
[40, 41], has been used successfully for many studies such as the adsorption
of hydrogen and deuterium [42], quantum effects in light and heavy water [43],
equation of state for solid neon [44] and transport properties of 4He at 80 K,
neon at 45 K and methane at 110 K [45]. This approach has its limitations.
Although it works well in moderate to high temperature regimes (> 80 K, for
the case of 4He), it breaks down in the temperature regime where the quantum
diffraction effects [45], explained in Sec. II B, are high (T < 40 K, for
4He). This is attributed to the FH method’s inability to correctly yield the
kinetic energy of the particles in the system due to the increased
delocalisation at lower temperatures [46]. A worthwhile future research
direction would be exploring different correction procedures to include
quantum effects in the classical simulation methods at cryogenic temperatures.

##  6 Acknowledgement

The authors sincerely thank the High-Performance Computing facility (HPC) at
IIT Kanpur, because of which this work has been made possible.

##  7 Author Declaration

###  7.1 Conflict of Interest

The authors declare no conflict of interest.

###  7.2 Data Availability

The data that supports the findings of this study are available within the
article.

Table 3: Viscosity ($\nu$) ($\mu Pa/s$) data obtained using numerical quantum calculations (QM), and molecular dynamics simulation (MDS) with the statistical error associated with the Green Kubo approach. Temperature (K) |  Pressure (MPa) | Viscosity ($\nu$)
---|---|---
|  |  QM |  MDS (GK)
10 |  0.101 |  $2.02$ |  $1.55$$\pm$0.01
30 |  |  $4.25$ |  $4.08$$\pm$0.01
40 |  |  $5.12$ |  $4.93$$\pm$0.02
50 |  |  $5.91$ |  $5.93$$\pm$0.03
70 |  |  $7.34$ |  $7.34$$\pm$0.02
100 |  |  $9.23$ |  $9.25$$\pm$0.03
150 |  |  $12.1$ |  $11.9$$\pm$0.1
30 |  0.140 |  $4.25$ |  $4.35$$\pm$0.03
|  0.160 |  |  $4.39$$\pm$0.01
|  0.180 |  |  $4.59$$\pm$0.02
|  0.210 |  |  $4.65$$\pm$0.03
100 |  0.140 |  $9.23$ |  $9.47$$\pm$0.04
|  0.160 |  |  $9.75$$\pm$0.07
|  0.180 |  |  $9.52$$\pm$0.05
|  0.210 |  |  $9.65$$\pm$0.04
Table 4: Thermal conductivity ($\kappa$) ($10^{-2}W/(K-m)$) data obtained using numerical quantum calculations (QM), and molecular dynamics simulation (MDS) with the statistical error associated with the Green Kubo approach. Temperature (K) |  Pressure (MPa) | Thermal conductivity ($\kappa$)
---|---|---
|  |  QM |  MDS (GK)
10 |  0.101 |  $1.57$ |  $1.02$$\pm$0.01
30 |  |  $3.31$ |  $3.14$$\pm$0.01
40 |  |  $3.99$ |  $3.84$$\pm$0.01
50 |  |  $4.61$ |  $4.73$$\pm$0.03
70 |  |  $5.72$ |  $5.69$$\pm$0.03
100 |  |  $7.19$ |  $7.07$$\pm$0.03
150 |  |  $9.40$ |  $9.76$$\pm$0.06
30 |  0.140 |  $3.31$ |  $3.34$$\pm$0.02
|  0.160 |  |  $3.32$$\pm$0.02
|  0.180 |  |  $3.50$$\pm$0.03
|  0.210 |  |  $3.57$$\pm$0.02
100 |  0.140 |  $7.19$ |  $7.5$$\pm$0.1
|  0.160 |  |  $7.50$$\pm$0.07
|  0.180 |  |  $7.4$$\pm$0.1
|  0.210 |  |  $7.54$$\pm$0.09
Table 5: Diffusion coefficients ($\mathcal{D}$) ($10^{-6}m^{2}/s$) data obtained using numerical quantum calculations (QM), and molecular dynamics simulation (MDS) with the statistical error associated with the Green Kubo approach. Temperature (K) |  Pressure (MPa) | Diffusivity ($\mathcal{D}$)
---|---|---
|  |  QM |  MSD |  MDS (GK)
10 |  0.101 |  $0.55$ |  $0.47$ |  $0.50$$\pm$0.001
30 |  |  $3.48$ |  $3.25$ |  $3.29$$\pm$0.005
40 |  |  $5.61$ |  $5.35$ |  $5.40$$\pm$0.003
50 |  |  $8.11$ |  $7.81$ |  $7.61$$\pm$0.006
70 |  |  $14.1$ |  $13.8$ |  $13.8$$\pm$0.01
100 |  |  $25.4$ |  $25.4$ |  $25.4$$\pm$0.02
150 |  |  $49.6$ |  $49.3$ |  $50.6$$\pm$0.1
30 |  0.14 |  $2.52$ |  $2.44$ |  $2.50$$\pm$0.002
|  0.16 |  $2.20$ |  $2.16$ |  $2.21$$\pm$0.002
|  0.18 |  $1.96$ |  $2.02$ |  $2.10$$\pm$0.004
|  0.21 |  $1.68$ |  $1.77$ |  $1.80$$\pm$0.001
100 |  0.14 |  $18.4$ |  $18.4$ |  $18.1$$\pm$0.02
|  0.16 |  $15.8$ |  $16.1$ |  $16.0$$\pm$0.02
|  0.18 |  $14.1$ |  $14.3$ |  $14.3$$\pm$0.02
|  0.21 |  $12.1$ |  $12.3$ |  $12.3$$\pm$0.01

## References

  * [1] H. M. Ismail, H. K. Ng, X. Cheng, S. Gan, T. Lucchini, G. D’Errico, Development of thermophysical and transport properties for the CFD simulations of in-cylinder biodiesel spray combustion, Energy Fuels 26 (8) (2012) 4857.
  * [2] A. Rodriguez, S. Lam, M. Hu, Thermodynamic and transport properties of lif and flibe molten salts with deep learning potentials, ACS Appl. Mater. Interfaces 13 (46) (2021) 55367.
  * [3] D. Bruno, C. Catalfamo, M. Capitelli, G. Colonna, O. De Pascale, P. Diomede, C. Gorse, A. Laricchiuta, S. Longo, D. Giordano, F. Pirani, Transport properties of high-temperature jupiter atmosphere components, Phys. Plasmas 17 (11) (2010) 112315.
  * [4] C. Wang, Y. Long, X.-T. He, J.-F. Wu, W.-H. Ye, P. Zhang, Transport properties of dense deuterium-tritium plasmas, Phys. Rev. E 88 (2013) 013106.
  * [5] D. Li, C. Wang, W. Kang, J. Yan, P. Zhang, Transport properties of hydrogen-helium mixtures at extreme density and temperature conditions, Phys. Rev. E 92 (2015) 043108.
  * [6] P. André, J. Aubreton, S. Clain, E. Duffour, M. E. Elchinger, B. Izrar, D. Rochette, R. Touzani, D. Vacher, Transport coefficients in thermal plasma. Applications to mars and titan atmospheres., Eur. Phys. J. D 57 (2010) 227.
  * [7] L. Biolsi, Transport properties in the jovian atmosphere, J. Geophys. Res. 83 (A3) (1978) 1125.
  * [8] J. Aubreton, M.-F. Elchinger, A. Hacala, U. Michon, Transport coefficients of typical biomass equimolar CO-H2 plasma, J. Phys. D 42 (9) (2009) 095206.
  * [9] S. Stauss, H. Muneoka, K. Terashima, Review on plasmas in extraordinary media: plasmas in cryogenic conditions and plasmas in supercritical fluids, Plasma Sources Sci. Technol. 27 (2) (2018) 023003.
  * [10] Y. Noma, J. Hyuk Choi, H. Muneoka, K. Terashima, Electron density and temperature of gas-temperature-dependent cryoplasma jet, J. Appl. Phys. 109 (5) (2011) 053303.
  * [11] D. Schaefer, M. Kohns, Molecular dynamics study of ion clustering in concentrated electrolyte solutions for the estimation of salt solubilities, Fluid Phase Equilibria 571 (2023) 113802.
  * [12] R. Takahashi, H. Nagashima, T. Tokumasu, S. Watanabe, S. Tsuda, Investigation of quantum effect of liquid hydrogen on homogeneous bubble nucleation using a density functional theory and molecular dynamics simulations, Fluid Phase Equilibria 553 (2022) 113300.
  * [13] M. Yiannourakou, X. Rozanska, B. Minisini, F. de Meyer, Molecular simulations for improved process modeling of an acid gas removal unit, Fluid Phase Equilibria 560 (2022) 113478.
  * [14] D. González-Salgado, E. Noya, E. Lomba, Simulation and theoretical analysis of the origin of the temperature of maximum density of water, Fluid Phase Equilibria 560 (2022) 113515.
  * [15] M. A. Galvani Cunha, M. O. Robbins, Determination of pressure-viscosity relation of 2,2,4-trimethylhexane by all-atom molecular dynamics simulations, Fluid Phase Equilibria 495 (2019) 28.
  * [16] S. S. Mishra, S. Bhattacharjee, Temperature guided behavioral transitions in confined helium: Gas-wall interaction effects on dynamics and transport in the cryogenic limit, Chemical Thermodynamics and Thermal Analysis 7 (2022) 100073\.
  * [17] L. d. S. Oliveira, P. A. Greaney, Method to manage integration error in the green-kubo method, Phys. Rev. E 95 (2017) 023308.
  * [18] S. Plumari, A. Puglisi, F. Scardina, V. Greco, Shear viscosity of a strongly interacting system: Green-kubo correlator versus chapman-enskog and relaxation-time approximations, Phys. Rev. C 86 (2012) 054902.
  * [19] H. Nagashima, S. ichi Tsuda, N. Tsuboi, M. Koshi, A. K. Hayashi, T. Tokumasu, An analysis of the quantum effect on the thermodynamic and transport properties of cryogenic hydrogen using molecular dynamics method, J. Phys. Conf. Ser 490 (1) (2014) 012160.
  * [20] M. Abbaspour, A. Keyvanloo, Many-body and quantum effects in some thermodynamic properties and infinite shear modulus of hfd-like fluid using the radial distribution function, J. Mol. Liq. 177 (2013) 1.
  * [21] M. S. Green, Markoff random processes and the statistical mechanics of time-dependent phenomena. II. irreversible processes in fluids, J. Chem. Phys. 22 (3) (1954) 398.
  * [22] R. Kubo, M. Yokota, S. Nakajima, Statistical-mechanical theory of irreversible processes. II. response to thermal disturbance, J. Phys. Soc. Japan 12 (11) (1957) 1203.
  * [23] J. Ferziger, H. Kaper, A. Brown, H. Kaper, Mathematical Theory of Transport Processes in Gases, North-Holland Publishing Company, 1972.
  * [24] B. Hess, Determining the shear viscosity of model liquids from molecular dynamics simulations, J. Chem. Phys. 116 (1) (2002) 209.
  * [25] D. Surblys, H. Matsubara, G. Kikugawa, T. Ohara, Methodology and meaning of computing heat flux via atomic stress in systems with constraint dynamics, J. Appl. Phys. 130 (21) (2021) 215104.
  * [26] M. Allen, D. Tildesley, Computer Simulation of Liquids, Oxford Science Publ, Clarendon Press, 1989.
  * [27] S. Chapman, T. Cowling, The Mathematical Theory of Non-uniform Gases: An Account of the Kinetic Theory of Viscosity, Thermal Conduction, and Diffusion in Gases, Cambridge University Press, 1958.
  * [28] J. De Boer, Chapter XVIII Transport properties of gaseous helium at low temperatures, Vol. 1 of Progress in Low Temperature Physics, Elsevier, 1955, p. 381.
  * [29] J. O. Hirschfelder, C. F. Curtiss, R. B. Bird, Molecular theory of gases and liquids., J. Polym. Sci. 17 (83) (May 1955).
  * [30] V. D. Arp, M. R. D, Thermophysical properties of Helium-4 from 0.8 to 1500 K with pressures to 2000 MPa, NIST Technical Note 1334, NIST, Boulder, 1989\.
  * [31] L. Monchick, E. A. Mason, R. J. Munn, F. J. Smith, Transport properties of gaseous He3 and He4, Phys. Rev. 139 (1965) A1076.
  * [32] V. J. Emery, Spin diffusion in gases at low temperatures, Phys. Rev. 133 (1964) A661.
  * [33] A. K. Varshneya, J. C. Mauro, Chapter 14 - Permeation, diffusion, and ionic conduction in glass, in: Fundamentals of Inorganic Glasses, Elsevier, 2019, pp. 383–424.
  * [34] M. P. Lautenschlaeger, H. Hasse, Transport properties of the lennard-jones truncated and shifted fluid from non-equilibrium molecular dynamics simulations, Fluid Phase Equilibria 482 (2019) 38.
  * [35] J. Hurly, M. Moldover, Ab initio values of the thermophysical properties of helium as standards, J. Res. Natl. Inst. Stand. Technol. (105) (2000) 667.
  * [36] A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, S. J. Plimpton, LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales, Comp. Phys. Comm. 271 (2022) 108171.
  * [37] A. J. H. McGaughey, M. Kaviany, Quantitative validation of the boltzmann transport equation phonon thermal conductivity model under the single-mode relaxation time approximation, Phys. Rev. B 69 (2004) 094303.
  * [38] M. Pillai, J. Goglio, T. G. Walker, Matrix Numerov method for solving Schrodinger’s equation, Am. J. Phys. 80 (11) (2012) 1017.
  * [39] J. J. Sakurai, Modern Quantum Mechanics, Pearson, 1982.
  * [40] R. Feynman, A. Hibbs, D. Styer, Quantum Mechanics and Path Integrals, Dover Books on Physics, Dover Publications, 2010.
  * [41] L. M. Sesé, Thermodynamic and structural properties of the path-integral quantum hard-sphere fluid, J. Chem. Phys. 108 (21) (1998) 9086.
  * [42] A. V. A. Kumar, H. Jobic, S. K. Bhatia, Quantum effects on adsorption and diffusion of hydrogen and deuterium in microporous materials, J. Phys. Chem. B 110 (33) (2006) 16666.
  * [43] B. Guillot, Y. Guissani, Quantum effects in simulated water by the feynman–hibbs approach, J. Chem. Phys. 108 (24) (1998) 10162.
  * [44] M. Abbaspour, Z. Borzouie, Equation of state, elastic constants, and melting curve of solid neon using an effective two-body potential including quantum corrections, Fluid Phase Equilibria 379 (2014) 167.
  * [45] N. Tchouar, F. Ould-Kaddour, D. Levesque, Computation of the properties of liquid neon, methane, and gas helium at low temperature by the feynman-hibbs approach, J. Chem. Phys. 121 (15) (2004) 7326.
  * [46] L. M. Sesé, Path Integrals and Effective Potentials in the Study of Monatomic Fluids at Equilibrium, John Wiley & Sons, Ltd, pp. 49–158.


After some elementary calculations, we obtain as $n\rightarrow\infty,$

$L_{1}-L_{0}\rightarrow-Mc+\sum_{k=1}^{K}d_{k}(\phi_{k}-\log\phi_{k}-1),$
$\varsigma_{l}\rightarrow\sqrt{(\alpha_{x}+1)(-\log(1-c)-c)},$
$\breve{\varsigma}_{l}-\sqrt{(\alpha_{x}+1)(-\log(1-c)-c)+\sum_{k=1}^{K}\frac{\left(\phi_{n}\left(\alpha_{k}\right)-1\right)^{2}}{n}}s_{k}^{2}\rightarrow
0.$

Therefore, we have as $n$ tends to infinity,

$\displaystyle
P_{L}-\Phi\left(\frac{\sum_{k=1}^{K}d_{k}\left(\phi_{n}(\alpha_{k})-\log\phi_{n}(\alpha_{k})\right)-M-Mc-
z_{\xi}\varsigma_{l}}{\sqrt{(\alpha_{x}+1)(-\log(1-c)-c)+\sum_{k=1}^{K}\frac{\left(\phi_{n}\left(\alpha_{k}\right)-1\right)^{2}}{n}}s_{k}^{2}}\right)\rightarrow
0,$ (7.35)

then the proof of Theorem 4.3 is finished.

###  7.6 Proof of Theorem 4.4

Since
$\xi=P\left(W>z_{\xi}\varsigma_{w}+c_{n}\left(p+\alpha_{x}+\beta_{x}\right)\right),$
for brevity, we use the notation $W_{0}=pc_{n}+c_{n}(\alpha_{x}+\beta_{x}),$
$W_{1}=(p-M)\int f_{W}(x)dF^{c_{nM},H_{2n}}+\breve{\mu}_{w}.$ Therefore, the
power to detect the hypothesis is

$\displaystyle P_{W}=P\left(W>z_{\xi}\varsigma_{w}+W_{0}\right)$
$\displaystyle=P\left(\frac{W-W_{1}}{\breve{\varsigma}_{w}}>\frac{z_{\xi}\varsigma_{w}+W_{0}-W_{1}}{\breve{\varsigma}_{w}}\right)$

Since $\frac{W-W_{1}}{\breve{\varsigma}_{w}}$ is asymptotically normal
distributed, then $P_{W}$ is approxiamted to
$\Phi(\frac{W_{1}-W_{0}}{\breve{\varsigma}_{w}}-z_{\xi}\frac{\varsigma_{w}}{\breve{\varsigma}_{w}}).$
Here
$W_{1}-W_{0}=(p-M)c_{nM}-pc_{n}+(\beta_{x}+\alpha_{x})c_{nM}-(\beta_{x}+\alpha_{x})c_{n}+\sum_{k=1}^{K}d_{k}(\phi_{n}\left(\alpha_{k}\right)-1)^{2}-Mc_{nM}^{2}.$
After some elementary calculations, we obtain, as $n\rightarrow\infty$,

$W_{1}-W_{0}\rightarrow\sum_{k=1}^{K}d_{k}\left(\phi_{n}\left(\alpha_{k}\right)-1\right)^{2}-Mc^{2}-2Mc,$
$\varsigma_{w}\rightarrow\sqrt{(\alpha_{x}+1)(4c^{3}+2c^{2})+4\beta_{x}c^{3}},$
$\breve{\varsigma}_{w}-\sqrt{(\alpha_{x}+1)(4c^{3}+2c^{2})+4\beta_{x}c^{3}+\sum_{k=1}^{K}\dfrac{4\phi_{n}^{2}\left(\alpha_{k}\right)\left(\phi_{n}\left(\alpha_{k}\right)-1\right)^{2}}{n}s_{k}^{2}}\rightarrow
0.$

Then we have as $n$ tends to infinity,

$\displaystyle
P_{W}-\Phi\left(\frac{\sum_{k=1}^{K}d_{k}\left(\phi_{k}-1\right)^{2}-Mc^{2}-2Mc-
z_{\xi}\varsigma_{w}}{\sqrt{(\alpha_{x}+1)(4c^{3}+2c^{2})+4\beta_{x}c^{3}+\sum_{k=1}^{K}\dfrac{4\phi_{n}^{2}\left(\alpha_{k}\right)\left(\phi_{n}\left(\alpha_{k}\right)-1\right)^{2}}{n}s_{k}^{2}}}\right)\rightarrow
0,$ (7.36)

then the proof is finished.

###  7.7 Proof of Theorem 4.5

From Jiang and Bai (2021), for spike $\alpha_{1}$, we eliminate the
multiplicity of it and then we have

$\displaystyle\sqrt{\frac{n\theta_{1}^{2}}{2\theta_{1}+\sum_{t=1}^{p}\left|u_{t1}\right|^{4}\beta_{x}\nu_{1}}}\frac{\lambda_{1}-\phi_{n}\left(\alpha_{1}\right)}{\phi_{n}\left(\alpha_{1}\right)}\stackrel{{\scriptstyle
d}}{{\longrightarrow}}N\left(0,1\right).$

Then the power of test R equals

$\displaystyle P_{R}$
$\displaystyle=P\left(\lambda_{1}>t_{\xi}\varsigma_{r}+\mu_{r}\right)$
$\displaystyle=P\left(\sqrt{\frac{n\theta_{1}^{2}}{2\theta_{1}+\sum_{t=1}^{p}\left|u_{t1}\right|^{4}\beta_{x}\nu_{1}}}\frac{\lambda_{1}-\phi_{n}\left(\alpha_{1}\right)}{\phi_{n}\left(\alpha_{1}\right)}>\right.$
$\displaystyle\left.\sqrt{\frac{n\theta_{1}^{2}}{2\theta_{1}+\sum_{t=1}^{p}\left|u_{t1}\right|^{4}\beta_{x}\nu_{1}}}\frac{t_{\xi}\varsigma_{r}+\mu_{r}-\phi_{n}\left(\alpha_{1}\right)}{\phi_{n}\left(\alpha_{1}\right)}\right)$

Since
$\sqrt{\frac{n\theta_{1}^{2}}{2\theta_{1}+\sum_{t=1}^{p}\left|u_{t1}\right|^{4}\beta_{x}\nu_{1}}}\frac{\lambda_{1}-\phi_{n}\left(\alpha_{1}\right)}{\phi_{n}\left(\alpha_{1}\right)}$
is asymptotically standard normal distributed, then $P_{R}$ is approximate to
$1-\Phi\left(\sqrt{\frac{n\theta_{1}^{2}}{2\theta_{1}+\sum_{t=1}^{p}\left|u_{t1}\right|^{4}\beta_{x}\nu_{1}}}\frac{t_{\xi}\varsigma_{r}+\mu_{r}-\phi_{n}\left(\alpha_{1}\right)}{\phi_{n}\left(\alpha_{1}\right)}\right),$
and it equals
$\Phi(-\sqrt{n}\frac{t_{\xi}\varsigma_{r}+\mu_{r}-\phi_{n}\left(\alpha_{1}\right)}{s_{1}\phi_{n}\left(\alpha_{1}\right)}),$
then the proof is finished.

##  8 Some deviations and calculations

This section contains proof of formulas stated in the proof of Theorems 4.1
and 4.2. We begin by deriving formula (7.32). First, we consider
$\oint_{\mathcal{C}}f_{L}\left(z\right)\frac{\underline{m}^{\prime}(z)}{\underline{m}(z)}dz$.

$\displaystyle\oint_{\mathcal{C}}f_{L}\left(z\right)\frac{\underline{m}^{\prime}(z)}{\underline{m}(z)}dz=\oint_{\mathcal{C}}f_{L}\left(z\right)d\log\underline{m}\left(z\right)=-\oint_{\mathcal{C}}f_{L}^{{}^{\prime}}\left(z\right)\log\underline{m}\left(z\right)dz$
$\displaystyle=$
$\displaystyle\int_{a(c)}^{b(c)}f_{L}^{{}^{\prime}}\left(z\right)\left[\log\underline{m}(x+i\varepsilon)-\log\underline{m}(x-i\varepsilon)\right]dx$
$\displaystyle=$ $\displaystyle
2i\int_{a(c)}^{b(c)}f_{L}^{{}^{\prime}}\left(z\right)\Im\log\underline{m}(x+i\varepsilon)dx$
(8.37)

Here, $a(c)=(1-\sqrt{c})^{2}$ and $b(c)=(1+\sqrt{c})^{2}$. Since

$\underline{m}\left(z\right)=-\frac{1-c}{z}+cm\left(z\right),$

under $H_{1}$, we have

$\underline{m}\left(z\right)=\frac{-(z+1-c)+\sqrt{(z-1-c)^{2}-4c}}{2z}.$

As $z\rightarrow x\in$ $[a(c),b(c)]$, we obtain

$\underline{m}\left(x\right)=\frac{-(x+1-c)+\sqrt{4c-(x-1-c)^{2}}i}{2x}.$

Therefore,

$\displaystyle\int_{a(c)}^{b(c)}f_{L}^{{}^{\prime}}\left(z\right)\Im\log\underline{m}(x+i\varepsilon)dx$
$\displaystyle=\int_{a(c)}^{b(c)}f_{L}^{{}^{\prime}}(x)\tan^{-1}\left(\frac{\sqrt{4c-(x-1-c)^{2}}}{-(x+1-c)}\right)dx$
$\displaystyle=\left[\left.\tan^{-1}\left(\frac{\sqrt{4c-(x-1-c)^{2}}}{-(x+1-c)}\right)f_{L}(x)\right|_{a(c)}^{b(c)}-\int_{a(c)}^{b(c)}f_{L}(x)d\tan^{-1}\left(\frac{\sqrt{4c-(x-1-c)^{2}}}{-(x+1-c)}\right)\right].$

It is easy to verify that the first term is $0$, and we now focus on the
second term,

$\displaystyle\int_{a(c)}^{b(c)}f_{L}(x)d\tan^{-1}\left(\frac{\sqrt{4c-(x-1-c)^{2}}}{-(x+1-c)}\right)$
$\displaystyle=\int_{a(c)}^{b(c)}\frac{\left(x-\log
x-1\right)}{1+\frac{4c-(x-1-c)^{2}}{(x+1-c)^{2}}}\cdot\frac{\sqrt{4c-(x-1-c)^{2}}+\frac{(x-1-c)(x+1-c)}{\sqrt{4c-(x-1-c)^{2}}}}{(x+1-c)^{2}}dx.$
(8.38)

By substituting $x=1+c-2\sqrt{c}\cos(\theta)$, we obtain

$\displaystyle(\ref{27})$
$\displaystyle=\frac{1}{2}\int_{0}^{2\pi}\left(1+c-2\sqrt{c}\cos(\theta)-\log\left(1+c-2\sqrt{c}\cos(\theta)\right)-1\right)\frac{c-\sqrt{c}\cos(\theta)}{1+c-2\sqrt{c}\cos(\theta)}d\theta$
$\displaystyle=\frac{1}{2}\int_{0}^{2\pi}\left[1-\frac{\log\left(1+c-2\sqrt{c}\cos(\theta)\right)+1}{1+c-2\sqrt{c}\cos(\theta)}\right]\left(c-\sqrt{c}\cos(\theta)\right)d\theta$
$\displaystyle=\frac{1}{2}\int_{0}^{2\pi}\left(c-\sqrt{c}\cos(\theta)\right)d\theta-\frac{1}{2}\int_{0}^{2\pi}\frac{\log\left(1+c-2\sqrt{c}\cos(\theta)\right)}{1+c-2\sqrt{c}\cos(\theta)}\left(c-\sqrt{c}\cos(\theta)\right)d\theta-$
(8.39)
$\displaystyle\frac{1}{2}\int_{0}^{2\pi}\frac{c-\sqrt{c}\cos(\theta)}{1+c-2\sqrt{c}\cos(\theta)}d\theta$

It is easy to obtain that the first term of (8.39) is $\pi c$; then, we
consider the second term. By substituting $\cos\theta=\frac{z+z^{-1}}{2}$, we
turn it into a contour integral on $\left|z\right|=1$

$\displaystyle\frac{1}{2}\int_{0}^{2\pi}\frac{\log\left(1+c-2\sqrt{c}\cos(\theta)\right)}{1+c-2\sqrt{c}\cos(\theta)}\left(c-\sqrt{c}\cos\theta\right)d\theta$
$\displaystyle=\frac{1}{2}\oint_{\left|z\right|=1}\log|1-\sqrt{c}z|^{2}\cdot\frac{c-\sqrt{c}\frac{z+z^{-1}}{2}}{1+c-2\sqrt{c}\cdot\frac{z+z^{-1}}{2}}\frac{dz}{iz}$
$\displaystyle=\frac{1}{4i}\oint_{\left|z\right|=1}\log|1-\sqrt{c}z|^{2}\cdot\frac{2cz-\sqrt{c}\left(z^{2}+1\right)}{(z-\sqrt{c})(-\sqrt{c}z+1)z}dz$

When $c<1$, $0$ and $\sqrt{c}$ are poles, by using the residue theorem, the
integral is $-\pi\log(1-c)$. The same argument also holds for the third term,
and the integral is $0$ after some calculation.

Therefore,

$\frac{M}{2\pi
i}\oint_{\mathcal{C}}f_{L}\left(z\right)\frac{\underline{m}^{\prime}(z)}{\underline{m}(z)}dz=-M(c+\log(1-c)),$

and the result is still valid if $c$ is replaces $c_{nM}$; therefore, formula
(7.32) holds.

Now, we prove (7.34). Since
$z=-\frac{1}{\underline{m}}+\frac{c}{1+\underline{m}}$, we have, for $c>1$,

$\displaystyle\oint_{\mathcal{C}}f_{W}\left(z\right)\frac{\underline{m}^{\prime}(z)}{\underline{m}(z)}dz=\oint_{\mathcal{C}_{1}}f_{W}\left(z\right)\frac{\underline{m}^{\prime}(z)}{\underline{m}(z)}dz+\oint_{\mathcal{C}_{2}}f_{W}\left(z\right)\frac{\underline{m}^{\prime}(z)}{\underline{m}(z)}dz,$

where $\mathcal{C}_{1}$ is a contour that includes the interval
$((1-\sqrt{c})^{2},(1+\sqrt{c})^{2})$, and $\mathcal{C}_{2}$ is a contour that
includes the origin. Using $\mathcal{C}_{\underline{m}}$ to denote the contour
of $\underline{m}$, we obtain

$\displaystyle\oint_{\mathcal{C}_{1}}f_{W}\left(z\right)\frac{\underline{m}^{\prime}(z)}{\underline{m}(z)}dz=\oint_{\mathcal{C}_{\underline{m}}}(-\frac{1}{\underline{m}}+\frac{c}{1+\underline{m}}-1)^{2}\frac{\underline{m}^{\prime}(z)}{\underline{m}(z)}\frac{dz}{d\underline{m}}d\underline{m}$
$\displaystyle=$
$\displaystyle\oint_{\mathcal{C}_{\underline{m}}}(-\frac{1+\underline{m}}{\underline{m}}+\frac{c}{1+\underline{m}})^{2}\frac{1}{\underline{m}}d\underline{m}=\oint_{\mathcal{C}_{\underline{m}}}(\frac{(1+\underline{m})^{2}}{\underline{m}^{3}}+\frac{c^{2}}{(1+\underline{m})^{2}\underline{m}}-\frac{2c}{\underline{m}^{2}})d\underline{m}$

Since the $z$ contour cannot enclose the origin, neither can the resulting
$\underline{m}$ contour. Thus, the only pole is $-1$, the residue is $-c^{2}$
by residue theorem, and we obtain the integral as $-2\pi ic^{2}.$

Then, we focus on the second integral
$\oint_{\mathcal{C}_{2}}f_{W}\left(z\right)\frac{\underline{m}^{\prime}(z)}{\underline{m}(z)}dz.$
When $z=0,$ we obtain $\underline{m}=\frac{1}{c-1}$; since $c>1,$
$\frac{1}{c-1}>0.$ Both $\underline{m}=0$ and $\underline{m}=-1$ are not in
the contour. Thus, the integrand
$(\frac{(1+\underline{m})^{2}}{\underline{m}^{3}}+\frac{c^{2}}{(1+\underline{m})^{2}\underline{m}}-\frac{2c}{\underline{m}^{2}})$
is analytic in the contour. The integral is $0$. Therefore, when $c>1$,
$\frac{M}{2\pi
i}\oint_{\mathcal{C}}f_{W}\left(z\right)\frac{\underline{m}^{\prime}(z)}{\underline{m}(z)}dz=-Mc^{2}$.
When $c<1,$, the contour integral
$\oint_{\mathcal{C}}f_{W}\left(z\right)\frac{\underline{m}^{\prime}(z)}{\underline{m}(z)}dz$
reduces to
$\oint_{\mathcal{C}_{1}}f_{W}\left(z\right)\frac{\underline{m}^{\prime}(z)}{\underline{m}(z)}dz$,
and the result is also the same as above. When $c=1$, the result is still true
by continuity in $c$. The results above are still valid if $c$ is replaced by
$c_{nM}.$ Therefore, the proof of (7.34) is complete.

##  9 Some useful lemmas

######  Lemma 9.1.

If ${\bf D}_{2}={\bf I}_{p-M}$, then the mean function $\mu_{1}$ and
$\kappa_{nst}$ in the covariance function of Theorem 3.2 can be simplified
from the results in Wang and Yao (2013) and Zheng et al. (2015), i.e.,

$\mu_{1}=\alpha_{x}I_{1}(f_{1})+\beta_{x}I_{2}(f_{1}),$
$\kappa_{nst}=(\alpha_{x}+1)J_{1}(f_{s},f_{t})+\beta_{x}J_{2}(f_{s},f_{t}),$
$\displaystyle I_{1}(f_{1})$ $\displaystyle=\lim_{r\downarrow 1}\frac{1}{2\pi
i}\oint_{|z|=1}f_{1}\left(\left|1+\sqrt{c_{nM}}z\right|^{2}\right)\left[\frac{z}{z^{2}-r^{-2}}-\frac{1}{z}\right]dz,$
$\displaystyle I_{2}(f_{1})$ $\displaystyle=\frac{1}{2\pi
i}\oint_{|z|=1}f_{1}\left(\left|1+\sqrt{c_{nM}}z\right|^{2}\right)\frac{1}{z^{3}}dz,$
$\displaystyle J_{1}(f_{s},f_{t})$ $\displaystyle=\lim_{r\downarrow
1}\frac{-1}{4\pi^{2}}\oint_{\left|z_{1}\right|=1}\oint_{\left|z_{2}\right|=1}\frac{f_{s}\left(\left|1+\sqrt{c_{nM}}z_{1}\right|^{2}\right)f_{t}\left(\left|1+\sqrt{c_{nM}}z_{2}\right|^{2}\right)}{\left(z_{1}-rz_{2}\right)^{2}}dz_{1}dz_{2},$
$\displaystyle J_{2}(f_{s},f_{t})$
$\displaystyle=-\frac{1}{4\pi^{2}}\oint_{\left|z_{1}\right|=1}\frac{f_{s}\left(\left|1+\sqrt{c_{nM}}z_{1}\right|^{2}\right)}{z_{1}^{2}}dz_{1}\oint_{\left|z_{2}\right|=1}\frac{f_{t}\left(\left|1+\sqrt{c_{nM}}z_{2}\right|^{2}\right)}{z_{2}^{2}}dz_{2}.$

######  Lemma 9.2.

Note that for any matrix ${\bf Z}$,

$\displaystyle{\bf Z}\left({\bf Z}^{\ast}{\bf Z}-\lambda{\bf
I}\right)^{-1}{\bf Z}^{\ast}={\bf I}+\lambda\left({\bf Z}{\bf
Z}^{\ast}-\lambda{\bf I}\right)^{-1}.$

##  10 Tables for simulation studies

In this section, we present addtional simulation tables regarding empirical
probability of rejecting alternative hypotheses in Section 5 in the main file.

Table 7: Empirical probability of rejecting $H_{1}$ at significance level $\xi=0.01$ under assumptions of Gaussian, Gamma, and Uniform distributions  |  | ($Dt_{1},H_{1}$) | ($Dt_{2},H_{1}$) | ($Dt_{3},H_{1}$)
---|---|---|---|---
test |  (p,n) | $\alpha_{1}=3$ | $\alpha_{1}=5$ | $\alpha_{1}=7$ | $\alpha_{1}=3$ | $\alpha_{1}=5$ | $\alpha_{1}=7$ | $\alpha_{1}=3$ | $\alpha_{1}=5$ | $\alpha_{1}=7$
$CLRT$ | (50,150) | 0.5142 | 0.9954 | 1 | 0.4896 | 0.9896 | 1 | 0.5130 | 0.9995 | 1
| (100,300) | 0.5056 | 0.9992 | 1 | 0.4954 | 0.9972 | 1 | 0.5134 | 0.9999 | 1
| (200,600) | 0.5178 | 0.9996 | 1 | 0.5085 | 0.9990 | 1 | 0.5210 | 1 | 1
$CNTT$ | (50,150) | 0.9814 | 1 | 1 | 0.9287 | 1 | 1 | 0.9984 | 1 | 1
| (100,300) | 0.9938 | 1 | 1 | 0.9662 | 1 | 1 | 0.9995 | 1 | 1
| (200,600) | 0.9985 | 1 | 1 | 0.9873 | 1 | 1 | 1 | 1 | 1
$RLRT$ | (50,150) | 0.9983 | 1 | 1 | 0.9947 | 1 | 1 | 1 | 1 | 1
| (100,300) | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
| (200,600) | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
Table 8: Empirical probability of rejecting $H_{2}$ at significance level $\xi=0.01$ under assumptions of Gaussian, Gamma, and Uniform distributions  |  | ($Dt_{1},H_{2}$) | ($Dt_{2},H_{2}$) | ($Dt_{3},H_{2}$)
---|---|---|---|---
test |  (p,n) | $\alpha_{1}=3$ | $\alpha_{1}=5$ | $\alpha_{1}=7$ | $\alpha_{1}=3$ | $\alpha_{1}=5$ | $\alpha_{1}=7$ | $\alpha_{1}=3$ | $\alpha_{1}=5$ | $\alpha_{1}=7$
$CLRT$ | (50,150) | 0.9387 | 1 | 1 | 0.9045 | 1 | 1 | 0.9574 | 1 | 1
| (100,300) | 0.9496 | 1 | 1 | 0.9330 | 1 | 1 | 0.9666 | 1 | 1
| (200,600) | 0.9624 | 1 | 1 | 0.9553 | 1 | 1 | 0.9710 | 1 | 1
$CNTT$ | (50,150) | 1 | 1 | 1 | 0.9998 | 1 | 1 | 1 | 1 | 1
| (100,300) | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
| (200,600) | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
$RLRT$ | (50,150) | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
| (100,300) | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
| (200,600) | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
Table 9: Empirical probability of rejecting $H_{4}$ at significance level $\xi=0.05$ under assumptions of Gaussian, Gamma, and Uniform distributions  |  | ($Dt_{1},H_{4}$) | ($Dt_{2},H_{4}$) | ($Dt_{3},H_{4}$)
---|---|---|---|---
test |  (p,n) | $\alpha_{1}=3$ | $\alpha_{1}=5$ | $\alpha_{1}=7$ | $\alpha_{1}=3$ | $\alpha_{1}=5$ | $\alpha_{1}=7$ | $\alpha_{1}=3$ | $\alpha_{1}=5$ | $\alpha_{1}=7$
$CLRT$ | (50,150) | 0.7288 | 0.9991 | 1 | 0.7074 | 0.9989 | 1 | 0.7269 | 0.9993 | 1
| (100,300) | 0.7397 | 0.9999 | 1 | 0.7364 | 0.9997 | 1 | 0.7411 | 0.9997 | 1
| (200,600) | 0.7568 | 1 | 1 | 0.7484 | 1 | 1 | 0.7649 | 1 | 1
$CNTT$ | (50,150) | 0.9917 | 1 | 1 | 0.9821 | 1 | 1 | 0.9974 | 1 | 1
| (100,300) | 0.9985 | 1 | 1 | 0.9951 | 1 | 1 | 0.9997 | 1 | 1
| (200,600) | 0.9997 | 1 | 1 | 0.9983 | 1 | 1 | 1 | 1 | 1
$RLRT$ | (50,150) | 0.9997 | 1 | 1 | 0.9992 | 1 | 1 | 0.9996 | 1 | 1
| (100,300) | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
| (200,600) | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
Table 10: Empirical probability of rejecting $H_{4}$ at significance level $\xi=0.01$ under assumptions of Gaussian, Gamma, and Uniform distributions  |  | ($Dt_{1},H_{4}$) | ($Dt_{2},H_{4}$) | ($Dt_{3},H_{4}$)
---|---|---|---|---
test |  (p,n) | $\alpha_{1}=3$ | $\alpha_{1}=5$ | $\alpha_{1}=7$ | $\alpha_{1}=3$ | $\alpha_{1}=5$ | $\alpha_{1}=7$ | $\alpha_{1}=3$ | $\alpha_{1}=5$ | $\alpha_{1}=7$
$CLRT$ | (50,150) | 0.5068 | 0.9967 | 1 | 0.5043 | 0.9952 | 1 | 0.5101 | 0.9966 | 1
| (100,300) | 0.5137 | 0.9996 | 1 | 0.5082 | 0.9987 | 1 | 0.5193 | 0.9996 | 1
| (200,600) | 0.5155 | 0.9995 | 1 | 0.5192 | 0.9997 | 1 | 0.5100 | 0.9998 | 1
$CNTT$ | (50,150) | 0.9782 | 1 | 1 | 0.9469 | 1 | 1 | 0.9934 | 1 | 1
| (100,300) | 0.9935 | 1 | 1 | 0.9803 | 1 | 1 | 0.9988 | 1 | 1
| (200,600) | 0.9984 | 1 | 1 | 0.9921 | 1 | 1 | 0.9997 | 1 | 1
$RLRT$ | (50,150) | 0.9978 | 1 | 1 | 0.9991 | 1 | 1 | 0.9981 | 1 | 1
| (100,300) | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
| (200,600) | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
Table 11: Empirical probability of rejecting $H_{5}$ at significance level $\xi=0.05$ under assumptions of Gaussian, Gamma, and Uniform distributions  |  | ($Dt_{1},H_{5}$) | ($Dt_{2},H_{5}$) | ($Dt_{3},H_{5}$)
---|---|---|---|---
test |  (p,n) | $\alpha_{1}=3$ | $\alpha_{1}=5$ | $\alpha_{1}=7$ | $\alpha_{1}=3$ | $\alpha_{1}=5$ | $\alpha_{1}=7$ | $\alpha_{1}=3$ | $\alpha_{1}=5$ | $\alpha_{1}=7$
$CLRT$ | (50,150) | 0.9803 | 1 | 1 | 0.9778 | 1 | 1 | 0.9819 | 1 | 1
| (100,300) | 0.9889 | 1 | 1 | 0.9864 | 1 | 1 | 0.9861 | 1 | 1
| (200,600) | 0.9928 | 1 | 1 | 0.9913 | 1 | 1 | 0.9908 | 1 | 1
$CNTT$ | (50,150) | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
| (100,300) | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
| (200,600) | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
$RLRT$ | (50,150) | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
| (100,300) | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
| (200,600) | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
Table 12: Empirical probability of rejecting $H_{5}$ at significance level $\xi=0.01$ under assumptions of Gaussian, Gamma, and Uniform distributions  |  | ($Dt_{1},H_{5}$) | ($Dt_{2},H_{5}$) | ($Dt_{3},H_{5}$)
---|---|---|---|---
test |  (p,n) | $\alpha_{1}=3$ | $\alpha_{1}=5$ | $\alpha_{1}=7$ | $\alpha_{1}=3$ | $\alpha_{1}=5$ | $\alpha_{1}=7$ | $\alpha_{1}=3$ | $\alpha_{1}=5$ | $\alpha_{1}=7$
$CLRT$ | (50,150) | 0.9321 | 1 | 1 | 0.9218 | 1 | 1 | 0.9418 | 1 | 1
| (100,300) | 0.9530 | 1 | 1 | 0.9521 | 1 | 1 | 0.9528 | 1 | 1
| (200,600) | 0.9610 | 1 | 1 | 0.9586 | 1 | 1 | 0.9638 | 1 | 1
$CNTT$ | (50,150) | 1 | 1 | 1 | 0.9998 | 1 | 1 | 1 | 1 | 1
| (100,300) | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
| (200,600) | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
$RLRT$ | (50,150) | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
| (100,300) | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
| (200,600) | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
Table 13: Empirical probability of rejecting $H_{6}$ at significance level $\xi=1\times 10^{-4}$ under assumptions of Gaussian, Gamma, and Uniform distributions  |  | ($Dt_{1},H_{6}$) | ($Dt_{2},H_{6}$) | ($Dt_{3},H_{6}$)
---|---|---|---|---
|  | $\alpha_{1}$
test |  (p,n) | $2.2$ | $2.5$ | $2.8$ | $2.2$ | $2.5$ | $2.8$ | $2.2$ | $2.5$ | $2.8$
$CLRT$ | (50,150) | 0.3876 | 0.8702 | 0.9943 | 0.6033 | 0.9347 | 0.9971 | 0.3881 | 0.8639 | 0.9946
| (100,300) | 0.3849 | 0.8870 | 0.9980 | 0.6069 | 0.9616 | 0.9993 | 0.3856 | 0.8917 | 0.9985
| (200,600) | 0.3761 | 0.9038 | 0.9978 | 0.6272 | 0.9698 | 0.9996 | 0.3842 | 0.9026 | 0.9990
$CNTT$ | (50,150) | 0.9974 | 1 | 1 | 0.9969 | 0.9999 | 1 | 0.9998 | 1 | 1
| (100,300) | 0.9996 | 1 | 1 | 0.9994 | 1 | 1 | 1 | 1 | 1
| (200,600) | 0.9998 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
$RLRT$ | (50,150) | 0.8845 | 0.9972 | 1 | 0.9111 | 0.9979 | 1 | 0.8718 | 0.9962 | 1
| (100,300) | 0.9833 | 1 | 1 | 0.9858 | 1 | 1 | 0.9802 | 1 | 1
| (200,600) | 0.9994 | 1 | 1 | 0.9995 | 1 | 1 | 0.9995 | 1 | 1

## References

  * Bai and Silverstein (2004) Zhidong Bai and Jack W. Silverstein.  CLT for linear spectral statistics of large-dimensional sample covariance matrices.  _The Annals of Probability_ , 32(1A):553 – 605, 2004.
  * Bai et al. (2009) Zhidong Bai, Dandan Jiang, Jianfeng Yao, and Shurong Zheng.  Corrections to LRT on large-dimensional covariance matrix by RMT.  _The Annals of Statistics_ , 37(6B):3822 – 3840, 2009.
  * Jiang and Bai (2021) Dandan Jiang and Zhidong Bai.  Generalized four moment theorem and an application to clt for spiked eigenvalues of high-dimensional covariance matrices.  _Bernoulli_ , 27(1):274–294, 2021\.
  * Liu et al. (2022) Zhijun Liu, Zhidong Bai, Jiang Hu and Haiyan Song.  A CLT for the LSS of large dimensional sample covariance matrices with diverging spikes.
  * Wang and Yao (2013) Qinwen Wang and Jianfeng Yao.  On the sphericity test with large-dimensional observations.  _Electronic Journal of Statistics_ , 7:2164–2192, 2013\.
  * Zheng et al. (2015) Shurong Zheng, Zhidong Bai, and Jianfeng Yao.  Substitution principle for CLT of linear spectral statistics of high-dimensional sample covariance matrices with applications to hypothesis testing.  _The Annals of Statistics_ , 43(2):546 – 591, 2015.


# An Optimal Constraint for QUBO Models

Clark Alexander111Chicago Quantum
email: the author

###### Abstract

A quadratic binary unconstrained optimization model, hereafter QUBO, by
definition is unconstrained. This, however, is not ideal if one needs to
select a model containing only a fixed size binary vector. In this work we
show how to add a constraint to a QUBO to force a particular size solution.

##  1 Introduction

As the industry of quantum computing has been making inroads in applications,
the first set of problems which are in-line for fast, efficient “solutions”
are combinatorial optimization problems [CQ1, DwaveI, STB] As of this writing,
the technique which appears to be furthest along is quantum annealing. This is
an appealing technique as it allows one with a reasonably powerful laptop and
a high level language compiler to benchmark the quantum annealing process.
Generally speaking, the most-straightforward way to apply quantum annealing is
to pose one’s problem in terms of a quadratic form. In particular, if one can
pose a problem in terms of a solution which marks a set of elements as “in” or
“out” then quantum annealing works well. Additionally, there are many
classical, heuristic, and probabilistic techniques which one can use to
benchmark a quantum annealer; for example, simulated annealing, genetic
algorithms, simulated bifurcation machines, branch-and-bound, etc. [DA, GKD,
SBM, EKU].

However, quadratic forms with “in”/“out” solutions require binary vectors and
depending on the particular type of quadratic form given, one can expect a
very specific set of solutions to show up. For example, if one provides a
positive definite matrix as a quadratic form then the known minimum is the
zero vector. There is no need to apply probabilistic techniques, all non zero
vectors will give positive solutions and in searching for a minimum, we cannot
go lower than zero. Additionally if the quadratic form is a well known random
matrix, a Gaussian ensemble, for instance, then the best solution will contain
roughly half “in” and half “out.” These solutions are expected, but maybe not
completely useful.

Recently, there has been a push for using quadratic forms in financial
optimization [CQ3, CQ2, MLO] One common strategy is to maximize expected
return while minimizing volatility. Volatility is, however, covariance, which
is positive semi-definite. In reality covariance is positive definite with
real data. Thus trying to minimize a positive definite matrix reveals the
“mathematically correct” solution of not investing as that reduces volatility
to zero. However, for an investor, this is not a useful solution. Perhaps an
investor wants to invest in an “optimal” set of 15 assets. An unconstrained
problem will never produce such an answer, thus one wishes to add a constraint
in a somewhat natural way so that a set of 15 assets becomes visible.

##  2 Elementary Properties of a Random QUBO

A quadratic unconstrained binary optimization (hereafter QUBO) can be thought
of as a square matrix. While this is a slight abuse of notation, there should
not be too much confusion as the actual optimization involves finding a vector
for the particular matrix. So where it is clear we will refer to the matrix
and the optimization problem both as QUBO.

######  Definition 1.

Given a real symmetric matrix $A\in\mathbb{R}^{N\times N}$ and a real vector
$B\in\mathbb{R}^{N}$, a quadratic unconstrained binary optimization is a
problem in which one seeks the vector $x\in\\{0,1\\}^{N}$ so that

$Q(A,B)=x^{t}Ax+B\cdot x$ (1)

achieves a minimum value. That is

$\min_{x}x^{t}Ax+B\cdot x$ (2)

It is also common to expect $A$ to be upper or lower triangular and traceless.
However by considering $(A+A^{t})/2$ one can trivially move back and forth
between triangular and symmetric matrices. We prefer the symmetric matrix as
it guarantees real eigenvalues and is more in-line with the principles of
quantum mechanics and therefore allows one to pose a matrix $A$ (or a small
variant thereof) as a Hamiltonian in an Ising model which is the current
technology used in quantum annealing and simulated bifurcation machines. This
also allows one the ability to easily translate between quantum, digital, and
simulated annealing for bench marking purposes.

######  Lemma 2.

While considering a binary vector $x\in\\{0,1\\}^{N}$ one can reduce a QUBO
into a single matrix (which need not be traceless). That is one can reduce

$x^{t}Ax+B\cdot x\rightarrow x^{t}\tilde{A}x$

where

$\tilde{A}=A+\text{diag}(B)$ (3)

with diag$(B)$ being a diagonal matrix with

$diag(B)_{ij}=B_{i}\delta_{ij}$

###### Proof.

Since each $x_{i}$ is 0 or 1 we trivially have $x_{i}=x_{i}^{2}$. This reveals

$B\cdot x=x^{t}\text{diag}(B)x$

Which allows us to factor $x^{t}$ on the left and $x$ on the right

$\displaystyle x^{t}Ax+B\cdot x$ $\displaystyle=$ $\displaystyle
x^{t}Ax+x^{t}\text{diag}(B)x$ (4) $\displaystyle=$ $\displaystyle
x^{t}(A+\text{diag}(B))x$ $\displaystyle=$ $\displaystyle x^{t}\tilde{A}x$

∎

From here forward we shall simply refer to $\tilde{A}$ as $A$.

In quantum annealing, the vector $x_{i}\in\\{0,1\\}$ is exchanged for a vector
$z_{i}\in\\{-1,1\\}$ with the simple transformation

$z_{i}=2x_{i}-1$

This affects our matrix $A$ and also gives one an offset vector where

$\displaystyle x^{t}Ax=z^{t}Jz+C\cdot z+\text{const}$ (5)

The matrix $J=A/4$ and offset vector $C$ are calculated by a simple change of
variables. For the purposes of this article we will remain in the space
$x_{i}\in\\{0,1\\}$ and only mention that can transform when necessary.

##  3 The Constraint Matrix

Given our model

$\min_{x}x^{t}Ax$

with $A$ a real symmetric matrix of size $N\times N$ we wish to add a
constraint matrix $C$ to $A$ so that our vector $x$ has norm $\|x\|^{2}=m$; or
more simply the $L^{1}$ norm is $M$, $\|x\|_{1}=M$ for some $1\leq M\leq N$

Thus our new model becomes

$\min_{x}x^{t}(A+C)x$

If we pass this quadratic form to a solver (whether a simulated annealer or
quantum annealing computer or a simulated bifurcation machine) we can expect
that a “good” solver will produce a vector $x$ of the required size.

######  Theorem 3.

Let $A$ be a real symmetric matrix of size $N\times N$. Then the addition of
matrix

$C=\alpha J_{N}+\beta I_{N}$ (6)

to $A$ will produce a QUBO with optimal vector $x$ of size $\|x\|_{1}=M$ when

$M=\frac{-\beta}{2\alpha}$ (7)

In particular there is a line of solutions. The larger $\alpha$ the stronger
the pull toward $\|x\|_{1}=M$.

###### Proof.

Since we wish to minimize $x^{t}Ax$ where $|x|_{1}=M$ we consider the addition
of the matrix $C=\alpha J_{N}+\beta I_{N}$ where $J_{N}$ is the $N\times N$
matrix of all ones and $I_{N}$ is the $N\times N$ identity matrix.

This reduces our calculation to

$\min_{x}x^{t}(A+C)x\implies\min_{|x|_{1}=M}x^{t}Cx$ (8)

Since $J_{N}$ and $I_{N}$ yield straightforward multiplications we have

$x^{t}(\alpha J_{N}+\beta I_{N})x=\alpha M^{2}+\beta M$ (9)

Now we see the results directly. We choose $\alpha$ and $\beta$ to minimize
$x^{t}Cx$ by

$M=\frac{-\beta}{2\alpha}$ (10)

Noting that $\alpha>0$ gives a minimum and $\alpha<0$ gives a maximum.

Further more we can reparameterize this to be a single parameter constraint as

$C(\alpha)=\alpha(J_{N}-2MI_{N})$ (11)

We can also see this as

$C(\alpha)=\alpha\begin{bmatrix}1-2M&1&1&\dots&1\\\ 1&1-2M&1&\dots&1\\\
1&1&1-2M&\dots&1\\\ \vdots&\vdots&\vdots&\ddots&\vdots\\\
1&1&1&\dots&1-2M\end{bmatrix}$

∎

##  4 Two Examples

######  Example 4.

Let’s take a quick look at how to simulate this in a modern computing
language. We’ll produce a random symmetric matrix of size $N\times N$ and
require a solution of size $M$ where $M$ is significantly different from
$N/2$.

In the code on the author’s local computer, we will run a simulated annealer
which produces two outputs, (a) the cost $x^{t}Ax$, (b) the solution vector.

Algorithm 1 Picking M asets; Julia/Octave style

$A=randn(N,N)$

$A=(A+A^{t})/2$

$C(\alpha)=\alpha*(\text{ones}(N,N)-2M*\text{eye}(N))$

cost, solution = anneal($A+C(0.5)$)

sum(solution)

Using five hundred runs with parameters $N=30,M=8$ we have a histograms for
$\alpha\in\\{0,0.1,0.2,0.5,1,2,10\\}$. When $\alpha=0$ we have the
unconstrained model, where ones expects a solution of size $15\pm\sqrt{30}$ so
we expect between 10 and 20 as the global best solution. We are running our
simulated annealer cooling very quickly and only a few trials per degree so as
to show the efficacy of the constraint. We can see the pull toward 8 assets as
$\alpha$ increases.

(a) Unconstrained

(b) $\alpha=.1$

(c) $\alpha=0.2$

(d) $\alpha=0.5$

Figure 1: Histograms with increasing constraints

(a) $\alpha=1$

(b) $\alpha=2$

(c) $\alpha=10$

(d) $\alpha=0.5$

Figure 2: Histograms with increasing constraints

The final frame in 2 is a stacked histogram showing all 7 histograms together.
We see the density of selections increasing as we move toward 8. Looking all
the way back to the first figure in 1 we see that the global minimum is likely
larger than 15, most likely at 17. If we were to have a larger Gaussian
ensemble or one in which the global minimum is at or below $N/2$ then the
shift toward a smaller number of selected assets is even clearer. In this
case, if we had selected the number of assets to be 20, we would see the
rightward shift in the histograms more strongly.

######  Example 5.

For our second example, we’ll use a positive semi-definite matrix of size
$30\times 30$ and again apply constraints to select 8 assets.

It is important to note that the simulated annealer on the author’s computer
tends to avoid picking exactly zero assets. Thus in our first figure in 3 we
see a split between zero and one, even though we know the mathematically sound
answer is exactly zero.

(a) Unconstrained

(b) $\alpha=.1$

(c) $\alpha=0.2$

(d) $\alpha=0.5$

Figure 3: Histograms with increasing constraints

(a) $\alpha=1$

(b) $\alpha=2$

(c) $\alpha=10$

(d) $\alpha=0.5$

Figure 4: Histograms with increasing constraints

We see the efficacy of the constraint here much more prominently. Simulated
annealers are not necessary with such small scale QUBO as $30\times 30$ as a
brute force solution can be obtained in a matter of minutes. Additionally,
small shifts in small QUBO models don’t have such a large effect. If one were
to repeat these examples with $N=5000$ and $M=50$ the effects of $\alpha=0.1$
would appear much more clearly. Nonetheless even with $\alpha$ as small as
$0.2$ we begin seeing solutions with 8 assets. At $\alpha=0.5$ the lion’s
share are solutions with 8 assets (442 out of 500 in this particular numerical
experiment). And at $\alpha=2$ and above we select 8 assets without fail. Thus
is the last figure in 4 we see an intense density at 8 assets. With all other
assets smaller. This is consistent with the first example in which 8 assets
becomes the minimum number of assets chosen as the unconstrained solution has
greater than 8 assets, in this case the unconstrained solution has fewer.

## References

  * [CQ1] Cohen, J., Alexander, C., _Picking Efficient Portfolios from 3,171 US Common Stocks with New Quantum and Classical Solvers_ https://arxiv.org/abs/2011.01308
  * [CQ2] Cohen, J., Khan, A., Alexander, C., _Portfolio Optimization of 60 Stocks Using Classical and Quantum Algorithms_ https://arxiv.org/abs/2008.08669
  * [CQ3] Cohen, J., Khan, A., Alexander, C., _Portfolio Optimization of 40 Stocks Using the DWave Quantum Annealer_ https://arxiv.org/abs/2007.01430
  * [DA] Fujitsu Research Lab, _Quantum Future - Quantum Present_ https://sp.ts.fujitsu.com/dmsp/Publications/public/wp-da-overview-ww-en.pdf
  * [DwaveI] https://docs.dwavesys.com/docs/latest/c_gs_3.html
  * [EKU] Eren,Y., Küçükdemiral, İ, Üstoğlu1, İ, _Optimization in Renewable Energy Systems, Chapter 2 - Introduction to Optimization_ , "https://doi.org/10.1016/B978-0-08-101041-9.00002-8" "http://www.sciencedirect.com/science/article/pii/B9780081010419000028"
  * [GKD] Glover, F., Kochenberger, G., Du, Y., _A Tutorial on Formulating and Using QUBO models_ , 2019 https://arxiv.org/pdf/1811.11538.pdf
  * [MLO] Miguel, S., Lizaso, E., Orus, R., _Use Cases of Quantum Optimization for Finance_https://arxiv.org/abs/2010.01312
  * [SBM] Goto, H., Tatsumura, K., & Dixon, A.R.(2019) Combinatorial optimization by simulating adiabaticbifurcations in nonlinear Hamiltonian systems, Science Advances, 5(4), DOI:10.1126/sciadv.aav2372
  * [STB] Şeker O., Tanoumand N., Bodur M., _Digital Annealer for quadratic unconstrained binary optimization: a comparative performance analysis_https://arxiv.org/abs/2012.12264
  * [Wiki] https://en.wikipedia.org/wiki/Quadratic_unconstrained_binary_optimization