# Strong Coulomb Coupling Influences Ion and Neutral Temperatures in
Atmospheric Pressure Plasmas

M. D. Acciarri1, C. Moore2, S. D. Baalrud1 1 Department of Nuclear Engineering
and Radiological Sciences, University of Michigan, Ann Arbor, MI 48109, USA  2
Sandia National Laboratories, Albuquerque, NM 87185, USA

###### Abstract

Molecular dynamics simulations are used to model ion and neutral temperature
evolution in partially-ionized atmospheric pressure plasma at different
ionization fractions. Results show that ion-ion interactions are strongly
coupled at ionization fractions as low as $10^{-5}$ and that the temperature
evolution is influenced by effects associated with the strong coupling.
Specifically, disorder-induced heating is found to rapidly heat ions on a
timescale of the ion plasma period ($\sim 10$s ps) after an ionization pulse.
This is followed by the collisional relaxation of ions and neutrals, which
cools ions and heats neutrals on a longer ($\sim$ns) timescale. Slight heating
then occurs over a much longer ($\sim 100$s ns) timescale due to ion-neutral
three-body recombination. An analytic model of the temperature evolution is
developed that agrees with the simulation results. A conclusion is that strong
coupling effects are important in atmospheric pressure plasmas.

CAPP, Atmospheric Pressure Plasmas, Strong Correlations, Strongly Coupled

##  I Introduction

Atmospheric pressure plasmas, in particular cold atmospheric pressure plasmas
(CAPP), have been widely tested for numerous applications including the
inactivation of various pathogens in medicine [1], food industry [2],
agriculture [3], and water purification [4], among others, showing promising
results and an increasing interest in those industries. CAPP technology
presents several advantages over other plasma technologies including
operational simplicity, low running cost, and environmental friendliness,
since no vacuum chamber is needed and the reactor is frequently open [1]. The
increasing interest shown in atmospheric pressure plasmas and the lack of
understanding about the main mechanisms involved in the plasma dynamics,
transport of reactive species and chemistry motivates a better understanding
of the physics involved in order to improve the development of plasma sources
that promote reactions of interest and make devices more efficient [5, 6].
Here, we show an example of such a fundamental physics effect: ion-ion
interactions are strongly coupled and the ion and neutral temperatures are
influenced by an associated disorder-induced heating (DIH) process.

Strong coupling refers to interactions in which the average potential energy
of interacting particles [$\phi_{ss^{\prime}}(r=a_{ss^{\prime}})$, where
$a_{ss^{\prime}}$ is the average distance between particles of species $s$ and
$s^{\prime}$] exceeds their average kinetic energy ($k_{B}T_{ss^{\prime}}$),
i.e. $\Gamma_{ss^{\prime}}\gtrsim 1$, where $s$ and $s^{\prime}$ are the
species involved in the interaction and,

$\centering\Gamma_{ss^{\prime}}=\frac{\phi_{ss^{\prime}}(r=a_{ss^{\prime}})}{k_{B}T_{ss^{\prime}}}.\@add@centering$
(1)

For example, the Coulomb potential for ion-ion interactions is
$\phi_{ii}=(Ze)^{2}/4\pi\epsilon_{0}a_{ii}$ and the ion-ion Coulomb coupling
parameter is $\Gamma_{ii}=(Ze)^{2}/(4\pi\epsilon_{0}a_{ii}k_{B}T_{i})$. In
most CAPPs, ions are expected to be in equilibrium with the neutral gas near
room temperature, while electrons have a much higher temperature on the order
of eV or several eV. A consequence is that the hotter electrons are
characterized by a weakly coupled regime ($\Gamma_{ee}\ll 1$), while the much
cooler ions can be strongly coupled ($\Gamma_{ii}\gtrsim 1)$.

Strongly coupled plasmas are influenced by physical effects that fundamentally
differ from those governing weakly coupled plasmas. An example is disorder-
induced heating. This arises in strongly coupled systems when the potential
energy landscape changes in such a way that particles move to a lower
potential energy configuration, releasing kinetic energy in the form of heat.
One way that this can occur is through the ionization of a neutral gas. When
short-range atomic interactions change to long-range Coulomb interactions,
ions reconfigure to a more ordered state due to their mutual repulsion. The
more ordered state has a lower potential energy than the initial state
(immediately after ionization) because ions spread further apart from one
another. The decrease in potential energy in this reconfiguration is
compensated by an increase in kinetic energy. Disorder-induced heating is not
important in weakly coupled plasmas because the kinetic energy gained in the
reconfiguration is small compared to the initial kinetic energy. However, it
can be a dominant effect in strongly coupled plasmas. For example, DIH caused
by ionization has been measured and studied in detail in ultracold neutral
plasmas [7, 8, 9, 10].

Here, we show that ionization can cause DIH in atmospheric pressure plasmas as
well, and that it can be a dominant effect that determines both the ion and
neutral temperature evolution. Using Molecular Dynamics (MD) simulations, we
show that immediately after ionization ions have a coupling strength much
larger than one, with $\Gamma_{ii}$ ranging from approximately 5 to 25 for
simulations with ionization fractions ranging from 0.01 to 0.7. Ions then
rapidly heat over a timescale of approximately one ion plasma period
($\omega_{pi}^{-1}$, where $\omega_{pi}=\sqrt{e^{2}n_{i}/\epsilon_{o}m_{i}}$
is the ion plasma frequency) to a condition where $\Gamma_{ii}\approx 1$. This
typically corresponds to a $\sim$ps timescale. This can raise the ion
temperature several times, depending on the ionization fraction. After DIH,
ions thermally equilibrate with neutrals via elastic scattering, causing them
to cool and the neutral gas to heat. The cooling causes ions to return to a
more strongly coupled state. Depending on the ionization fraction, the neutral
heating can be substantial, leading to neutral temperatures that are several
times room temperature. Ion-neutral thermal equilibration typically takes
hundreds to thousands of ion plasma periods, corresponding to a $\sim$ns
timescale. In parallel with these comparatively fast processes, a weaker
heating of both ions and neutrals is observed over a much longer timescale
($\sim 100$s ns) due to ion-neutral three-body recombination. In addition to
the simulations, a model is developed to describe the main features of the ion
and neutral temperature evolution.

An implication of these findings is that common approaches to modeling
atmospheric pressure plasmas need to be reassessed. Most modeling techniques,
including Particle-In-Cell (PIC) methods and solutions of multi-fluid
equations, are based on solving a Boltzmann equation or approximations of it
(such as multi-fluid models obtained from taking moments of the Boltzmann
equation). However, the Boltzmann equation is valid only for weakly coupled
systems, such as dilute gases and plasmas. At a fundamental level, it does not
treat strong coupling effects that can influence dynamics of atmospheric
pressure plasmas. Disorder-induced heating is an example of this. If strong
correlations are not accounted for, the ion and neutral temperatures could be
mistakenly underestimated.

Although ion temperatures are difficult to measure in atmospheric pressure
plasmas, measurements have been made of neutral gas temperatures that are
elevated well above the ambient room temperature conditions [11, 12, 13, 14,
15, 16, 17]. Some explanations for this have been postulated. One is that
energy transfer from electron-neutral elastic collisions heats neutrals [18].
Although some heating through this mechanism is certain, it is not clear if it
happens at a fast enough rate or if the resulting ion temperature is correct
to explain the measurements. Other possible explanations for a more rapid
transfer of energy from electrons to neutral atoms/molecules include energy
transfer between vibrationally excited states in molecules [19, 20], energy
transfer from vibrational to translational states [19, 20, 15], electron-ion
recombination processes [16] followed by ion-neutral elastic collisions,
electron impact dissociation reactions [13] and quenching of electronically
excited molecules by oxygen atoms [14]. Although each is a plausible
mechanism, there hasn’t been a sufficiently detailed study to test if they
explain the experimental measurements. Here, we propose that DIH is another
mechanism that should be considered. We compare the predicted neutral gas
temperature after ions and neutrals have thermally equilibrated to an
experiment [18], and find that the predictions are close to the measured
temperature. These results are relevant for a number of plasma sources that
achieve high ionization fractions at atmospheric and elevated pressures,
including microdischarges [21], spark filaments [22, 23, 24], laser-produced
plasmas [25, 26], arc filaments [27] and the anode region of electric arcs
[28], among others.

In the next section the coupling parameter associated with ion-ion, ion-
neutral and neutral-neutral interactions is shown in a pressure–ionization
fraction space, predicting that at atmospheric pressure and over a wide range
of ionization fractions ions are expected to be strongly coupled. In section
III the MD simulation setup is described. This is followed by a discussion of
the ion and neutral temperature evolution in section IV. In section V a
theoretical model for predicting the maximum and equilibrium ion temperature
as a function of the ionization fraction and pressure is described. Finally,
section VI shows a comparison of the observed heating in experiments conducted
at atmospheric pressure and the prediction for the temperature using our
model.

##  II Coupling Parameter Space

Figure 1: Coupling parameter space for ion-ion, ion-neutral and neutral-
neutral interactions for an Ar gas at room temperature and different pressures
and ionization fractions.

For a partially ionized plasma with one gas species, neutral-neutral, ion-
neutral and ion-ion interactions are modeled using the Lennard-Jones, charge
induced dipole and Coulomb potentials respectively [29]

$\phi_{\textrm{LJ}}(r)=4\epsilon\left[\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^{6}\right],$
(2)
$\phi_{\textrm{ind}}(r)=-\frac{q^{2}}{8\pi\epsilon_{0}}\frac{\alpha_{\textrm{R}}a^{3}_{0}}{r^{4}},$
(3)

and

$\phi(r)=\frac{q^{2}}{4\pi\epsilon_{0}}\frac{1}{r},$ (4)

where $\epsilon=120k_{B}$, $\sigma=0.34$ nm, and $\alpha_{\textrm{R}}=11.08$
for Ar [29] and $a_{0}$ is the Bohr radius. Considering an Ar plasma at room
temperature and variable ionization fraction and pressure, the coupling
parameter associated with each interaction can be computed from equations
(1)–(4), using $a_{in}=(3/4\pi n_{in})^{1/3}$ where $n_{in}\approx
x_{i}n_{i}+x_{n}n_{n}$ to estimate the average interparticle spacing between
ions and neutrals. Here, $x_{i}=n_{i}/n$ and $x_{n}=n_{n}/n$ are the ion and
neutral concentrations, $n_{i}$ and $n_{n}$ are the ion and neutral densities,
and $n=n_{i}+n_{n}$ is the total density. The pressure and ionization fraction
at which the transition from a weakly to a strongly coupled regime occurs is
estimated from the limit $\Gamma=0.1$, which is the condition where the
Boltzmann equation is known to break down [30, 31].

As shown in figure 1, at small ionization fractions and pressures below
atmospheric pressure, none of the interactions are strongly coupled. However,
by increasing the ionization fraction or the pressure, it is possible to find
a strong ion-ion coupling region. In fact, the figure shows that ion-ion
interactions at many CAPP conditions are expected to be strongly coupled. If
the pressure is increased above $\approx$ 10 atm, a strong ion-neutral
coupling regime is reached. Further increases in pressure can lead to a strong
neutral-neutral coupling regime at around 1000 atm.

Since the ion-ion Coulomb coupling parameter exceeds 0.1 in many CAPPs
applications, the binary ion-ion collision picture and hence the Boltzmann
equation are not expected to apply. In this strongly coupled regime, many-body
collisions dominate the interactions between ions and a first-principles
approach is required. Hence, we used MD simulations, which provide a first-
principles solution of Newton’s equation of motion for the dynamics of $N$
interacting particles.

##  III Molecular Dynamics Simulations

Molecular dynamics simulations were carried out using the open-source software
LAMMPS [32]. Since electrons are much hotter than ions and are weakly coupled,
they are treated as a background non-interacting species when modeling ion and
neutral dynamics. Thus, they were not included in the simulation. This is
similar to the one-component plasma model [33], which is known to provide an
accurate description of ions in the presence of weakly coupled (comparatively
hot) electrons, such as in ultracold neutral plasmas [34]. Short (neutral-
neutral), medium (ion-neutral) and long (ion-ion) range interactions were
modeled using the potentials defined in equations (2)–(4).

Since the charge induced dipole potential is attractive, particles can
interact at short spatial scales. In order to avoid the rare occurrence of
close interactions that require a very short timestep to resolve, a repulsive
core term was added to the charge-induced dipole potential from expression (3)

$\centering\phi_{\textrm{ind}}(r)=\frac{q^{2}}{8\pi\epsilon_{0}}\frac{\alpha_{R}a^{3}_{0}}{r^{4}}\biggl{(}\frac{r_{\phi}^{8}}{r^{8}}-1\biggr{)},\@add@centering$
(5)

where $r_{\phi}$ is the radius at which the repulsive core acts. It was
desired to choose $r_{\phi}$ to be small enough to minimize the occurrence of
non-physical force values, but large enough to decrease the computational cost
due to the smaller timestep requirement. This is similar to what has been done
in MD simulations of ultracold neutral plasmas [35, 36]. The effect of the
repulsive core also depends on the simulation setup. For simulations conducted
at thermodynamic equilibrium, we found that the value of $r_{\phi}$ sets the
number of ions that attach to neutrals; as shown in the Appendix A. However,
this work concentrates on non-equilibrium simulations where $r_{\phi}$ is only
a numerical convergence parameter needed to avoid closely orbiting particles.
The reason why $r_{\phi}$ is not important in the non-equilibrium simulations
is that ion-neutral three-body recombination is slow compared to the DIH and
ion-neutral relaxation timescales that we focus on. For the non-equilibrium
simulations we chose $r_{\phi}=0.133a_{\textrm{in}}$ (here $a_{\textrm{in}}$
was estimated as $\approx a_{\textrm{nn}}$ for $x_{i}<0.5$ and $\approx
a_{\textrm{ii}}$ for $x_{i}>0.5$). This is the smallest value at which we can
avoid closely orbiting particles and non-physical force values; see Appendix A
and figure 2. A convergence test was conducted varying $r_{\phi}$ from
0.133$a_{\textrm{in}}$ to 0.5$a_{\textrm{in}}$ at an ionization fraction of
$x_{i}=0.5$ (figure 2) showing that the ion and neutral temperature at the end
of the simulation (after ion-neutral relaxation) is independent of $r_{\phi}$
when $r_{\phi}<0.25a_{\textrm{in}}$.

Figure 2: Ion and neutral final temperatures for different $r_{\phi}$ values.

In order to study the evolution of a non-equilibrium discharge, a neutral Ar
gas at room temperature and atmospheric pressure was simulated until
equilibrium was reached. This stage of the simulation was run with a Nosé-
Hoover thermostat (NVT ensemble) applied [37]. Then, a fraction of the
particles were instantly ionized and a NVE (microcanonical) simulation was run
including the ion-neutral, neutral-neutral and ion-ion interactions. This
simulation setup was repeated for different ionization fractions. The timestep
used was $5\times 10^{-4}\omega_{pi}$ and the NVE simulation was run until the
equilibrium was reached. The number of particles was varied such that the
minimum number of ions in the system was 2500 and the number of neutral atoms
was scaled in order to achieve the desired ionization fraction in each
simulation. The simulation domain size was scaled with the total number of
particles in order to maintain the desired gas density. The simulation domain
was a three-dimensional box with periodic boundary conditions. The ionization
fraction, total number of particles, time step and length of simulation used
for each simulation are specified in table 1.

Table 1: Parameters used for the Molecular Dynamics simulations. The time step value was $\Delta t\times\omega_{pi}=5\times 10^{-4}$ in all cases.  | Ionization
---
Fraction
($x_{i}$)
| Number of
---
Particles
| Plasma
---
Frequency
($\omega_{pi}$)
$\times$ $10^{12}$ (rad/s)
| Length of
---
Simulation
($t_{f}$ $\times$ $\omega_{pi}$)
0.01 | 250000 | 0.1046 | 150
0.1 | 25000 | 0.3308 | 500
0.3 | 10000 | 0.5731 | 650
0.5 | 10000 | 0.7398 | 800
0.7 | 10000 | 0.8754 | 1000

##  IV Molecular Dynamics Results

As shown in figure 3, the evolution of the ion temperature can be divided into
three stages. First, a rapid increase in the ion temperature was observed over
the first ion plasma period of the simulation. This is thought to be due to
disorder induced heating. This stage generates fluctuations in the ion
temperature that persist for several plasma periods; as shown in figures 3 and
7. Secondly, ion-neutral temperature relaxation was observed with a timescale
characterized by the ion-neutral collision frequency. Finally, a gradual
heating of both ions and neutrals is observed over a much longer timescale due
to three-body recombination of ions and neutrals (this is difficult to view in
figure 3, but is demonstrated more clearly below).

Figure 3: Evolution of the ion and neutral temperatures during the discharge
at atmospheric pressure for different ionization fractions.

###  IV.1 Disorder Induced Heating

Before the ionization pulse is applied, the initial distribution of positions
of neutral atoms corresponds to the equilibrium state of a gas interacting via
the short-range Lennard-Jones potential. Therefore, ionization creates many
ions separated by distances that are much smaller than the average distance
between ions $a_{ii}=(3/4\pi n_{i})^{1/3}$ where $n_{i}$ is the ion density.
Since the interaction between ions is governed by the Coulomb potential, which
is a long range potential, this leads to a large repulsion force between ion
pairs that brings ions apart, as illustrated in figure 4.

The separation of ions corresponds to the formation of a correlated state.
This can be quantified by the radial distribution function $[g(r)]$, which is
defined by setting $n_{o}g(r)4\pi r^{2}dr=N(r)$, where $N(r)$ is the total
number of particles in a spherical shell of radius $r$ and thickness $dr$
centered on a chosen particle. Here, $n_{o}=N/V$ is the average number of
particles ($N$) in a volume ($V$); i.e., it is the uniform background number
density. When $g(r)=1$ for all distances $r$, the system is in an uncorrelated
state. Figure 5 shows that just after ionization the ion positions correspond
to the weakly correlated state of the neutral gas they were formed from.
Correlations quickly develop as ions spread apart over the timescale of an ion
plasma period. This is indicated by the void of particles that forms from
distances $r=0$ to approximately $r=a_{ii}$. Such a void, which is sometimes
referred to a Coulomb hole in plasmas, is a characteristic property of a
strongly coupled system. Another characteristic is a peak near the average
nearest neighbor distance ($a_{ii}$), which is also observed at times later
than one ion plasma period.

Figure 4: After an instant ionization pulse, ions can be at distances much
smaller compared to $a_{ii}$ leading to large repulsion forces.

A simple estimate for the magnitude and timescale of DIH can be obtained by
considering the motion of a typical ion. In the limit that the initial
configuration is random in space, and the final distribution is perfectly
ordered (a lattice), ions will move an average distance of approximately
$a_{ii}/2$. The electrostatic potential change in moving from a position
$r_{1}=a_{ii}/2$ to $r_{2}=a_{ii}$ is
$\Delta\phi=\phi_{2}-\phi_{1}=-eZ/(4\pi\epsilon_{0}a_{ii})$. This will
decrease the total potential energy of the ions by a factor of $\Delta$U
$\approx-Z^{2}e^{2}/(4\pi\epsilon_{0}a_{ii})$. Since energy is conserved
during the NVE simulation, this leads to a corresponding increase in the ion
kinetic energy $\Delta$K $=-\Delta$U $\approx
Z^{2}e^{2}/(4\pi\epsilon_{0}a_{ii})$. For ionization fractions above
approximately $10^{-4}$, the initial value of $\Gamma_{ii}$ based on room
temperature (when the ionization pulse is applied) is considerably larger than
1, so the change in the ion kinetic energy due to DIH is much greater than the
kinetic energy before the pulse. Estimating the temperature after DIH as
$\frac{3}{2}k_{B}T\approx\Delta$K, the Coulomb coupling parameter from
equation (1) is $\Gamma_{ii}\approx 1.5$. In all cases simulated, DIH is
observed to increase the ion temperature until a critical value of
approximately $\Gamma_{ii}=1.9$ is reached; see figure 6. This is quite close
to the estimate of 1.5. The value of the temperature directly after DIH is the
maximum value observed, $T^{\max}_{i}$, as ions subsequently cool due to
collisional relaxation with neutrals.

Figure 5: Ion-ion radial distribution function $g(r)$ at different timesteps
during the simulation of a partially ionized Ar plasma with an ionization
fraction of $x_{i}=0.01$. The ion-ion radial distribution function corresponds
to an OCP at the same equilibrium $\Gamma_{ii}$ showing that the ions are
strongly coupled and the ion-ion interactions are not screened by the presence
of neutral atoms. In each simulation, the $g(r)$ was obtained after computing
an average over all the ions in the simulation and 100 consecutive time steps
that corresponded to a time window of $t\times\omega_{pi}=0.05$. An ensemble
average was also computed using data from 20 different simulations.

The timescale for DIH can be estimated from the time it takes an ion to move
from a position $r_{1}=a_{ii}/2$ to $r_{2}=a_{ii}$ due to the Coulomb
repulsion of a nearby ion. Considering Newton’s equation of motion,
$md^{2}x/dt^{2}=eE$ and approximating $d^{2}x/dt^{2}\sim\Delta x/\Delta
t^{2}\sim(a_{ii}/2)/\Delta t^{2}$, and $eE\sim
e^{2}/(4\pi\epsilon_{o}a_{ii}^{2})$, the characteristic timescale for this
process is the ion plasma period $\Delta t\approx\omega_{pi}^{-1}$. This also
agrees well with the simulations, where the maximum ion temperature is reached
approximately $1.5\omega_{pi}^{-1}$ after the ionization pulse; see figure 6.

Disorder induced heating has been previously observed in ultracold neutral
plasmas formed by photoionizing laser-cooled atoms [7]. Immediately after
ionization ions have a very small kinetic energy since the neutral gas they
are formed from was at millikelvin temperature. However, they have significant
excess potential energy since the ionization has changed the potential energy
landscape to that associated with long-range interactions. As the excess
potential energy is converted into kinetic energy the ions heat to a state
where $\Gamma_{ii}\approx 1$ over a timescale of $1\omega_{pi}^{-1}$, just as
described above. The process observed here is essentially the same. One
important difference is that the ionization pulse in an atmospheric pressure
plasma typically only partially ionizes the gas, whereas near total ionization
is common in an ultracold neutral plasma.

An implication is that DIH is only expected to be important if the ionization
fraction is high enough. Heating occurs only when the initial ion state (just
after ionization) satisfies $\Gamma_{ii}>1$. Otherwise, the kinetic energy
gained by DIH is small compared to the initial kinetic energy. Furthermore,
the basic mechanism isn’t expected to apply since the ions remain in a weakly
coupled disordered state even after ionization. Assuming ions are born from
neutral gas at room temperature, $\Gamma_{ii}>1$ requires that $x_{i}\gtrsim
10^{-4}$. Thus, at room temperature atmospheric pressure gas conditions, DIH
is expected to occur only if the ionization fraction is larger than one part
in ten thousand $x_{i}\gtrsim 10^{-4}$.

After DIH, ions overshoot their equilibrium positions leading to oscillations
of the Coulomb potential energy near the ion plasma frequency. Since the total
energy is conserved during the NVE simulation and the ions are strongly
coupled after the DIH, those oscillations translate to observable kinetic
energy fluctuations. This is shown in figure 7 from the simulation with an
ionization fraction of 0.01. Such fluctuations are only noticeable in the
strongly coupled regime, where the ion kinetic energy is comparable to the
Coulomb potential energy ($\Gamma\approx 1$). These fluctuations were observed
to damp over a period of time between 2$\omega^{-1}_{p}$ and
50$\omega^{-1}_{p}$ depending on the ionization fraction; see figure 3. The
presence of large fluctuations in the ion temperature after DIH has also been
observed in ultracold neutral plasmas experiments [7] and MD simulations [35].

Figure 6: Evolution of the Coulomb coupling parameter $\Gamma_{ii}$ during the
discharge at different ionization fractions. Figure 7: Ion kinetic energy
($K_{i}$), neutral kinetic energy ($K_{n}$) and total energy $E_{T}$ for a
discharge at atmospheric pressure and an ionization fraction of $x_{i}=0.01$.
Each energy is normalized by the number of ions ($N_{i}$), neutral atoms
($N_{n}$) and total number of particles $N$ respectively and the kinetic
energy at room temperature $T_{0}=300\;K$. Large fluctuations were observed in
the ion kinetic energy due to the oscillations in the Coulomb potential energy
after the DIH.

###  IV.2 Ion-Neutral Temperature Relaxation

After disorder-induced heating, ion and neutral temperatures equilibrate due
to ion-neutral collisions. This causes the ions to cool and the neutrals to
heat by an amount and at a rate that depends on the ionization fraction. Ion
cooling increases the ion-ion coupling strength, as shown in figure 6. This
leads to an ion-ion coupling strength that is larger than one,
$\Gamma_{ii}>1$, after the ion and neutral temperatures relax for all
ionization fractions simulated. It is also noteworthy that fast neutral
heating can be significant, especially when the ionization fraction is high.
Neutral heating is a common observation in atmospheric pressure plasma
experiments [16, 14, 11, 13].

Ion-neutral temperature relaxation was modeled using standard methods based on
the Boltzmann equation [38]. The cross section for ion-neutral interactions
was computed based on two-body collisions interacting through the charge-
induced dipole potential from equation (5). The Boltzmann-based approach is
expected to be valid since the ion-neutral interaction is in a weakly coupled
regime; see figure 1. The resulting temperature relaxation rate is described
by

$\frac{dT_{i}}{dt}=-\frac{3}{2}\nu_{in}(T_{i}-T_{n}),$ (6)
$\frac{dT_{n}}{dt}=\frac{3}{2}\nu_{ni}(T_{i}-T_{n}),$ (7)

and

$\nu_{ss^{\prime}}=\frac{4n_{s^{\prime}}\bar{v}_{ss^{\prime}}}{3}\int_{0}^{\infty}dg\;\;Q^{(1)}_{ss^{\prime}}(g)g^{5}e^{-g^{2}}.$
(8)

Here, $T_{i}$ and $T_{n}$ are ion and neutral temperatures respectively,
$\nu_{ss^{\prime}}$ is the energy transfer collision frequency between the
species $s$ and $s^{\prime}$, $g=u/\bar{v}_{ss^{\prime}}$ where $u$ is the
relative velocity and
$\bar{v}_{ss^{\prime}}^{2}=2k_{B}T_{s}/m_{s}+2k_{B}T_{s^{\prime}}/m_{s^{\prime}}$,
where $m_{s}$ and $m_{s^{\prime}}$ are both the Ar mass, and
$Q^{(1)}_{ss^{\prime}}$ is the momentum transfer cross section

$Q^{(1)}_{ss^{\prime}}=2\pi\int^{\infty}_{0}[1-\cos(\chi)]bdb,$ (9)

where

$\chi=\pi-2b\int^{\infty}_{r_{0}}\frac{dr/r^{2}}{\sqrt{1-\frac{b^{2}}{r^{2}}-\frac{2\phi_{in}(r)}{m_{ss^{\prime}}u^{2}}}}$
(10)

is the scattering angle. Here, $b$ is the impact parameter, $r_{0}$ is the
distance of closest approach obtained from the largest root of the denominator
in equation (10), $r$ is the radial distance, $m_{ss^{\prime}}$ is the reduced
mass, $u$ is the relative velocity and $\phi_{in}$ is the charge induced
dipole potential from equation (5) with $r_{\phi}=0.133a_{in}$ computed from
equation (9). The ion-neutral momentum transfer cross section obtained is
shown in figure 8.

The evolution of the ion temperature obtained using equations (6)–(8) is shown
in figure 9 for the discharge at an ionization fraction of 0.01. Here, the
initial temperature used in the model was taken from the value of the MD
simulations after DIH. The predictions show good agreement with the MD
simulations in both the equilibrium temperature as well as the relaxation
time. At other ionization fractions, the equilibrium temperature increases
with ionization fraction due to smaller neutral atom densities and due to the
larger ion temperature achieved due to the DIH. The equilibrium temperature
was higher than room temperature and the corresponding coupling parameter was
smaller compared to the the initial value of $\Gamma_{ii}$ with the exception
of the smallest ionization fraction of 0.01, where the large neutral density
led to an equilibrium temperature similar to room temperature; see figure 3.

Figure 8: Ion-neutral momentum transfer cross section at different relative
velocities, obtained for the charge induced dipole potential used in the
simulations. Figure 9: Evolution of the ion temperature for a discharge with
$x_{i}=0.01$ along with the model described in equations (6)–(8).

###  IV.3 Three-Body Recombination

As described in section III, $r_{\phi}$ was used as a numerical convergence
parameter in the non-equilibrium simulations and the temperature evolution
over the timescale of DIH and ion-neutral equilibration was independent of
$r_{\phi}$. However, as shown in the Appendix A, for $r_{\phi}/a_{in}<0.133$
neutral particles start to orbit ions, which forms bound states that can cause
heating over a much longer timescale than the DIH or ion-neutral relaxation.
In order to compare these time scales, we repeated the non-equilibrium
simulation on a longer time scale with different values of $r_{\phi}/a_{in}$.
Simulations at an ionization fraction of $x_{i}=0.01$ and $r_{\phi}/a_{in}$
values of $0.100$ and $0.090$ were carried out with a total simulation time of
$t\times\omega_{pi}=4000$ and compared with a simulation at the same
ionization fraction but using $r_{\phi}/a_{in}=0.133$. As shown in figure 10,
the simulations show a similar ion temperature evolution immediately after the
ionization pulse, including DIH followed by ion-neutral temperature relaxation
through collisions. However, at $t\times\omega_{pi}\approx 50$, bound states
start to form when $r_{\phi}/a_{in}<0.133$, which increases the ion
temperature when compared to the same simulation with $r_{\phi}/a_{in}=0.133$.
This is due to particles becoming trapped in the potential well of the charge-
induced dipole potential, which exchanges a reduction in the potential energy
with an increase in the kinetic energy. The increase of the fraction of ion-
neutral bound states over time is shown in figure 11 for both
$r_{\phi}/a_{in}$ values. The recombination fraction was determined from a
histogram of nearest neighbor positions, as described in appendix A. The
additional heating due to the ion-neutral three-body recombination causes the
temperature to increase by approximately 75 K and 125 K, reaching a value near
400 K and 425 K at $t\times\omega_{pi}=4000$ for $r_{\phi}/a_{in}$ values of
$0.100$ and $0.090$ respectively. This corresponds to a time at which $97\%$
of ions have at least one neutral atom orbiting them. It was observed that the
timescale of the ion-neutral three-body recombination is much longer than ion-
neutral temperature relaxation; $\sim 1000\omega_{pi}^{-1}$ corresponds to
$\sim 100$s ns.

Figure 10: Evolution of the temperature for a discharge with a ionization
fraction $x_{i}=0.01$ for different $r_{\phi}$ values.  Figure 11: Ion-neutral
three-body recombination rate for a discharge with a ionization fraction
$x_{i}=0.01$ for different $r_{\phi}$ values.

###  IV.4 Radial Distribution Function

The radial distribution gives information of how the density varies with the
distance from a reference particle and it can be used to look for signatures
of strong correlations. It also determines the non-ideal components of
thermodynamic parameters [39]. Figure 5 shows that the radial distribution for
the interactions between ions corresponds to that of the one-component plasma
(OCP) at the same equilibrium coupling parameter. Here, the ionization
fraction was 0.01. This result shows that the ion-ion interactions are in a
strongly coupled regime and the radial distribution function shows signatures
of these strong correlations. Since the radial distribution function
corresponds to an OCP at the same $\Gamma_{ii}$ the ion-ion interactions are
not screened by the presence of neutral atoms. This behaviour was observed for
all the ionization fractions simulated. The ion-ion radial distribution
function also exhibits a Coulomb hole and correlation peak, which are common
features of strongly coupled plasmas [40].

##  V Model for Maximum and Equilibrium Ion Temperature.

The maximum and equilibrium ion temperatures can be estimated using simple
energy conservation arguments. Since the increase in the ion temperature is
due to DIH, where the ions travel a distance $\approx\;a_{ii}$ over a plasma
period $\omega^{-1}_{p}$, the maximum ion temperature corresponds to that
which makes the ion-ion coupling parameter approximately unity. Utilizing the
simulation result that $\Gamma_{ii}\approx 1.91$ at the peak temperature and
equation (1), the maximum ion temperature is estimated to be

$T^{\max}_{i}=\frac{1}{1.91}\frac{Z^{2}e^{2}}{4\pi\epsilon_{0}k_{B}}\frac{1}{a_{ii}}$
(11)

where $a_{ii}=(3/4\pi x_{i}n)^{1/3}$ is the average interparticle spacing
between ions. The maximum ion temperature estimated with the equation (11)
shows good agreement with the results obtained from the MD simulations, as
shown in figure 12.

If we assume that the total kinetic energy is conserved after DIH until
thermodynamic equilibrium is reached, it is possible to write a simple energy
equation to obtain the equilibrium temperature $T_{\textrm{eq}}$, i.e., the
temperature after ion-neutral thermal equilibration,

$T_{\textrm{eq}}=x_{i}\;T^{\max}_{i}+(1-x_{i})\;T_{n}(t=0)$ (12)

where $T_{n}$ is the neutral atom temperature at the beginning of the
simulation, $x_{i}$ is the ionization fraction and $T^{\max}_{i}$ is the
maximum ion temperature obtained with equation (11). As shown in figure 12,
the values of $T_{\textrm{eq}}$ and $T^{\max}_{i}$ obtained using the model
show good agreement with MD simulations over a broad range of ionization
fractions. Furthermore, the expected coupling parameter with and without
accounting for the DIH can be calculated by using the correct equilibrium
temperature and the room temperature respectively, as shown in figure 13. The
values predicted for $\Gamma^{\textrm{eq}}_{ii}$ show good agreement with the
results from the MD simulations. It is noticeable how the coupling parameter
at equilibrium is smaller than what would be predicted using room temperature,
since the DIH increases the equilibrium temperature of the system at large
ionization fractions. However, at small ionization fractions the values for
$\Gamma_{ii}(T_{0},x_{i})$ and
$\Gamma^{\textrm{eq}}_{ii}(T_{\textrm{eq}},x_{i})$ match due to the large
neutral density compared to the ion density.

Figure 12: Variation of the maximum ion temperature and equilibrium
temperature with the ionization fraction from the MD simulations and the
model. Figure 13: The Coulomb coupling parameter at equilibrium is smaller
than the coupling parameter at room temperature without accounting for DIH. At
small ionization fractions $\Gamma^{\textrm{eq}}_{ii}(T^{\textrm{eq}},x_{i})$
converges to $\Gamma_{ii}(T_{0},x_{i})$ where $T_{0}$ is the room temperature.

##  VI Comparison with Experiment

Previous measurements have observed considerable neutral gas heating in ns
pulsed spark discharges at atmospheric pressure. For example van der Horst _et
al_ measured a neutral gas temperature of 750 K at at time 1 $\mu$s after
ignition in a N2/H2O mixture [18]. They also showed that the pressure due to
the discharge increased to 3 bar and the corresponding ionization fraction was
$x_{i}=0.16$. Using the model for the equilibrium temperature from equation
(12) and the pressure and ionization fraction measured in [18], we calculate
an expected ion-neutral equilibrium temperature of 765 K with a relaxation
time of 2.6 ns. This agrees quite closely with the experimental measurement.
The original reference suggested that the gas heating may be due to elastic
collisions between electrons and $N_{2}$. However, for the conditions of this
experiment, the energy transfer time due to elastic collisions is about 0.5
$\mu$s. This is quite a slow process where the relaxation time is near the
time at which the temperature was measured.

The heating mechanism proposed here occurs over the first 13 ps after ignition
for the pressure and ionization fraction measured in [18]. The subsequent ion-
neutral temperature relaxation occurs over 2.6 ns after the ignition. The
noticeable separation of timescales shows that a better time resolution for
experiments could help to identify the responsible heating mechanisms.

The results presented here suggest that the DIH and subsequent ion-neutral
relaxation should be considered as a possible mechanism for the observed
neutral heating.

##  VII Conclusions

This work demonstrates that ions are strongly coupled in atmospheric pressure
plasmas and that strong correlations influence ion and neutral temperature
evolution. First-principles MD simulations reveal that the ion and neutral
temperatures in CAPPs are influenced by an associated disorder induced heating
process. Disorder-induced heating is not important in weakly coupled plasmas,
but can be a dominant effect in strongly coupled plasmas. It causes ions to
rapidly heat to temperatures that can reach several times the background
neutral gas temperature. After DIH, ions thermally equilibrate with neutrals
causing them to cool and the neutrals to heat by an amount and at a rate that
depends on the ionization fraction. The cooling causes ions to return to a
more strongly coupled state. Due to DIH and ion-neutral temperature
relaxation, the final neutral gas temperature depends on the ionization degree
achieved at plasma formation. We show that at atmospheric pressure, DIH is
significant for ionization fractions larger than $x_{i}\gtrsim 10^{-4}$ and
thus, it could be important for experiments that involve spark discharges,
laser-produced plasmas, or other sources that achieve high ionization
fractions. In parallel with these comparatively fast processes, a weaker
heating of both ions and neutrals was observed over a much longer timescale
due to ion-neutral three-body recombination. A model is presented to describe
the main features of the ion and neutral temperature evolution, showing a good
agreement with the neutral gas temperature measured in an experiment at
atmospheric pressure [18]. Therefore, we show that DIH is an additional
heating mechanism that can occur in CAPP applications and must be accounted
for.

One implication of this work is that modeling techniques for CAPPs may need to
be reassessed. Current methods are largely based on solutions of a Boltzmann
equation, such as PIC, or moments of the Boltzmann equation, such as multi-
fluid models. These implicitly assume that the plasma is weakly coupled, since
the Boltzmann equation only treats weakly coupled gases and plasmas. This may
be limiting in the context of CAPPs, because fundamentally different physical
behaviors arise in strongly coupled fluids. For instance, DIH could not be
obtained by solving the Boltzmann equation, as in PIC simulations. Using such
models may lead to a significant underestimation of the ion and neutral
temperatures. This motivates the need for further development of theory and
computational models that can account for strong coupling effects in cold
atmospheric pressure plasmas.

##  VIII Acknowledgements

This work was supported in part by the US Department of Energy under award no.
DE-SC0022201, and in part by Sandia National Laboratories. Sandia National
Laboratories is a multi mission laboratory managed and operated by National
Technology and Engineering Solutions of Sandia, LLC., a wholly owned
subsidiary of Honeywell International, Inc., for the U.S. Department of
Energy’s National Nuclear Security Administration under Contract No. DE-
NA0003525. This article describes objective technical results and analysis.
Any subjective views or opinions that might be expressed in the paper do not
necessarily represent the views of the U.S. Department of Energy or the United
States Government.

##  Appendix A Equilibrium MD Simulations.

An analysis was realized varying the radius $r_{\phi}$ of the charge induced
dipole potential. For this analysis the simulation setup consisted of an NVT
(canonical) simulation using the Nosé–Hoover thermostat, where the temperature
was set to 293 K, over 1500 plasma periods. This was followed by an NVE
(microcanonical) simulation for 1500 plasma periods. The timestep used varied
with $r_{\phi}$ since smaller $r_{\phi}$ led to larger attractive forces and
smaller timesteps were required. The timesteps used varied from $5\times
10^{-4}\omega_{p}^{-1}$ to $1\times 10^{-4}\omega_{p}^{-1}$ where
$\omega_{p}^{-1}$ is the ion plasma period. The total number of particles was
10000 for all the simulations. Since the numerical charge induced dipole
potential has both a repulsive and an attractive part, it presented a
potential well, as shown in the figure 14. Different radius and, therefore,
minimum values of the charge induced dipole potential were explored. Once the
equilibrium was reached in the NVE simulation, the distribution of minimum
distances between ion-neutral and neutral-neutral pairs of particles was
computed in order to study the presence of bound states. These bound states
consisted of ion-neutral molecules with a spatial scale characterized by
$r_{\phi}$ and $r_{LJ}\approx 2^{1/6}\;\sigma$ which corresponds to the
minimum of the charge induced dipole and the Lennard Jones potentials
respectively.

Figure 14: Numerical charge induced dipole potential for different values of
$r_{\phi}$ compared to the charge induced dipole potential.

It was observed that for relatively small values of $r_{\phi}$ the formation
of bound states was more prominent. As it is shown in the distributions of
minimum distances between ion-neutral and neutral-neutral pairs of particles
after 1500 plasma periods in figure 15, for $r_{\phi}=0.046a_{in}$ and a
ionization fraction of 0.5, large peaks were observed corresponding to the
distances $r_{\phi}$ and $r_{LJ}$.

Figure 15: a) Ion-neutral and b) neutral-neutral minimum distance distribution
for $r_{\phi}=0.046a_{\textrm{in}}$ at an ionization fraction of 0.5.

As the radius of the charge induced dipole potential was increased, the
population of ion-neutral and neutral-neutral bound states decreased until it
was negligible for $r_{\phi}\approx 0.133a_{in}$ as shown in figure 16. The
observed behaviour suggests that in order to avoid the formation of bound
states an $r_{\phi}\approx 0.133a$ must be used.

Figure 16: Fraction of ion-neutral and neutral-neutral bound states for
different $r_{\phi}$ values.

## References

  * [1] Mária Domonkos, Petra Tichá, Jan Trejbal, and Pavel Demo.  Applications of cold atmospheric pressure plasma technology in medicine, agriculture and food industry.  Applied Sciences, 11(11), 2021.
  * [2] Patrick J. Cullen, James Lalor, Laurence Scally, Daniela Boehm, Vladimir Milosavljević, Paula Bourke, and Kevin Keener.  Translation of plasma technology from the lab to the food industry.  Plasma Processes and Polymers, 15(2):1700085, 2018.
  * [3] N. N. Misra, Oliver Schlüter, and P. J. Cullen, editors.  Cold plasma in food and agriculture: fundamentals and applications.  Elsevier/AP, Academic Press is an imprint of Elsevier, Amsterdam ; Boston, 2016.  OCLC: ocn934606423.
  * [4] I Adamovich, S Agarwal, E Ahedo, L L Alves, S Baalrud, N Babaeva, A Bogaerts, A Bourdon, P J Bruggeman, C Canal, E H Choi, S Coulombe, Z Donkó, D B Graves, S Hamaguchi, D Hegemann, M Hori, H-H Kim, G M W Kroesen, M J Kushner, A Laricchiuta, X Li, T E Magin, S Mededovic Thagard, V Miller, A B Murphy, G S Oehrlein, N Puac, R M Sankaran, S Samukawa, M Shiratani, M Šimek, N Tarasenko, K Terashima, E Thomas Jr, J Trieschmann, S Tsikata, M M Turner, I J van der Walt, M C M van de Sanden, and T von Woedtke.  The 2022 plasma roadmap: low temperature plasma science and technology.  Journal of Physics D: Applied Physics, 55(37):373001, jul 2022.
  * [5] E C Neyts and A Bogaerts.  Understanding plasma catalysis through modelling and simulation—a review.  Journal of Physics D: Applied Physics, 47(22):224010, may 2014.
  * [6] Annemie Bogaerts, Xin Tu, J Christopher Whitehead, Gabriele Centi, Leon Lefferts, Olivier Guaitella, Federico Azzolina-Jury, Hyun-Ha Kim, Anthony B Murphy, William F Schneider, Tomohiro Nozaki, Jason C Hicks, Antoine Rousseau, Frederic Thevenet, Ahmed Khacef, and Maria Carreon.  The 2020 plasma catalysis roadmap.  Journal of Physics D: Applied Physics, 53(44):443001, aug 2020.
  * [7] Thomas C. Killian and Steven L. Rolston.  Ultracold neutral plasmas.  Physics Today, 63(3):46–51, 2010.
  * [8] T. Pohl, T. Pattard, and J. M. Rost.  Relaxation to nonequilibrium in expanding ultracold neutral plasmas.  Phys. Rev. Lett., 94:205003, May 2005.
  * [9] D.O. Gericke and M.S. Murillo.  Disorder-induced heating of ultracold plasmas.  Contributions to Plasma Physics, 43(5-6):298–301, 2003.
  * [10] S. G. Kuzmin and T. M. O’Neil.  Numerical simulation of ultracold plasmas.  Physics of Plasmas, 9(9):3743–3751, 2002.
  * [11] David Z Pai, Deanna A Lacoste, and Christophe O Laux.  Nanosecond repetitively pulsed discharges in air at atmospheric pressure—the spark regime.  Plasma Sources Science and Technology, 19(6):065015, nov 2010.
  * [12] D L Rusterholtz, D A Lacoste, G D Stancu, D Z Pai, and C O Laux.  Ultrafast heating and oxygen dissociation in atmospheric pressure air by nanosecond repetitively pulsed discharges.  Journal of Physics D: Applied Physics, 46(46):464010, oct 2013.
  * [13] N A Popov.  Fast gas heating in a nitrogen–oxygen discharge plasma: I. kinetic mechanism.  Journal of Physics D: Applied Physics, 44(28):285201, jun 2011.
  * [14] N A Popov.  Pulsed nanosecond discharge in air at high specific deposited energy: fast gas heating and active particle production.  Plasma Sources Science and Technology, 25(4):044003, may 2016.
  * [15] Ryo Ono and Tetsuji Oda.  Measurement of gas temperature and OH density in the afterglow of pulsed positive corona discharge.  Journal of Physics D: Applied Physics, 41(3):035204, jan 2008.
  * [16] E I Mintoussov, S J Pendleton, F G Gerbault, N A Popov, and S M Starikovskaia.  Fast gas heating in nitrogen–oxygen discharge plasma: II. energy exchange in the afterglow of a volume nanosecond discharge at moderate pressures.  Journal of Physics D: Applied Physics, 44(28):285202, jun 2011.
  * [17] S Hofmann, A F H van Gessel, T Verreycken, and P Bruggeman.  Power dissipation, gas temperatures and electron densities of cold atmospheric pressure helium and argon RF plasma jets.  Plasma Sources Science and Technology, 20(6):065010, nov 2011.
  * [18] R M van der Horst, T Verreycken, E M van Veldhuizen, and P J Bruggeman.  Time-resolved optical emission spectroscopy of nanosecond pulsed discharges in atmospheric-pressure n2 and n2/h2o mixtures.  Journal of Physics D: Applied Physics, 45(34):345201, aug 2012.
  * [19] Igor V. Adamovich.  Three-dimensional analytic model of vibrational energy transfer in molecule-molecule collisions.  AIAA Journal, 39(10):1916–1925, 2001.
  * [20] Vasco Guerra, Antonio Tejero del Caz, Carlos D Pintassilgo, and Luís L Alves.  Modelling n2-o2 plasmas: volume and surface kinetics.  Plasma Sources Science and Technology, 28(7):073001, jul 2019.
  * [21] Lifang Dong, Junxia Ran, and Zhiguo Mao.  Direct measurement of electron density in microdischarge at atmospheric pressure by stark broadening.  Applied Physics Letters, 86(16):161501, 2005.
  * [22] E V Parkevich, M A Medvedev, G V Ivanenkov, A I Khirianova, A S Selyukov, A V Agafonov, Ph A Korneev, S Y Gus’kov, and A R Mingaleev.  Fast fine-scale spark filamentation and its effect on the spark resistance.  Plasma Sources Science and Technology, 28(9):095003, sep 2019.
  * [23] N Minesi, S Stepanyan, P Mariotto, G D Stancu, and C O Laux.  Fully ionized nanosecond discharges in air: the thermal spark.  Plasma Sources Science and Technology, 29(8):085003, aug 2020.
  * [24] A. Bataller, S. Putterman, S. Pree, and J. Koulakis.  Observation of shell structure, electronic screening, and energetic limiting in sparks.  Phys. Rev. Lett., 117:085001, Aug 2016.
  * [25] A. Bataller, G. R. Plateau, B. Kappus, and S. Putterman.  Blackbody emission from laser breakdown in high-pressure gases.  Phys. Rev. Lett., 113:075001, Aug 2014.
  * [26] Alexander Bataller, Alexandra Latshaw, John P. Koulakis, and Seth Putterman.  Dynamics of strongly coupled two-component plasma via ultrafast spectroscopy.  Opt. Lett., 44(23):5832–5835, Dec 2019.
  * [27] T Verreycken, R M van der Horst, A H F M Baede, E M Van Veldhuizen, and P J Bruggeman.  Time and spatially resolved LIF of OH in a plasma filament in atmospheric pressure he-h2o.  Journal of Physics D: Applied Physics, 45(4):045205, jan 2012.
  * [28] J Heberlein, J Mentel, and E Pfender.  The anode region of electric arcs: a survey.  Journal of Physics D: Applied Physics, 43(2):023001, dec 2009.
  * [29] Michael A. Lieberman and Allan J. Lichtenberg.  Atomic Collisions, chapter 3, pages 43–85.  John Wiley ‘I&’ Sons, Ltd, 2005.
  * [30] Brett Scheiner and Scott D. Baalrud.  Mean force kinetic theory applied to self-diffusion in supercritical lennard-jones fluids.  The Journal of Chemical Physics, 152(17):174102, 2020.
  * [31] Jérôme Daligault.  Diffusion in ionic mixtures across coupling regimes.  Phys. Rev. Lett., 108:225004, May 2012.
  * [32] 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, and S. J. Plimpton.  LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales.  Comp. Phys. Comm., 271:108171, 2022.
  * [33] Marc Baus and Jean-Pierre Hansen.  Statistical mechanics of simple coulomb systems.  Physics Reports, 59(1):1–94, 1980.
  * [34] T.C. Killian, T. Pattard, T. Pohl, and J.M. Rost.  Ultracold neutral plasmas.  Physics Reports, 449(4):77–130, 2007.
  * [35] S. G. Kuzmin and T. M. O’Neil.  Numerical simulation of ultracold plasmas: How rapid intrinsic heating limits the development of correlation.  Phys. Rev. Lett., 88:065003, Jan 2002.
  * [36] Sanat Kumar Tiwari, Nathaniel R. Shaffer, and Scott D. Baalrud.  Thermodynamic state variables in quasiequilibrium ultracold neutral plasma.  Phys. Rev. E, 95:043204, Apr 2017.
  * [37] Daan Frenkel and Berend Smit.  Chapter 6 - molecular dynamics in various ensembles.  In Daan Frenkel and Berend Smit, editors, Understanding Molecular Simulation (Second Edition), pages 139–163. Academic Press, San Diego, second edition edition, 2002.
  * [38] J. H. Ferziger and H. G. Kaper.  Mathematical theory of transport processes in gases.  American Journal of Physics, 41(4):601–603, 1973.
  * [39] Jean-Pierre Hansen and Ian R. McDonald.  Chapter 2 - statistical mechanics.  In Jean-Pierre Hansen and Ian R. McDonald, editors, Theory of Simple Liquids (Fourth Edition), pages 13–59. Academic Press, Oxford, fourth edition edition, 2013.
  * [40] T. Ott, M. Bonitz, L. G. Stanton, and M. S. Murillo.  Coupling strength in coulomb and yukawa one-component plasmas.  Physics of Plasmas, 21(11):113704, 2014.


# Simple and realistic models for efficient epidemic control on multiplex
networks

Minsuk Kim  Soon-Hyung Yook<EMAIL_ADDRESS>Department of Physics and
Research Institute for Basic Sciences, Kyung Hee University, Seoul 130-701,
Korea

## summary paragraph

The real epidemic spreading has generally two different types of transmission
routes. One is the random anonymous infection [1] and the other one is the
transmission through regular and fixed contacts [2]. If the infectious disease
has high mortality and there is no available vaccine or medicine, then many
health authorities rely on the non-pharmaceutical interventions (NPIs) along
the traceable fixed contacts, such as isolation of the infected [3]. In our
study, such realistic situations are implemented by the susceptible-infected-
recovered (SIR) model [4] with isolation on multiplex networks [5]. The
multiplex networks are composed of random interaction layer and fixed
interaction layer. We quantitatively compare the efficiency of two isolation
protocols imposed on the SIR dynamics. One of them is the most popular
protocol adopted by many health organizations over the globe. From the
numerical simulations we find that the isolation of the second nearest
neighbors of the hospitalized individuals on the layer for regular and fixed
contacts significantly reduces both the final epidemic size and the number of
the isolated per unit time. Our finding suggests a better NPI for any type of
epidemic even though the contact tracing is only partially available.

## Introduction

An outbreak of a new disease, such as the bubonic plague pandemic in the 14th
century [6], the 1918 influenza pandemic [7], and the recent outbreak of
severe acute respiratory syndrome [8], has been a large threat throughout
human history. Despite great advances in medical science and pharmacology,
immediate use of an effective vaccine or antiviral drug is not always possible
when new infectious diseases emerge. For example, due to the absence of
vaccines or antiviral drugs for new severe acute respiratory coronavirus 2
(SARS-CoV-2) during the early stage of the coronavirus disease 2019 (COVID-19)
pandemic, more than 172 million people have been infected and has caused more
than 3.7 million deaths until March 2021 [9]. Thus, finding an efficient NPI
is crucial to mitigate the pandemic situation for new emerging diseases.

The best NPI for a new disease is a perfect lockdown, under which all
individuals are strictly isolated. For example, during the early stage of the
COVID-19 pandemic, strict lockdown measures had been successfully applied in
mainland China and many European countries [10, 11]. However, the strict
lockdown policy is not sustainable if the pandemic period continues long
enough to cause a severe recession of economic activity and to increase social
fatigue [12, 13]. Thus, it is necessary to find NPIs that minimize the impact
on the social and economic systems. The efficacy of various NPIs has been
intensively studied based on real data and theoretical models to alleviate the
recent pandemic situations [14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 1, 25,
26, 27, 28].

Among the various NPIs, the quarantine of the infected individuals and their
contacts is one of the most intuitive measures and commonly shared by many
health authorities over the world [3]. Thus, the isolation of the infected
[28, 23] and tracing the contacts [29, 30, 14] are two important factors.
However, if the infection from asymptomatic and pre-symptomatic patients are
potential transmission routes like the COVID-19 case [31, 32], then finding
the contacts with such patients are not trivial. Furthermore, when airborne
transmission is another important route of spreading [33], the tracing becomes
much harder due to the random anonymous contacts through the publicly opened
environments [1]. In this study, such random anonymous transmission is
implemented by the double-layered multiplex networks (DLMNs) [5]. To model the
situation with the pre-symptomatic and asymptomatic transmissions, we assume
that individuals in our models have one of the following disease states,
susceptible ($S$), infected ($I$), and recovered ($R$) [4]. Furthermore, the
contact tracing probability and isolation states are introduced to account for
more realistic situations in our model. As we will show we first model the
most popular protocol for NPI adopted by many health authorities, and also
introduce a reinforced protocol model. From the quantitative comparison of the
two models, we suggest a simple and efficient NPI strategy for epidemic
controls of any emerging infectious disease by using only the known topology
of the fixed interaction layer.

## Model

Individuals are denoted by nodes and the interactions between them are
represented by links in the DLMN. Let $F$ and $W$ be the two layers in the
DLMN (see Fig. 1(a)). On $F$ each node is connected with randomly chosen $k$
nodes drawn from a given degree distribution $P_{F}(k)$. The topology of
networks on $F$ does not change in time. At the same time, each node interacts
with $k^{\prime}(t)$ random nodes on $W$, where $k^{\prime}(t)$ is drawn from
another degree distribution $P_{W}(k^{\prime})$ at each time $t$. Thus, the
interaction topology on $W$ changes at each $t$. Under severe epidemic
situations, the government tries to cordon people off public facilities, and
each individual refrains from social activities. Thus, the number of contacts
of each individual is significantly restricted and homogeneous. To generate
such homogeneous interaction structures, we use the Poisson distribution for
both $P_{F}(k)$ and $P_{W}(k^{\prime})$ [5, 34, 35] (see ’Underlying networks’
in Methods).

The state of each node at time $t$ in the DLMN is described by a two-component
variable ${\bm{\sigma}}=(\sigma_{1},\sigma_{2})$. $\sigma_{1}$ has one of the
three disease states: $S$, $I$, and $R$. $\sigma_{2}$ denotes the state for
isolation measure. When the disease causes severe symptoms, two isolated
states are possible: i) self-isolation at home when an individual feels mild
symptoms or recognizes a suspicious contact but has not been confirmed yet,
and ii) hospitalization by the health authority when the patient is confirmed
to be infected. If an individual is not isolated then it is in the unisolated
state. Thus, $\sigma_{2}$ can be one of the following states: self-isolated
($X$), hospitalized ($H$), and unisolated ($U$) (see Fig. 1(b)). Since we
cannot trace the contacts on $W$ due to the random anonymity, the self-
isolation for the suspicious contacts can be applied only to $F$.

Depending on the range of the self-isolation, we introduce two intervention
protocols, the basic isolation protocol (BIP) and the reinforced isolation
protocol (RIP). Under the BIP only the confirmed patients and the one who has
direct contact with the confirmed patient are isolated. This is the most
popular quarantine protocol adopted by many health authorities [36]. However,
the fraction of the household or workplace infection cannot be ignored in some
infectious diseases, for example, the household infection is more than 15% for
COVID-19 [2]. Such household and workplace infection can be caused by a self-
isolated node. Thus controlling such local contacts is another important
factor to mitigate the transmission. For the preemptive protection of the
susceptible, in the RIP model if a node is hospitalized, then its first and
also second nearest neighbors on $F$ are isolated.

In our models, each infected node transmits the disease to the connected
susceptible nodes with the probability $\beta_{F}$ ($\beta_{W}$) on $F$ ($W$).
With the probability $\theta_{X}$ ($\theta_{H}$) the nodes are self-isolated
(hospitalized) for the isolation (hospitalization) period $t_{X}$ ($t_{H}$).
The infected nodes are recovered after the recovery time $t_{R}$ (see ’Basic
Isolation Protocol’ and ’Reinforced Isolation Protocol’ in Methods). We set
the size of each layer as $N=100,000$. Since we are interested in the control
of the severe epidemic outbreak, we use $\beta_{F}=\beta_{W}\equiv\beta(=0.2)$
and $t_{R}=6$ to guarantee that the whole system becomes infected without any
intervention. In our model, the strength of the intervention measures is
controlled by four parameters, $\theta_{H}$, $t_{H}$, $\theta_{X}$, and
$t_{X}$. For simplicity, we assume that
$\theta_{H}=\theta_{X}\equiv\theta^{*}$ and $t_{H}=t_{X}\equiv t^{*}$. Thus,
we use only two control parameters $\theta^{*}$ and $t^{*}$ in the following
simulations.

## Results

### The fraction of nodes in each state

Let $\rho_{m}(t)$ ($m\in\\{S,I,R\\}$) be the fraction of nodes whose disease
states is $\sigma_{1}=m$ at time $t$, regardless of $\sigma_{2}$. If
$m\in\\{U,H,X\\}$, then $\rho_{m}(t)$ represents the fraction of nodes with
$\sigma_{2}=m$. The peak of $\rho_{m}(t)$ for each state $m$ is denoted by
$\rho^{peak}_{m}$. In Figs. 2(a)-(d) we show $\\{\rho_{m}(t)\\}$’s under the
BIP for the various values of parameters $(\theta^{*},t^{*})$. For small
$\theta^{*}(=0.3)$, we find that $\rho^{peak}_{I}(t)>0.9$ followed by
$\rho^{peak}_{H}(\gtrsim 0.8)$, regardless of $t^{*}$ (Figs. 2(a), (b)).
$\rho_{S}(t)$ rapidly decreases and reaches $\rho_{S}\approx 0$ for $t>11$.
The value of $\rho^{peak}_{X}$ ($<0.2$) is relatively small. Thus,
$\rho_{R}(t\rightarrow\infty)\simeq 1$ when $\theta^{*}$ is small. On the
other hand, as $\theta^{*}$ increases, $\rho^{peak}_{I}$ is drastically
suppressed as well as $\rho^{peak}_{H}$ (Figs. 2(c), (d)). As a result,
$\rho_{R}(t\rightarrow\infty)$ is reduced to
$\rho_{R}(t\rightarrow\infty)\simeq 0.6\sim 0.7$ for $\theta^{*}=0.9$. For
both values of $\theta^{*}$ displayed in Fig. 2, $t^{*}$ only affects the
behavior of $\rho_{H}$ and $\rho_{X}$ (the population of the isolated nodes),
i.e., $\rho_{H}^{peak}(t;t^{*}=4)<\rho_{H}^{peak}(t;t^{*}=10)$ and
$\rho_{X}^{peak}(t;t^{*}=4)<\rho_{X}^{peak}(t;t^{*}=10)$. Note that
$\rho^{peak}_{H}$ is comparable with $\rho^{peak}_{I}$ for all values of
$(\theta^{*},t^{*})$, and $\rho_{H}(t)$ becomes wider as $t^{*}$ increases.
This means that the hospitalized period becomes longer without any significant
change in the final epidemic size, $\rho_{R}(t\rightarrow\infty)$, as $t^{*}$
increases for all $\theta^{*}$. Thus, increasing $t^{*}$ without the
improvement of traceability causes an overload on the medical system by making
patients be hospitalized for a longer period.

Figs. 2(e)-(h) show $\\{\rho_{m}(t)\\}$’s for the RIP model. When
$\theta^{*}\lesssim 0.3$, $\\{\rho_{m}(t)\\}$’s for the RIP model show almost
the similar behavior with those for the BIP, but
$\rho_{S}(t\rightarrow\infty)$ for the RIP is slightly larger than that for
the BIP. Since additional nodes are self-isolated under the RIP,
$\rho_{X}^{peak}$ increases compared with that for the BIP with the same
$(\theta^{*}$, $t^{*}$). However, we find that $\rho_{H}$ for the RIP becomes
much smaller than that for the BIP. This effect becomes more drastic for
$\theta^{*}>0.3$. For example, $\rho_{R}(t\rightarrow\infty)$ and
$\rho_{I}^{peak}(t)$ significantly decrease to
$\rho_{R}(t\rightarrow\infty)=0.5\sim 0.6$ and $\rho_{I}^{peak}\approx 0.1$
for the RIP with $\theta^{*}=0.9$. $\left\\{\rho_{m}(t)\right\\}$’s for other
values of $\left(\theta^{*},t^{*}\right)$ are presented in the Extended Data
Fig. 2 and Extended Data Fig. 3. The results indicate that the collapse of the
medical systems can be avoided under the RIP if we trace the contacts with
sufficiently high accuracy. In addition, we find $\rho_{I}$, $\rho_{H}$, and
$\rho_{X}$ oscillate with decreasing amplitude under the RIP as $\theta^{*}$
increases. This suggests that even though there is a rapid decrease in
$\rho_{I}(t)$ after its first peak when $\theta^{*}$ is sufficiently large, it
is still possible to be followed by successive multiple peaks of $\rho_{I}$
(see also Extended Data Fig. 8).

### The effective reproduction number

To quantify the efficacy of intervention measures, we estimate the
instantaneous effective reproduction number, $R_{e}(t)$, at $t$. For a
practical purpose, we define $R_{e}(t)$ as

$R_{e}(t)=\frac{N^{new}_{I}(t)}{N_{I}(t-1)},$ (1)

where $N^{new}_{I}(t)$ is the number of new infected nodes at $t$ and
$N_{I}(t)$ is the number of infected nodes at $t$ [37, 38]. Thus $R_{e}(t)$
represents a metric to quantify how many nodes are newly infected by the
existing infected nodes at each $t$.

In Figs. 3(a),(b), we show $R_{e}(t)$ for the BIP and RIP with
$\theta^{*}=0.9$ and $t^{*}=2\sim 12$. The dashed line denotes $R_{e}(t)$
without intervention. As shown in Fig. 3(a), $R_{e}(t)$ for the BIP rapidly
decreases when $t\lesssim 4$ and shows a plateau followed by another rapid
drop, regardless of $t^{*}$. $R_{e}(t)=1$ at $t\approx 11$ and $R_{e}(t)<1$
for $t>11$. When $t>15$ $R_{e}(t)$ approaches to $R_{e}(t)\approx 0$. On the
other hand, $R_{e}(t)$ for $\theta^{*}=0.9$ under the RIP decreases more
drastically and $R_{e}(t)<1$ for $t\gtrsim 5$ as shown in Fig. 3(b). When
$t>20$, $R_{e}(t)$ oscillates with decreasing amplitudes and approaches
$R_{e}\approx 0$ under the RIP.

To investigate how $\theta^{*}$ affects the epidemic spreading, we also
measure $R_{e}(t)$’s for various $\theta^{*}$ when $t^{*}$ is fixed. In Figs.
3(c),(d), as an example, we display $R_{e}(t)$’s for $t^{*}=10$. Since
$\theta^{*}$ denotes traceability, $R_{e}(t)$ should decrease as $\theta^{*}$
increases for both protocols as shown in Figs. 3(c), (d). Note that when
$\theta^{*}<0.3$, the difference between the BIP and RIP is not noticeable.
However, if $\theta^{*}>0.3$, then $R_{e}(t)$ for the RIP becomes much smaller
than those for the BIP. From the data in Fig. 3, we find that increasing
$\theta^{*}$ is more important than increasing $t^{*}$ (See also Extended Data
Figs. 4 and 5).

The rapid drop of $R_{e}(t)$ under both protocols has two different origins
depending on $\theta^{*}$. For $\theta^{*}<0.3$ due to the large infection of
the early stage, there does not remain a sufficient number of susceptible
nodes for $t>11$ (see Figs. 2(a), (b), (e), (f)). On the other hand, if
$\theta^{*}>0.3$ then a significant amount of the susceptible is self-
isolated, which protects the susceptible nodes before contact with the
patients for $t>10$ (see Figs. 2(c), (d), (g), (h)).

### The final epidemic size and the number of isolated nodes per unit time

For a direct comparison between the two intervention protocols, we measure the
difference, $\Delta_{R}$, of the final epidemic sizes between the BIP and RIP
(see ‘Final epidemic size’ in Methods). Here, if $\Delta_{R}>0$ then
$\rho_{R}$ for the BIP is larger than that for the RIP. The data in Fig. 4(a)
clearly shows that $\Delta_{R}$ rarely depends on $t^{*}$. However,
$\Delta_{R}$ strongly depends on $\theta^{*}$. $\Delta_{R}\leq 0.05$ for all
$t^{*}$ when $\theta^{*}<0.3$, while $\Delta_{R}>0.1$ for $\theta^{*}\gtrsim
0.4$ and $\Delta_{R}$ increases as $\theta^{*}$ increases. This means that the
RIP significantly reduces the final epidemic size compared to the BIP when
$\theta^{*}\geq 0.3$. For the maximal traceability, $\theta^{*}=1$, we find
that $\rho_{R}(t\rightarrow\infty)$ for the RIP is reduced by 67% compared to
that for the BIP. This corresponds to the 100% increase of
$\rho_{S}(t\rightarrow\infty)$ under the RIP compared to the BIP (see Extended
Data Fig. 6). Thus, isolation of the possible suspicious contacts in advance
by applying the RIP significantly reduces the final epidemic size when
$\theta^{*}\geq 0.3$.

In epidemic control, reducing the number of isolated individuals at each time
step becomes another crucial factor to minimize social and economic recession.
For the quantitative analysis, we define the difference in the fractions of
the isolated nodes per unit time, $\Delta_{XH}$ (see ’Fraction of isolated
nodes per unit time’ in Methods). By definition, if $\Delta_{XH}>0$ then more
nodes are isolated under the BIP than the RIP. As shown in Fig. 4(b), for all
values of $\theta^{*}$, we find that $\Delta_{XH}\approx 0$ for $t^{*}\leq 6$.
However, we find that $\Delta_{XH}>0$ when $t^{*}>6$ and $\Delta_{XH}$
increases as $t^{*}$ increases. Thus, $t^{*}$ affects only $\rho_{X}$ and
$\rho_{H}$ per unit time for both models. Note that, even though
$\Delta_{XH}\leq 0$ for $t^{*}\leq 6$, $\Delta_{R}$ increases with
$\theta^{*}$ and $\Delta_{R}\geq 0$ as shown in Fig. 4(a)(see also Fig. 2).
Therefore, the RIP more effectively controls the disease spreading through the
preemptive isolation of suspicious contacts with fewer isolated nodes per unit
time than the BIP.

## Discussion

In summary, we model the NPI adopted by many health authorities over the
world, and introduce a model for reinforced NPI. In these models the state of
each individual is characterized by three disease states with additional
isolation states. Two different types of transmission routes observed in real
world are implemented by the multiplex networks. By using numerical
simulations, we compare the efficacy of the two models, BIP and RIP models,
and find that the RIP controls the spreading of disease more efficiently by
reducing both the epidemic size and the average number of isolated individuals
per unit time, despite its simplicity. Especially, when the traceability is
maximal, the final fraction of the susceptible nodes under the RIP increases
by almost 100% (almost doubled) compared to that under the BIP. This indicates
that the RIP significantly and efficiently protect the susceptible nodes
through the preemptive isolation of the possible contacts. Furthermore, since
we do not assume any characteristic property of a specific disease, we expect
that the suggested models can be used as a general framework for modeling
disease control for any real disease outbreak.

## References

  * [1] Sneppen, K., Nielsen, B. F., Taylor, R. J., Simonsen, L. Overdispersion in COVID-19 increases the effectiveness of limiting nonrepetitive contacts for transmission control. Proc. Nat. Acad. Sci. 118 (2021).
  * [2] Park, Y. J. et al. Contact tracing during coronavirus disease outbreak, South Korea, 2020. Emerg. Infect. Dis. 26, 2465 (2020).
  * [3] Parmet, W. E., Sinha, M. S. Covid-19—the law and limits of quarantine. N. Engl. J. Med. 382, e28 (2020).
  * [4] Hethcote, H. W. The mathematics of infectious diseases. SIAM Rev. 42, 599 (2000).
  * [5] Newman, M. E. Networks (Oxford university press, 2018).
  * [6] Gould, G. M. Anomalies and curiosities of medicine (Blacksleet River, 1966).
  * [7] Taubenberger, J. K., Morens, D. M. 1918 Influenza: the mother of all pandemics. Rev. Biomed. 17, 69 (2006).
  * [8] Chan-Yeung, M.,Xu, R.-H. SARS: epidemiology. Respirology 8, S9 (2003).
  * [9] WHO, WHO Coronavirus(COVID-19) Dashboard, https://covid19.who.int/ (2020).
  * [10] Ren, X. Pandemic and lockdown: a territorial approach to COVID-19 in China, Italy and the United States. Eurasian. Geogr. Econ. 61, 423 (2020).
  * [11] Lu, G. et al. COVID-19 in Germany and China: mitigation versus elimination strategy. Global health action 14, 1875601 (2021).
  * [12] Harvey, N. Behavioral fatigue: real phenomenon, naïve construct, or policy contrivance?. Front. Psychol. 11, 2960 (2020).
  * [13] Nicola, M. et al. The socio-economic implications of the coronavirus pandemic (COVID-19): A review. Int. J. Surg. 78, 185 (2020).
  * [14] Ferretti, L. et al. Quantifying SARS-CoV-2 transmission suggests epidemic control with digital contact tracing. Science 368 (2020).
  * [15] Flaxman, S. et al. Estimating the effects of non-pharmaceutical interventions on COVID-19 in Europe. Nature 584, 257 (2020).
  * [16] Perra, N. Non-pharmaceutical interventions during the COVID-19 pandemic: A review. Phys. Rep. (2021).
  * [17] Maier, B. F., Brockmann, D. Effective containment explains subexponential growth in recent confirmed COVID-19 cases in China. Science 368, 742 (2020).
  * [18] Chan, L. Y. H., Yuan, B., Convertino, M. COVID-19 non-pharmaceutical intervention portfolio effectiveness and risk communication predominance. Sci. Rep. 11, 1 (2021).
  * [19] Thurner, S., Klimek, P., Hanel, R. A network-based explanation of why most COVID-19 infection curves are linear. Proc. Nat. Acad. Sci. 117, 22684 (2020).
  * [20] Nimmagadda, V., Kogan, O., Khain, E. Path-dependent course of epidemic: Are two phases of quarantine better than one?. Europhys. Lett 132, 28003 (2020).
  * [21] Mukhamadiarov, R. I. et al. Social distancing and epidemic resurgence in agent-based Susceptible-Infectious-Recovered models, Sci. Rep. 11, 1 (2021).
  * [22] Choi, K., Choi, H., Kahng, B. Covid-19 epidemic under the K-quarantine model: Network approach. Preprint at https://arxiv.org/abs/2010.07157 (2020).
  * [23] Arenas, A. et al. Modeling the spatiotemporal epidemic spreading of COVID-19 and the impact of mobility and social distancing interventions. Phys. Rev. X 10, 041055 (2020).
  * [24] Schlosser, F. et al. COVID-19 lockdown induces disease-mitigating structural changes in mobility networks. Proc. Nat. Acad. Sci. 117, 32883 (2020).
  * [25] Pastor-Satorras, R., Vespignani, A. Epidemic dynamics and endemic states in complex networks. Phys. Rev. E 63, 066117 (2001).
  * [26] Newman, M. E. Spread of epidemic disease on networks. Phys. Rev. E 66, 016128 (2002).
  * [27] Balcan, D. et al. Multiscale mobility networks and the spatial spreading of infectious diseases. Proc. Nat. Acad. Sci. 106, 21484 (2009).
  * [28] Zuzek, L. A., Stanley, H. E., Braunstein, L. A. Epidemic model with isolation in multilayer networks. Sci. Rep. 5, 1 (2015).
  * [29] Fraser, C., Riley, S., Anderson, R. M., Ferguson, N. M. Factors that make an infectious disease outbreak controllable. Proc. Nat. Acad. Sci. 101, 6146 (2004).
  * [30] Peak, C. M., Childs, L. M., Grad, Y. H., Buckee, C. O. Comparing nonpharmaceutical interventions for containing emerging epidemics. Proc. Nat. Acad. Sci. 114, 4023 (2017).
  * [31] Tong, Z.-D. et al. Potential presymptomatic transmission of SARS-CoV-2, Zhejiang province, China, 2020. Emerg. Infect. Dis. 26, 1052 (2020).
  * [32] Bai, Y. et al. Presumed asymptomatic carrier transmission of COVID-19. Jama 323, 1406 (2020).
  * [33] Zhang, R. et al. Identifying airborne transmission as the dominant route for the spread of COVID-19. Proc. Nat. Acad. Sci. 117, 14857 (2020).
  * [34] Solomonoff, R., Rapoport, A. Connectivity of random nets. Bull. Math. Biol. 13, 107 (1951).
  * [35] Erdös, P., Rényi, A. On the evolution of random graphs. Publ. Math. Inst. Hung. Acad. Sci 5, 17 (1960).
  * [36] Webster, R. K. et al. How to improve adherence with quarantine: rapid review of the evidence. Public Health 182, 163 (2020).
  * [37] Fraser, C. Estimating individual and household reproduction numbers in an emerging epidemic. PloS One 2, e758 (2007).
  * [38] Cori, A., Ferguson, N. M., Fraser, C., Cauchemez, S. A new framework and software to estimate time-varying reproduction numbers during epidemics. Am. J. Epidemiol. 178, 1505 (2013).

Fig. 1 $|$ Schematic diagram for the SIR model with isolation on DLMNs. (a)
The schematic diagram for interaction on the DLMN. $W$ and $F$ layers are
composed of the same set of nodes, i.e., the nodes at the ends of each dotted
line are identical. The yellow infected node chooses random partners on $W$ at
each time step (dashed arrows). Thus, if $t_{1}\neq t_{2}$ then it interacts
with different nodes on $W$, while its interacting partner does not change on
$F$ (black solid lines). (b) The change of states under the BIP or RIP. Red,
yellow, and green boxes denote the states for $\sigma_{1}$, and white boxes
represent the states for $\sigma_{2}$.

Fig. 2 $|$ The fraction of nodes in each state.

(a)-(d) are {$\rho_{m}(t)$}’s ($m\in\left\\{S,I,R,H,X\right\\}$) under the BIP
with (a) $\theta^{*}=0.3$, $t^{*}=4$, (b) $\theta^{*}=0.3$, $t^{*}=10$, (c)
$\theta^{*}=0.9$, $t^{*}=4$, (d) $\theta^{*}=0.9$, $t^{*}=10$. (e)-(h) shows
{$\rho_{m}(t)$}’s under the RIP with (e) $\theta^{*}=0.3$, $t^{*}=4$, (f)
$\theta^{*}=0.3$, $t^{*}=10$, (g) $\theta^{*}=0.9$, $t^{*}=4$, (h)
$\theta^{*}=0.9$, $t^{*}=10$. $\beta=0.2$ and $t_{R}=6$ are used in common.
Each data is obtained from 500 independent simulations by averaging the
surviving samples at time $t$.

Fig. 3 $|$ The effective reproduction number. Plot of $R_{e}(t)$ under the (a)
BIP and (b) RIP when $t^{*}=2\sim 12$ and $\theta^{*}=0.9$ with $\beta=0.2$,
$t_{R}=6$. Plot of $R_{e}(t)$ for the (c) BIP and (d) RIP when $t^{*}=10$ and
$\theta^{*}=0,3\sim 1.0$ with $\beta=0.2$, $t_{R}=6$. The black dashed curve
denotes the case when the isolation protocol is absent. The red dotted
horizontal line depicts $R_{e}(t)=1$.

Fig. 4 $|$ Difference in the final epidemic size and the number of isolated
nodes per unit time. Plot of (a) $\Delta_{R}$ and (b) $\Delta_{XH}$ against
$\theta^{*}$ for $t^{*}=2\sim 12$ when $\beta=0.2$ and $t_{R}=6$. The black
dashed horizontal line depicts (a) $\Delta_{R}=0$ and (b)$\Delta_{XH}=0$.

## Methods

## Underlying Networks

Let $N$ be the number of nodes in a network. Each node belongs to the networks
defined on both layers, $F$ and $W$ at the same time. Thus we investigate the
epidemics on multiplex networks. To construct a fixed random network on $F$,
we randomly select two nodes among $N$ nodes and connect them if they are not
linked. This process continues until we have $L=N\left<k\right>/2$ links on
the $F$ layer. Here $\left<k\right>$ is the mean degree of the network on $F$.
The degree distribution of the obtained network on $F$, $P_{F}(k)$, is known
to be the Poisson distribution,
$P(k)=\left<k\right>^{k}e^{-\left<k\right>}/k!$ [5, 34, 35]. The links on $F$
do not change in time.

On the other hand, the topology of the network on $W$ changes in time.
Therefore, at each time $t$, a node $i$ with $\sigma_{i,1}=I$ randomly chooses
$k_{i}^{\prime}$ neighbors on $W$. $k_{i}^{\prime}$ is drawn from the Poisson
distribution,
$P_{W}(k^{\prime})=\left<k^{\prime}\right>^{k^{\prime}}e^{-\left<k^{\prime}\right>}/k^{\prime}!$.
We use the mean degrees $\left<k\right>=\left<k^{\prime}\right>=8$. The value
of the mean degree only affects the epidemic threshold[5], and does not change
the main conclusion of our study.

## Basic Isolation Protocol (BIP)

To specify the update rule for each protocol, we introduce additional
parameters $T^{I}$, $T^{X}$, and $T^{H}$ which correspond to the time of
infection, self-isolation, and hospitalization, respectively. The state of
each node at time $t$ in the DLMN is described by a two-component variable
${\bm{\sigma}}=(\sigma_{1},\sigma_{2})$.
$\sigma_{1}(\in\left\\{S,I,R\right\\})$ represents the disease state and
$\sigma_{2}(\in\left\\{U,X,H\right\\})$ denotes the isolation state.

In the BIP model, all nodes are initially in the state ${\bm{\sigma}}=(S,U)$.
Then a node $i$ is randomly selected and set to be ${\bm{\sigma}}_{i}=(I,U)$
and $T^{I}_{i}=0$. At each time step $t$, three processes are repeated for all
infected nodes whose time of infection is $T^{I}<t$: (1) infection, (2)
isolation, and (3) unisolation and recovery. Each process is composed of the
following sub-processes. Infection: (1-i) Each node $i$ with
${\bm{\sigma}}_{i}=(I,U)$ or $(I,X)$ transmits the disease to $k_{i}$
connected nodes on $F$ with the probability $\beta_{F}$ if the state of the
connected node $j$ is ${\bm{\sigma}}_{j}=(S,U)$ or $(S,X)$. (1-ii) If
$\sigma_{i,2}=U$, then it randomly chooses $k^{\prime}_{i}$ nodes on $W$ and
infects with probability $\beta_{W}$ when the randomly chosen node $j$ is in
the state ${\bm{\sigma}}_{j}=(S,U)$. $T^{I}_{j}$ for all new infected nodes
$j$ is set to be $T^{I}_{j}=t$. Isolation: (2-i) Each infected node $i$ with
$T^{I}_{i}<t$ is hospitalized with the probability $\theta_{H}$, i.e.,
$\sigma_{i,2}=H$ and $T^{H}_{i}=t$. (2-ii) Let $\Gamma_{i}$ be the set of
nodes connected to the hospitalized node $i$ on $F$. Then we set
$\sigma_{j,2}=X$ and $T^{X}_{j}=t$ for all $j\in\Gamma_{i}$ with the
probability $\theta_{X}$, if $\sigma_{j,2}=U$ at $t-1$. This corresponds to
the self-isolation. (2-iii) If the state of node $i$ was
${\bm{\sigma}}_{i}=(I,X)$ at $t-1$, then it becomes ${\bm{\sigma}}_{i}=(I,H)$
and $T^{H}_{i}=t$ for all $i(=1,2,\cdots,N)$. Unisolation and recovery: (3-i)
For all nodes $i$ with $\sigma_{i,1}=I$ becomes $\sigma_{i,1}=R$, if
$t>T^{I}_{i}+t_{R}$. Here $t_{R}$ is a constant representing a recovery time
from the infection. (3-ii) For all nodes $i$ with
$\sigma_{i,1}\in\left\\{S,R\right\\}$ and $\sigma_{i,2}=X$ are unisolated if
$t>T_{i}^{X}+t_{X}$, where $t_{X}$ is the duration time for self-isolation.
(3-iii) For all nodes $i$ with ${\bm{\sigma}}_{i}=(R,H)$, if
$t>T^{H}_{i}+t_{H}$ then the node $i$ is released from hospitalization, and
its state becomes ${\bm{\sigma}}_{i}=(R,U)$. Here $t_{H}$ denote the duration
time for hospitalization. These processes are repeated until there left no
infected node. Under the BIP only the confirmed patients and the one who has
direct contact with the confirmed patient are isolated as shown in Extended
Data Fig. 1. This is the most popular quarantine protocol adopted by many
health authorities over the globe [36].

## Reinforced Isolation Protocol (RIP)

However, household or workplace infection is significant in some infectious
diseases, for example, it is more than 15% for COVID-19 [2]. Such transmission
frequently occurs during when the patients are self-isolated or refrained from
the unnecessary social activities by themselves. Furthermore, such
transmissible paths can become a part of super-spreading events if the pre-
symptomatic or asymptomatic infection is possible. Thus controlling such local
contacts is another important factor to mitigate the transmission. For this
purpose, we introduce the RIP by adding the sub-process (2-iv) to the end of
the isolation process of the BIP: (2-iv) Let $j$ be the node whose state is
changed into $\sigma_{j,2}=X$ at $t$ and $\Gamma_{j}$ be the set of nodes
connected to $j$ on $F$. Then the state of node $n(\in\Gamma_{j})$ with
$\sigma_{n,2}=U$ also becomes $\sigma_{n,2}=X$ with the same probability
$\theta_{X}$ for all $n$. Thus, in the RIP model if a node is hospitalized,
then its first and second nearest neighbors on $F$ are isolated with the given
probability (see Extended Data Fig. 1).

## Final epidemic size

The final epidemic size under each protocol $Y$($=BIP$ or $RIP$) is defined as
$\rho^{*}_{R,Y}\equiv\rho_{R,Y}(t\rightarrow\infty)$. The difference in final
epidemic size under each protocol is defined as
$\Delta_{R}\equiv\left<\rho_{R,BIP}^{*}\right>_{s}-\left<\rho_{R,RIP}^{*}\right>_{s}$.
Here $\left<...\right>_{s}$ denotes the sample average over independent runs.
In order to obtain the average final epidemic size, we used $10,000$ samples.
The average final epidemic sizes, $\left<\rho_{R,BIP}^{*}\right>_{s}$ and
$\left<\rho_{R,RIP}^{*}\right>_{s}$ with various parameter sets,
$(\theta^{*},t^{*})$, are plotted in Extended data Fig. 6.

## Fraction of isolated nodes per unit time

To compare the number of isolated nodes we define the fractions of the
isolated nodes per unit time (with $\sigma_{2}=X$ or $H$) for protocol $Y$ as,

$\left<\rho_{XH,Y}\right>_{t}=\frac{\int_{0}^{T_{final}}(\rho_{X}(t)+\rho_{H}(t))dt}{T_{final}},$
(2)

where $T_{final}$ represents the time at which $\rho_{I}(t)$ becomes zero.
Then we define the difference between two protocols as
$\Delta_{XH}\equiv\left<\left<\rho_{XH,BIP}\right>_{t}\right>_{s}-\left<\left<\rho_{XH,RIP}\right>_{t}\right>_{s}$.
Here $\left<...\right>_{s}$ denotes the sample average over independent runs.
The fraction of isolated nodes per unit time,
$\left<\left<\rho_{XH,BIP}\right>_{t}\right>_{s}$ and
$\left<\left<\rho_{XH,RIP}\right>_{t}\right>_{s}$, are obtained from 500
independent trajectories and plotted in Extended data Fig. 7.

## Code availability

All code used in this work is available from the corresponding authors upon
reasonable request.

###### Acknowledgements.

This research was supported by Basic Science Research Program through the
National Research Foundation of Korea (NRF) funded by the Ministry of
Education (Republic of Korea) (grant number: NRF-2019R1F1A1058549).

Extended Data Fig. 1 $|$ Schematic diagram of the intervention protocol on $F$
layer. (a) On a network with size $N=8$ the node $A$ is infected while nodes
$B\sim H$ are susceptible. (b) Node $A$ is hospitalized with probability
$\theta_{H}$. (c) BIP and RIP: Node $B$ and $C$ are self-isolated with
probability $\theta_{X}$ while node $D$ remains unisolated with probability
$1-\theta_{X}$. (d) RIP: In addition to (c) node $E$ and $G$ are self-isolated
with probability $\theta_{X}$ while node $F$ remains unisolated with
probability $1-\theta_{X}$. Node $H$ are not included in the self-isolation
candidate since node $D$ is not self-isolated in step (c).

Extended Data Fig. 2 $|$ Plot of $\\{\rho_{m}(t)\\}$’s
($m\in\left\\{S,I,R,X,H\right\\}$) under the BIP. $t^{*}$ increases from
$t^{*}=2$ to $t^{*}=12$ (from left to right). $\theta^{*}$ changes from
$\theta^{*}=0.1$ to $\theta^{*}=1.0$ (from top to bottom). Each data is
obtained from 500 independent simulations by averaging the surviving samples
at time $t$.

Extended Data Fig. 3 $|$ Plot of $\\{\rho_{m}(t)\\}$’s
($m\in\left\\{S,I,R,X,H\right\\}$) under the RIP. $t^{*}$ increases from
$t^{*}=2$ to $t^{*}=12$ (from left to right). $\theta^{*}$ changes from
$\theta^{*}=0.1$ to $\theta^{*}=1.0$ (from top to bottom). When both
$\theta^{*}$ and $t^{*}$ are large, oscillatory behaviors are observed.

Extended Data Fig. 4 $|$ Effective reproduction number under the BIP. Each
plot shows the obtained $R_{e}(t)$’s for various $t^{*}(=2\sim 12)$ under the
BIP when $\theta^{*}$ is fixed. $\theta^{*}$ increases by $0.1$ from
$\theta^{*}=0.1$ to $\theta^{*}=1.0$ (from top-left to bottom-right). As shown
in the data, $t^{*}$ rarely affect the behavior of $R_{e}(t)$. However, the
increase of $\theta^{*}$ leads to a drastic drop of $R_{e}(t)$.

Extended Data Fig. 5 $|$ Effective reproduction number under the RIP. Plot of
$R_{e}(t)$’s under the RIP for $t^{*}=2\sim 12$ and $\theta^{*}=0.1\sim 1.0$.
Like $R_{e}(t)$’s under the BIP in the Extended Data Fig. 4, $t^{*}$ does
affect the behavior of $R_{e}(t)$, while $\theta^{*}$ significantly changes
the behavior of $R_{e}(t)$. In addition, we find that there are multiple peaks
in $R_{e}(t)$ for $\theta^{*}>0.3$ when we the RIP is applied. This implies
that, even though the preemptive isolation of the susceptible nodes
effectively reduces the infection rate, the instantaneous effective
reproduction number can increase again when the isolated susceptible nodes are
released to be free.

Extended Data Fig. 6 $|$ Final epidemic size. We plot the average final
epidemic size under the (a) BIP and (b) RIP for $t^{*}=2\sim 12$ and
$\theta^{*}=0.1\sim 1.0$. The increment of the isolation period ($t^{*}$) has
negligible effect on the final epidemic size while the traceability
($\theta^{*}$) significantly changes the final epidemic size.

Extended Data Fig. 7 $|$ The fraction of the isolated nodes per unit time. We
plot the average fraction of the isolated nodes per unit time under the (a)
BIP and (b) RIP for $t^{*}=2\sim 12$ and $\theta^{*}=0.1\sim 1.0$. Under the
BIP, $\left<\left<\rho_{XH,BIP}\right>_{t}\right>_{s}$ hardly change as
$\theta^{*}$ increases for a fixed $t^{*}$. However, when the RIP is applied,
$\left<\left<\rho_{XH,RIP}\right>_{t}\right>_{s}$ drastically decrease when
$\theta^{*}>0.3$. and $t^{*}>6$. This suggests that with high traceability,
the RIP can effectively mitigate the epidemic spreading while minimizing the
damage on the social and economic system.

Extended Data Fig. 8 $|$ Oscillatory behaviors under the RIP. Plot of the
fractions $\rho_{H}$, $\rho_{I}$ and $\rho_{f}$ under the RIP. $\rho_{f}(t)$
is defined as the fraction of the susceptible nodes which are just released
from self-isolation at time $t$. The epidemic parameters are $\beta=0.2$ and
$t_{R}=6$, and the intervention parameters are $\theta^{*}=0.9$ and
$t^{*}=12$. The period of oscillation $\tau$ for each fraction with the given
parameters is estimated as $\tau\approx 21$. The peak position of each curve
indicates that the increase of the infected causes an increase of the
hospitalized individuals. Due to the hospitalization and self-isolation, the
number of the infected rapidly decreases. However, after $t^{*}$ the isolated
nodes are set to be free which increases the number of unisolated susceptible
nodes. Thus it increases the number of infected individuals again. This
pattern is repeated with decreasing amplitude due to the depletion of the
susceptible nodes until there is no more infected node left. This oscillatory
behavior is observed only for the case of large $\theta^{*}$ and $t^{*}$ in
the RIP.