# Multiple peaks in
network SIR epidemic models

Martina Alutto  Leonardo Cianfanelli  Giacomo Como  Member, IEEE  Fabio
Fagnani Some of the results appeared in a preliminary form in [9992408].The
authors are with the Department of Mathematical Sciences “G.L. Lagrange,”
Politecnico di Torino, 10129 Torino, Italy. E-mail:
<EMAIL_ADDRESS>G. Como is also with the Department of Automatic Control, Lund University,
22100 Lund, Sweden.This work was partially supported by a MIUR grant
“Dipartimenti di Eccellenza 2018–2022” [CUP: E11G18000350001], a MIUR Research
Project PRIN 2017 “Advanced Network Control of Future Smart Grids”
(http://vectors.dieti.unina.it) and by the Compagnia di San Paolo.

###### Abstract

We study a susceptible-infected-recovered (SIR) epidemic model on a network of
interacting subpopulations and analyze the dynamical behavior of the fraction
of infected agents in each node of the network. In contrast to the classical
scalar epidemic SIR model, in which the fraction of infected is known to have
a unimodal behavior (decreasing over time or initially increasing until
reaching a peak and then decreasing), we show the possible occurrence of new
multimodal behaviors in the network SIR model. We focus on the special case of
rank-1 network matrices, which model subpopulations of homogeneously mixing
agents with different interaction levels. We provide an upper bound on the
number of changes of monotonicity of the fraction of infected at the single
node level and provide sufficient conditions under which such multimodal
behavior occurs. We then conduct a numerical analysis, revealing that, in case
of more general network matrices, the dynamics may exhibit complex behaviors
with multiple peaks of infection in each node.

Epidemic models, Susceptible-Infected-Recovered model.

##  1 Introduction

The pandemic emergence of last years has generated a renewed huge interest on
epidemic compartment models that have proven to be effective tools for the
forecasting of virus spreading but also in assisting in the design of policy
containment rules such as distanciation and lockdown.

The simplest and most popular among these models is the SIR epidemic model
introduced almost one century ago [Kermack.McKendrick:1927] and since then
well studied in the literature [Hethcote2000TheMO, MKendrickApplicationsOM,
Diekmann2000MathematicalEO]. According to it, a population is split into three
categories: the susceptible agents, who have not yet been infected and can
catch the disease, the infected agents, who are currently carrying the
pathogen and allow the transmission of the disease, and the recovered agents,
who have healed from the infection and are for ever immune. The model assumes
that the rate of new infections is proportional to the product between the
number of susceptible and infected agents due to pairwise interactions,
implicitly assuming the homogeneous mixing of the population. The crucial
index in the SIR model is the so called _reproduction number_ $R(t)$, a time-
dependent scalar quantity which describes the average new infections that an
infected individual is producing. If $R(t)<1$, then the fraction of infected
agents is decreasing at time $t$. On the other hand, if $R(t)>1$, then the
fraction of infected agents increases. As $R(t)$ is monotonically decreasing
in $t$ and eventually becomes less than $1$, an epidemic wave, as modeled by
SIR, has necessarily a _unimodal_ behavior. Precisely, if $R(0)<1$, no spread
will occur, the number of infected individuals will be decreasing and it will
approach $0$ as time gets large. If instead $R(0)>1$, the curve of infected
will be increasing up to a maximum value (the _peak_) and then will start to
decrease. The unimodal behavior of the SIR dynamics has been shown to hold
also for more general interaction mechanisms [Alutto2021OnSE] and is at the
basis of several control strategies, including some recently proposed in the
context of the COVID-19 pandemic, see, e.g., [Cianfanelli.ea:2021] and
[Miclo.ea:2022].

The classical SIR model relies on a number of homogeneity assumptions on the
population regarding their mixing, their aptitude to contract the infection as
well the time needed to recover that can hardly be met in realistic scenarios.
This has motivated the introduction of networked versions of the SIR model
[hethcote1978immunization, Fall:2007, Nowzari.ea:2016, Mei.ea:2017,
Ogura2016StabilityOS]. In such models (that will be referenced as _network SIR
models_), each node of the network represents a subpopulation of
indistinguishable agents and may describe, depending on the applicative
contexts, geographical areas or rather population categories like age, life
style, etc… Interactions between agents of different nodes is coded into a
matrix $A$: given subpopulations $i$ and $j$, $A_{ij}$ is the rate of new
infections in node $i$ due to the presence of infected agents in node $j$ and
may incorporate the peculiar susceptibility of individuals in $i$, the rate of
interactions among members of the two subpopulations, and, possibly, the
effect of targeted containment policies. We shall refer to $A$ as the
_interaction matrix_ of the SIR network model.

Several recent papers [acemoglu2021optimal, NBERw27063] use calibrated network
SIR models to examine the impact of age-targeted mitigation policies for the
COVID-19 pandemic showing how such policies (even with just two age groups)
can outperform uniform intervention policies in terms of both mortality rates
and economic productivity.

While most of the studies on the network SIR model are empiric, there are two
notable exceptions. In [Mei.ea:2017] the authors discover a novel network
reproduction number that is a decreasing function of time converging to $0$
and that plays an equivalent role as the one in the scalar SIR model. When
this reproduction number is below one, then a certain aggregated infection
index (a linear combination of the number of infected agents in the various
subpopulations) decreases. When instead this number is above one, this
aggregated infection index will first increase and, once the reproduction
number becomes smaller than 1, it will start decreasing to $0$. However, this
aggregated infection index is defined through weights that depend on the
initial condition and are possibly time varying and this limits its possible
applications. In [ellison2020implications] the authors analyze a network SIR
model with symmetric rank-$1$ network matrices, which results from assuming
that agents have different interaction levels but there is no homophily in the
society (the analysis is then generalized to include homophily). It is shown
that the heterogeneity affects the final size of the outbreak, namely it is
possible to reach herd immunity with an aggregate fraction of infected smaller
than in the scalar SIR. At the best of our knowledge, there is no analysis in
the literature on the behavior of the curve of infected at single nodes which
is relevant in understanding the effectiveness of targeted interventions.

This paper gives a novel theoretical contribution to understanding the network
SIR model with $n$ nodes in the special case when the interaction matrix has
rank $1$. This is a relevant case previously analyzed in [acemoglu2021optimal,
ellison2020implications]. Our results are threefold. First, we provide $n$
invariants of motion and a new aggregated infection index with fixed weights
only depending on the matrix $A$, which always exhibits a unimodal behavior as
function of time, analogously to the curve of infected in the scalar case.
Second, we carry on a node level analysis proving that the curve of infected
at every node can undertake at most two changes of monotonicity before the
reproduction number gets below $1$, and since then it is monotonically
decreasing. Third, we exhibit a class of network SIR model with just two nodes
where the curve of infected at one of the two nodes effectively does not
present an unimodal behavior, namely an initial decreasing behavior with a
subsequent local minimum point, leading to an increasing behavior to a peak
and then decreasing to zero.

We are aware that the phenomenon of multiple waves of infection cannot be
totally explained by the heterogeneity introduced by the network and is also
largely determined by the adaptive behavior and endogenous response of agents
to the epidemic, as well as by the phenomenon of waving immunity
[Ehrhardt2019]. For example, some papers have studied models that take into
account how agents adapt their behavior, resulting in a modification of the
parameters of the model at macroscopic level [CAPASSO197843, Alutto2021OnSE],
or in relation to loss of immunity over time [Jin2007, Li2014].

We also acknowledge the many possible extensions of the network SIR model to
more than three compartments to keep track for instance of the many forms of
infection and possibly vaccination [Li1995, Giordano_2020, Birge.ea:2020,
parino2021modelling]. The possibility to extend our results to more complex
models is left as future research.

The rest of the paper is organized as follows. In Section 2 we describe the
network SIR model and summarize the results known in the literature. In
Section 3 we state our main result, which characterizes all the possible
behaviors that the dynamics may exhibit at single node levels when the
interaction matrix has rank 1. In Section 4, we provide sufficient conditions
for the existence of multimodal behaviors. In Section 5 we illustrate
numerical simulations on more generale networks. Finally, in Section 6 we
summarize our work and discuss future research lines.

###  1.1 Notation

Here we briefly gather some notational conventions adopted throughout the
paper. We denote by $\mathbb{R}$ and $\mathbb{R}_{+}$ the sets of real and
nonnegative real numbers, respectively, while $\mathbb{R}_{+}^{n\times n}$
indicates the set of real matrices with dimension $n\times n$ and nonnegative
entries. The all-1 vector and the all-0 vector are denoted by $\bm{1}$ and
$\bm{0}$ respectively, where the size of them may be deduced from the context.
Given a vector $x$, we let $x^{T}$ denote its transpose, and $[x]$ indicate
the diagonal matrix whose entries coincide with vector $x$. For an irreducible
nonnegative matrix $A$, we let $\lambda_{max}(A)$ and $v_{max}(A)$ denote
respectively the dominant eigenvalue of $A$ and the corresponding left
eigenvector, which is unique and has positive entries due to the Perron-
Frobenius theorem. Inequalities between two vectors $x$ and $y$ in
$\mathbb{R}^{n}$ are meant to hold true entry-wise, i.e., e.g., $x\leq y$
means that $x_{i}\leq y_{i}$ for every $i$, whereas $x<y$ means that
$x_{i}<y_{i}$ for every $i$, and whereas $x\lneq y$ means that $x_{i}\leq
y_{i}$ for every $i$ and $x_{j}<y_{j}$ for some $j$ .

##  2 Network SIR epidemic model

In this section, we introduce the network SIR epidemic model and gather some
known results that will prove useful in the sequel.

We shall model networks as finite weighted directed graphs
$\mathcal{G}=(\mathcal{V},\mathcal{E},A)$, where
$\mathcal{V}=\\{1,2,\ldots,n\\}$ is the set of nodes,
$\mathcal{E}\subseteq\mathcal{V}\times\mathcal{V}$ is the set of directed
links, and $A$ in $\mathbb{R}_{+}^{n\times n}$ is a nonnegative weight matrix,
to be referred as the interaction matrix, with the property that $A_{ij}>0$ if
and only if there exists a link $(i,j)$ in $\mathcal{E}$ directed from node
$i$ to node $j$. A network is referred to as strongly connected if its
interaction matrix $A$ is irreducible.

In a network SIR epidemic model, a set of interacting subpopulations
$\mathcal{V}$ are identified with the nodes of a network $\mathcal{G}$. For
every subpopulation $i$ in $\mathcal{V}$ the time-varying variables $x_{i}$,
$y_{i}$, and $z_{i}$ represent the fractions of susceptible, infected, and
recovered individuals, respectively, so that sum

$x_{i}+y_{i}+z_{i}=1$

remains constant in time. The entries $A_{ij}$ of the interaction matrix
represent the product between the infection rate and the contact frequency
between agents of subpopulation $i$ and agents of subpopulation $j$. Finally,
a positive scalar parameter $\gamma$ models the recovery rate, which is
assumed to be homogeneous across the network.

The network SIR epidemic model with interaction matrix $A$ and recovery rate
$\gamma$ is then the autonomous system of ordinary differential equations

$\begin{cases}\dot{x}_{i}=-x_{i}\sum_{j}A_{ij}y_{j}\\\
\dot{y}_{i}=x_{i}\sum_{j}A_{ij}y_{j}-\gamma y_{i}\\\ \dot{z}_{i}=\gamma
y_{i}\,,\end{cases}$ (1)

for every $i=1,\ldots,n$. Notice that the third equation is redundant since
$z_{i}(t)=1-x_{i}(t)-y_{i}(t)$ for every subpopulation $i$ in $\mathcal{V}$
and time $t\geq 0$. The network SIR epidemic model (1) can then be more
compactly rewritten in the following vectorial form

$\dot{x}=-[x]Ay\,,\qquad\dot{y}=[x]Ay-\gamma y\,,$ (2)

where $x$ and $y$ in $\mathbb{R}^{n}_{+}$ denote the vectors of susceptible
and infected individuals, respectively, in all the different subpopulations.

The following result gathers some basic properties of the network SIR model.

######  Proposition 1.

Consider the network SIR epidemic model (2) with irreducible interaction
matrix $A$ and recovery rate $\gamma>0$. Then,

  1. (i)

the set $\mathcal{S}=\\{(x,y)\in[0,1]^{2n}:x+y\leq\bm{1}\\}$ is invariant;

  2. (ii)

the set of equilibrium points in $\mathcal{S}$ is

$\mathcal{S}^{*}=\\{(x^{*},\bm{0}):\,x^{*}\in[0,1]^{n}\\}\,;$

  3. (iii)

an equilibrium point $(x^{*},\bm{0})$ in $\mathcal{S}^{*}$ is stable if and
only if

$\lambda_{\max}([x^{*}]A)<\gamma\,.$

Moreover, for every initial condition $(x(0),y(0))$ in $\mathcal{S}$:

  1. (iv)

for every $i=1,\ldots,n$, $x_{i}(t)$ is non-increasing for $t\geq 0$, and
$x_{i}(0)>0$ if and only if $x_{i}(t)>0$ for every $t\geq 0$;

  2. (v)

if $y(0)\gneq\bm{0}$, then $y(t)>\bm{0}$ for every $t>0$;

  3. (vi)

there exists $\bm{0}\leq x^{*}\leq x(0)$ such that

$\lim_{t\to+\infty}x(t)=x^{*}\,,\qquad\lim_{t\to+\infty}y(t)=\bm{0}\,.$

###### Proof.

See [Nowzari.ea:2016] and [Mei.ea:2017]. ∎

In the special case when $n=1$, so that the interaction matrix reduces to a
positive scalar value $A=\beta>0$ representing the contagion rate, the network
SIR epidemic model (2) reduces to the classical scalar SIR epidemic model

$\dot{x}=-\beta xy\,,\qquad\dot{y}=(\beta x-\gamma)y\,.$ (3)

For the scalar SIR epidemic model (3), a more refined analysis is available.
In particular, the following fundamental result is known to hold true.

######  Proposition 2.

For the scalar SIR epidemic model (3), with contagion rate $\beta>0$, recovery
rate $\gamma>0$, and initial condition $(x(0),y(0))$ such that $0<x(0)\leq
1-y(0)\leq 1$,

  1. (i)

the quantity $\beta(x+y)-\gamma\log x$ is an invariant of motion.

Moreover, if $y(0)>0$, then:

  1. (ii)

if $\beta x(0)\leq\gamma$, then $y(t)$ is strictly decreasing for $t\geq 0$;

  2. (iii)

if $\beta x(0)>\gamma$, then there exists a peak time $\hat{t}>0$ such that
$y(t)$ is strictly increasing for $t$ in $[0,\hat{t}]$ and strictly decreasing
for $t$ in $[\hat{t},+\infty)$;

  3. (iv)

the limit value $x^{*}=\lim\limits_{t\to+\infty}x(t)$ is the unique solution
of the equation

$\beta x^{*}-\gamma\log x^{*}=\beta(x(0)+y(0))-\gamma\log x(0)\,,$ (4)

in the interval $[0,\gamma/\beta]$.

###### Proof.

See [Brauer2019MathematicalMI, Chapter 2.4]. ∎

We now provide a simple example of a network SIR epidemic model with just two
nodes where the curve of infected agents at single node level has multiple
peaks, for a range of initial conditions.

######  Example 1.

Figure 1: Numerical simulation of the network SIR epidemic model with $n=2$
nodes with interaction matrix $A=\bm{1}\bm{1}^{\prime}$, recovery rate
$\gamma=1$, and initial condition $y_{1}(0)=1-x_{1}(0)=0.2$ and
$y_{2}(0)=1-x_{2}(0)=0$ satisfying (5)-(6).

Consider the network SIR epidemic model (1) with $n=2$ subpopulations,
interaction matrix $A=\bm{1}\bm{1}^{\prime}$, unitary recovery rate
$\gamma=1$, and initial condition

$y_{1}(0)=1-x_{1}(0)=\varepsilon\,,\quad y_{2}(0)=1-x_{2}(0)=0\,,$ (5)

for some $\varepsilon>0$ is such that

$\frac{1-\varepsilon}{2-\varepsilon}(1-\log(2-\varepsilon))>\varepsilon\,.$
(6)

Notice that a range of such values of $\varepsilon$ always exist since the
function
$g(\varepsilon)=\frac{1-\varepsilon}{2-\varepsilon}(1-\log(2-\varepsilon))-\varepsilon$
is continuous in the interval $[0,1]$ and such that $g(0)=1/2(1-\log 2)>0$.

Observe that with these initial conditions we have

$\dot{y}_{1}(0)=x_{1}(0)(y_{1}(0)+y_{2}(0))-y_{1}(0)=-\varepsilon^{2}<0\,,$
(7)

which implies that $y_{1}(t)$ is strictly decreasing for sufficiently small
$t>0$. We will now show that $y_{1}(t)$ cannot remain decreasing for all
values of $t\geq 0$, but will necessary become increasing in a certain time
range, before eventually starting to decrease again and vanish as $t$ grows
large.

Towards this goal, first observe that the aggregate variables
$\overline{x}=x_{1}+x_{2}$ and $\overline{y}=y_{1}+y_{2}$ satisfy an
autonomous scalar SIR epidemic model

$\dot{\overline{x}}=-\overline{x}\,\overline{y}\,,\qquad\dot{\overline{y}}=(\overline{x}-1)\overline{y}\,.$
(8)

Then, since $\dot{\overline{y}}(0)=(\overline{x}(0)-1)\overline{y}(0)>0\,,$
Proposition 2 (iii) implies that there exists a peak time $\hat{t}>0$ at which
$\dot{\bar{y}}(\hat{t}\,)=0$, i.e., $\overline{x}(\hat{t}\,)=1$. This,
Proposition 2(i) and (5) imply that

$\overline{y}(\hat{t}\,)=\overline{x}(0)+\overline{y}(0)-\overline{x}(\hat{t}\,)+\log\frac{\overline{x}(\hat{t}\,)}{\overline{x}(0)}\\\
=1-\log(2-\varepsilon)\,.$ (9)

Since $\dot{x}_{2}=-x_{2}\overline{y}$ and $x_{2}(0)=1$, we have that

$x_{2}(\hat{t}\,)=\exp\Big{(}-\int_{0}^{\hat{t}}\overline{y}(t)\mathrm{d}t\Big{)}=\frac{\overline{x}(\hat{t}\,)}{\overline{x}(0)}=\frac{1}{2-\varepsilon}\,,$
(10)

where the second equality follows from integrating the first equation in (8)
and the last one follows from (5). It then follows from (9) and (10) that

$\dot{y}_{2}(\hat{t}\,)=x_{2}(\hat{t}\,)\overline{y}(\hat{t}\,)-y_{2}(\hat{t}\,)=\frac{1-\log(2-\varepsilon)}{2-\varepsilon}-y_{2}(\hat{t}\,)\,.$
(11)

Now, assume by contradiction that $\dot{y}_{1}(t)\leq 0$ for all $t\geq 0$. In
particular, this would imply that $y_{1}(\hat{t}\,)\leq y_{1}(0)=\varepsilon$,
so that

$y_{2}(\hat{t}\,)=\overline{y}(\hat{t}\,)-y_{1}(\hat{t}\,)\geq
1-\log(2-\varepsilon)-\varepsilon\,.$

by (9). Recalling that $\dot{\overline{y}}(\hat{t}\,)=0$, substituting the
above in the righthand side of (11), and using (6) we would then get

$\dot{y}_{1}(\hat{t}\,)=\dot{\overline{y}}(\hat{t}\,)-\dot{y}_{2}(\hat{t}\,)\geq\frac{1-\varepsilon}{2-\varepsilon}\left(1-\log(2-\varepsilon)\right)-\varepsilon>0\,.$

It then follows that the must exist some values of $t\geq 0$ such that
$\dot{y}_{1}(t)>0$. Together with (7) and the fact
$\lim_{t\to+\infty}y_{1}(t)=0$ by Proposition 1(vi), this implies that
$y_{1}(t)$ has a multimodal behavior.

In fact, the results in Section 3 will imply that such behavior is necessarily
as illustrated in Figure 1: $y_{1}(t)$ is strictly decreasing in an interval
$[0,\check{t}_{1}]$, until reaching a positive local minimum point
$\check{t}_{1}>0$, it is then strictly increasing in an interval
$[\check{t}_{1},\hat{t}_{1}]$ until reaching a second peak at some time
$\hat{t}_{1}>\check{t}_{1}$, and is eventually strictly decreasing for
$t\geq\hat{t}_{1}$.

Somewhat surprisingly, the network SIR epidemic model considered in this
example can actually be interpreted as a scalar SIR epidemic model where the
population of agents has simply been split into two equally sized
subpopulations, which have distinct initial conditions. Specifically, the
first population has all the initially infected people while the second
population is totally constituted by susceptible individuals. The parameters
chosen are such that if the two subpopulations were isolated, the first
subpopulation would undertake an exponential decrease to a disease-free state.
However, because of the presence of the second subpopulation, the infection
can further spread and eventually hit back the first subpopulation that thus
suffers a second wave of infection with a second peak (the first one being at
time $0$).

In the next section we study the network SIR epidemic model with rank-$1$
interaction matrices and provide results on the changes of monotonicity of the
curve of the infected fraction of individuals in each subpopulation. In
particular, we show that, for arbitrary rank-$1$ interaction matrices, the
number of such changes of monotonicity can never exceed two.

##  3 The network SIR model with rank-$1$ interaction matrices

In this section, we study the network SIR model in the special case when the
interaction matrix $A$ has rank $1$, as per the following equivalent
assumption.

######  Assumption 1.

The interaction matrix $A$ satisfies

$A=ab^{T}\,,$ (12)

for two vectors $a>\bm{0}$ and $b>\bm{0}$ in $\mathbb{R}^{n}$.

######  Remark 1.

Notice that this case encompasses the one studied in [ellison2020implications]
where authors impose the extra condition that $A$ is symmetric.

Let us define the weighted sums of susceptible fraction of individuals

$\bar{x}=\sum_{j=1}^{n}b_{j}x_{j}\,,$ (13)

and, respectively, of infected fraction of individuals

$\bar{y}=\sum_{j=1}^{n}b_{j}y_{j}\,,$ (14)

across the network. Notice that, for rank-1 interaction matrices $A=ab^{T}$,
the network SIR epidemic model’s equations (2) can be rewritten as

$\dot{x}_{i}=-a_{i}x_{i}\bar{y}\,,$ (15a)
$\dot{y}_{i}=a_{i}x_{i}\bar{y}-\gamma y_{i}\,,$ (15b)

for every $i=1,\ldots,n$. Moreover, let

$\tilde{x}=\sum_{j=1}^{n}a_{j}b_{j}x_{j}\,,$ (16)

and

$w_{i}=\tilde{x}-\gamma-a_{i}\overline{y}\,,$ (17)

for every $i=1,\ldots,n$.

###  3.1 Invariants of motion and unimodality in the aggregate

We have the following technical result.

######  Lemma 1.

Consider the rank-$1$ network SIR epidemic model (15a)–(15b). Then,

$\dot{\overline{x}}=-\overline{y}\tilde{x}\,,$ (18)

and

$\dot{\overline{y}}=\overline{y}\left(\tilde{x}-\gamma\right)\,.$ (19)

Moreover,

$\ddot{y}_{i}=a_{i}x_{i}\overline{y}w_{i}-\gamma\dot{y}_{i}\,.$ (20)

for every $i=1,\ldots,n$.

###### Proof.

See Appendix 7. ∎

Our next result generalizes Proposition 2(i) and determines $n$ invariants of
motions for the network SIR epidemic model with rank-$1$ interaction matrix.

######  Proposition 3.

Consider the rank-$1$ network SIR epidemic model (15a)–(15b). Then, for every
$i=1,\ldots,n$, the quantity

$h_{i}(x,y)=a_{i}(\overline{x}+\overline{y})-\gamma\log x_{i}\,,$

is an invariant of motion.

###### Proof.

It follows from (18), (19), and (15a) that

$\begin{array}[]{rcl}\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}h_{i}(x(t),y(t))&=&\displaystyle
a_{i}(\dot{\overline{x}}(t)+\dot{\overline{y}}(t))-\gamma\frac{\dot{x}_{i}(t)}{x_{i}(t)}\\\\[7.0pt]
&=&\displaystyle-a_{i}\gamma\overline{y}(t)+\gamma
a_{i}\overline{y}(t)\\\\[7.0pt] &=&0\,,\end{array}$

thus implying that $h_{i}(x(t),y(t))$ remains constant along the solutions of
the network SIR epidemic model (15a)–(15b). ∎

We can now prove the following result, establishing that, for the network SIR
model on rank-1 interaction networks, the weighted sum of the infected
fractions of individuals $\overline{y}$ defined in (14) has a unimodal
behavior.

######  Proposition 4.

Consider the rank-$1$ network SIR epidemic model (15a)–(15b). Then, for every
initial condition $(x(0),y(0))$ such that

$0\lneq x(0)\leq\bm{1}-y(0)\lneq\bm{1}\,,$ (21)

we have that:

  1. (i)

$\tilde{x}(t)$ is strictly decreasing for $t\geq 0$ and

$\lim_{t\to+\infty}\tilde{x}(t)<\gamma\,;$ (22)

  2. (ii)

if

$\tilde{x}(0)\leq\gamma\,,$ (23)

then $\overline{y}(t)$ is strictly decreasing for $t\geq 0$;

  3. (iii)

if

$\tilde{x}(0)>\gamma\,,$ (24)

then there exists a peak time $\hat{t}>0$ such that $\overline{y}(t)$ is
strictly increasing on $[0,\hat{t}]$ and strictly decreasing on
$[\hat{t},+\infty)$.

###### Proof.

(i) By taking the derivative of both sides of equation (16) and substituting
equation (15a), we get

$\dot{\tilde{x}}=\sum_{j=1}^{n}a_{j}b_{j}\dot{x}_{j}=-\overline{y}\sum_{j=1}^{n}a_{j}^{2}b_{j}x_{j}\,.$
(25)

Now, the rightmost inequality in (21) and Proposition 1(v) imply that $y(t)>0$
for every $t>0$, whereas the leftmost inequality in (21) and Proposition 1(iv)
imply that there exists some $i$ in $\\{1,\ldots,n\\}$ such that $x_{i}(t)>0$
for every $t\geq 0$. It then follows from (25) and the assumption that
$a>\bm{0}$ and $b>\bm{0}$ that

$\dot{\tilde{x}}(t)=-\overline{y}(t)\sum_{j=1}^{n}a_{j}^{2}b_{j}x_{j}(t)\leq-
a_{i}^{2}b_{i}x_{i}(t)\overline{y}(t)<0\,,$

for every $t>0$, which implies that $t\mapsto\tilde{x}(t)$ is strictly
decreasing for $t\geq 0$. Now, let

$\tilde{x}(\infty)=\lim_{t\to+\infty}\tilde{x}(t)\,.$

Clearly, if $\tilde{x}(0)\leq\gamma$, then
$\tilde{x}(\infty)<\tilde{x}(0)\leq\gamma$, so that inequality (22) is
satisfied. On the other hand, if $\tilde{x}(0)>\gamma$ and
$\tilde{x}(\infty)\geq\gamma$, then $\tilde{x}(t)>\gamma$ for every $t\geq 0$,
so that, by equation (19),

$\dot{\overline{y}}(t)=\overline{y}(t)\left(\tilde{x}(t)-\gamma\right)\geq
0\,,\qquad\forall t\geq 0\,.$

The above would then imply that $t\mapsto\overline{y}(t)$ is nondecreasing for
$t\geq 0$, so that

$\lim_{t\to+\infty}\overline{y}(t)\geq\overline{y}(0)>0\,,$

thus contradicting Proposition 1(vi). Therefore, also when
$\tilde{x}(0)>\gamma$ we must have $\tilde{x}(\infty)<\gamma$, thus completing
the proof of point (i) of the claim.

(ii) If $\tilde{x}(0)\leq\gamma$, by point (i) we have that
$\tilde{x}(t)<\gamma$ for every $t>0$. Hence, equation (19) implies that

$\dot{\overline{y}}(t)=\overline{y}(t)\left(\tilde{x}(t)-\gamma\right)<0\,\qquad\forall
t>0\,,$

thus showing that $\overline{y}(t)$ is strictly decreasing for $t\geq 0$.

(iii) If $\tilde{x}(0)>\gamma$, by point (i), $\tilde{x}(t)$ is strictly
decreasing for $t\geq 0$ and

$\lim_{t\to+\infty}\tilde{x}(t)<\gamma\,.$

Then, there necessarily exits a time $\hat{t}>0$ such that
$\tilde{x}(t)>\gamma$ for every $0\leq t<\hat{t}$,
$\tilde{x}(\hat{t})=\gamma$, and $\tilde{x}(t)<\gamma$ for every $t>\hat{t}$.
It then follows from equation (19) that $t\mapsto\overline{y}(t)$ is strictly
increasing for $t$ in $[0,\hat{t}]$ and strictly decreasing for $t$ in
$[\hat{t},+\infty)$, thus proving the claim. ∎

######  Remark 2.

The result in [Mei.ea:2017] suggests a sort of unimodal behavior of the
infection curve analogous to the scalar case. In particular, if
$\lambda_{\max}([x]A)(\tau)<\gamma$ for some $\tau\geq 0$, then the aggregate
curve of infected $v_{max}(\tau)^{T}y(t)$ is monotonically decreasing to 0.
Instead, if $\lambda_{\max}([x]A)(0)<\gamma$, then for small times
$v_{max}(0)^{T}y(t)$ growa exponentially fast. However, notice how the
aggregate index $v_{max}(\tau)^{T}y(t)$ explicitly depends on $x(\tau)$ and is
not clear that $v_{max}(0)^{T}y(t)$ is indeed unimodal. In our case study with
rank-1 interaction matrices $A=ab^{T}$, the left eigenvector of the $[x]A$ is
precisely $b$, so the aggregate is constant in this case and does not depend
on the time instant.

###  3.2 Dynamic behavior of the single populations

Let

$\hat{t}=\inf\\{t\geq 0:\,\tilde{x}(t)\leq\gamma\\}\,,$ (26)

and observe that Proposition 4(i) implies that $\hat{t}<+\infty$. Also, for
every $i=1,\ldots,n$, let

$\overline{t}_{i}=\inf\\{t\\!\geq\\!0:\,w_{i}(t)\\!\leq\\!0\\}=\inf\\{t\\!\geq\\!0:\,\tilde{x}(t)\\!\leq\\!\gamma+a_{i}\overline{y}(t)\\}\,,$
(27)

and notice that

$\overline{t}_{i}\leq\hat{t}\,,$ (28)

and $\overline{t}_{i}\leq\overline{t}_{j}$ if and only if $a_{j}\leq a_{i}$.
Hence, it is possible to order these time instants from the entries of vector
$a$.

We now present the following technical results that will prove useful in
deriving our main result.

######  Lemma 2.

Consider the rank-$1$ network SIR epidemic model (15a)–(15b) and every initial
condition $(x(0),y(0))$ such that $y(0)\gneq 0$. Then,

  1. (i)

for every $t\geq 0$,

$\dot{w}_{i}(t)<-a_{i}\overline{y}w_{i}(t)\,;$

  2. (ii)

$w_{i}(t)$ is strictly decreasing for $0\leq t\leq\overline{t}_{i}$;

  3. (iii)

for every $t>\overline{t}_{i}$,

$w_{i}(t)<0\,;$

  4. (iv)

if $\dot{y}_{i}(t)=0$ for some $t\geq\overline{t}_{i}$, then $t$ cannot be a
local minimum point of $y_{i}(t)$.

###### Proof.

See Appendix 8. ∎

We can now state our main result, characterizing the dynamic behavior of the
fraction of infected individuals in the single populations of the network SIR
epidemic model with rank-$1$ interaction matrix.

######  Theorem 1.

Consider the rank-$1$ network SIR epidemic model (15a)–(15b). Let
$i\in\\{1,\ldots,n\\}$ be such that $y_{i}(0)>0$. Then,

  1. (i)

$y_{i}(t)$ admits at most one local minimum time $\check{t}_{i}\geq 0$.

Moreover, if such local minimum time $\check{t}_{i}$ exists,

  1. (ii)

it satisfies

$0\leq\check{t}_{i}\leq\overline{t}_{i}\,,$ (29)

with $\check{t}_{i}=\overline{t}_{i}=0$ if and only if $w_{i}(0)\leq 0$ and
$\dot{y}_{i}(0)>0$;

  2. (iii)

it cannot occur after any stationary local maximum point of $y_{i}(t)$.

###### Proof.

If $w_{i}(0)\leq 0$ holds true, then $\overline{t}_{i}=0$ and Lemma 2(iv)
implies that no stationary point $t\geq 0$ of $y_{i}(t)$ can be a local
minimum point. It follows that the only local minimum point of $y_{i}(t)$ can
possibly be $\check{t}_{i}=0$ (which is the case if and only if
$\dot{y}_{i}(0)>0$). On the other hand, if $w_{i}(0)>0$, then the interior
extremum theorem and Lemma 2(iv) imply that there cannot be any minimum points
of $y_{i}(t)$ in the interval $[\overline{t}_{i},+\infty)$. This proves point
(ii).

We are then left with studying local minimum points of $y_{i}(t)$ in the
interval $[0,\overline{t}_{i})$. Let $s\geq 0$ be a stationary local maximum
point of $y_{i}(t)$, and let $u$ in $(s,\overline{t}_{i})$ be a (necessarily
stationary) local minimum point of $y_{i}(t)$. Then, we have that

$\dot{y}_{i}(s)=\dot{y}_{i}(u)=0\,,$ (30)

and

$\ddot{y}_{i}(s)\leq 0\,,\qquad\ddot{y}_{i}(u)\geq 0\,.$ (31)

Since $y_{i}(t)>0$ for all $t\geq 0$, we have that

$\begin{array}[]{rcl}0&\geq&\ddot{y}_{i}(s)/y_{i}(s)\\\\[5.0pt]
&=&a_{i}x_{i}(s)\overline{y}(s)w_{i}(s)/y_{i}(s)-\gamma\dot{y}_{i}(s)/y_{i}(s)\\\\[5.0pt]
&=&\gamma w_{i}(s)\\\\[5.0pt] &>&\gamma w_{i}(u)\\\\[5.0pt]
&=&a_{i}x_{i}(u)\overline{y}(u)w_{i}(u)/y_{i}(u)-\gamma\dot{y}_{i}(u)/y_{i}(u)\\\\[5.0pt]
&=&\ddot{y}_{i}(u)/y_{i}(u)\\\\[5.0pt] &\geq&0\,,\end{array}$ (32)

where the first and the last inequalities above follow from (31), the first
and the last identities follow from (20), the other two identities from (30)
and the fact that, by (15b), $a_{i}x_{i}(t)\bar{y}(t)=\gamma y_{i}(t)$ when
$\dot{y}_{i}(t)=0$, and the strict inequality in the middle holds true because
of Lemma 2(ii). As (32) is a contradiction, this shows that a local minimum
point $u<\hat{t}_{i}$ of $y_{i}(t)$ cannot follow any stationary local maximum
point $s\geq 0$ of $y_{i}(t)$, thus proving point (iii).

Finally, to prove point (i), assume by contradiction that there exist two
distinct local minimum points $r<u$ of $y_{i}(t)$ in the interval
$[0,\hat{t})$. Then, there would necessarily exist a local maximum point of
$y_{i}(t)$ in the interval $(r,u)$. But, since $s>r\geq 0$, such local maximum
point would also be stationary, thus violating point (iii). Therefore, there
cannot exist two distinct local minimum points $r<u$ of $y_{i}(t)$ in the
interval $[0,\hat{t}_{i})$, thus completing the proof of point (i). ∎

######  Remark 3.

Notice that (29) and (28) imply that, if it exists, the local minimum point of
$y_{i}$ cannot occur after the peak of $\bar{y}$.

As a consequence of Theorem 1 we get the following result classifying the
possible behaviors of the fraction of infected individuals in the single
populations of the network SIR epidemic model with rank-$1$ interaction
matrix. This classification is based on the study of the sign of two
quantities:

$\dot{y}_{i}(0)=a_{i}x_{i}(0)\overline{y}(0)-\gamma
y_{i}(0),\,w_{i}(0)=\tilde{x}(0)-\gamma-a_{i}\overline{y}(0).$

######  Theorem 2.

Consider the rank-$1$ network SIR epidemic model (15a)–(15b). Then, for every
$i=1,\ldots,n$,

  1. (i)

if

$\dot{y}_{i}(0)\leq 0\,,$ (33)

and

$w_{i}(0)\leq 0\,,$ (34)

then $y_{i}(t)$ is strictly decreasing for $t\geq 0$;

  2. (ii)

if

$\dot{y}_{i}(0)>0\,,$ (35)

or if

$\dot{y}_{i}(0)=0\,,$ (36)

and

$w_{i}(0)>0\,,$ (37)

then there exists a peak time $\hat{t}_{i}>0$ such that $y_{i}(t)$ is strictly
increasing on $[0,\hat{t}_{i}]$ and strictly decreasing on
$[\hat{t}_{i},+\infty)$;

  3. (iii)

if

$\dot{y}_{i}(0)<0\,,$ (38)

and

$w_{i}(0)>0\,,$ (39)

then either $y_{i}(t)$ is strictly decreasing for $t\geq 0$ or there exist a
local minimum time $\check{t}_{i}$ and a peak time $\hat{t}_{i}$ such that
$0<\check{t}_{i}<\hat{t}_{i}$ and $y_{i}(t)$ is strictly decreasing on
$[0,\check{t}_{i}]$, strictly increasing on $[\check{t}_{i},\hat{t}_{i}]$, and
strictly decreasing on $[\hat{t}_{i},+\infty)$;

###### Proof.

(i) If (34) holds true, then $\overline{t}_{i}=0$. On the other hand, (33) and
Theorem 1(ii) rule out the possibility that there exist any minimum point for
$y_{i}(t)$. Therefore, $y_{i}(t)$ is strictly decreasing for $t\geq 0$.

(ii) If condition (35) holds true, then Theorem 1(i) implies that
$\check{t}_{i}=0$ is the only minimum point of $y_{i}(t)$. On the other hand,
if equation (36) and condition (39) both hold true, then it follows from (20)
that

$\ddot{y}_{i}(0)=a_{i}x_{i}(0)\overline{y}(0)w_{i}(0)-\gamma\dot{y}_{i}(0)=\gamma\overline{y}^{2}(0)w_{i}(0)>0\,,$

(where the second identity follows from the fact that, by (15b),
$a_{i}x_{i}(t)\bar{y}(t)=\gamma y_{i}(t)$ when $\dot{y}_{i}(t)=0$) thus
implying that also in this case $\check{t}_{i}=0$ is a local minimum point for
$y_{i}(t)$.

Since, by Proposition 1(vi),

$\lim_{t\to+\infty}y_{i}(t)=0<y_{i}(0)\,,$ (40)

and, by Theorem 1(i), $y_{i}(t)$ cannot have another local minimum points
besides $\check{t}_{i}=0$, it follows that exists a peak time $\hat{t}_{i}>0$
such that $y_{i}(t)$ is strictly increasing for $t$ in $[0,\hat{t}_{i}]$ and
strictly decreasing for $t$ in $[\hat{t}_{i},+\infty)$.

(iii) From condition (38), $0$ is a nonstationary local maximum point of
$y_{i}(t)$. Since, by Theorem 1(i), $y_{i}(t)$ can have at most one local
minimum point, and (40) holds true, it follows that either $y_{i}(t)$ is
strictly decreasing for $t\geq 0$ (in case there is no local minimum point)
or, if a local minimum point $\check{t}_{i}>0$ exists, then there exists also
a peak time $\hat{t}_{i}>\check{t}_{i}$ so that $y_{i}(t)$ is strictly
increasing on $[0,\hat{t}_{i}]$ and strictly decreasing on
$[\hat{t}_{i},+\infty)$. ∎

Figure 2: Conceptual illustration of Theorem 2.

The previous result provides a classification of the behavior of the fraction
of infected individuals in the single populations and in Figure 2 there is a
conceptual illustration of this result. In particular, observe that if
$\dot{y}_{i}(0)=0$, then the behavior of the single infected fraction depends
only on the sign of the quantity $w_{i}(0)$, thus, in this case, the Theorem 2
provides a tight condition.

In Theorem 2, each condition is considered from the perspective of the
individual population $i$. However, note that if (33) is true for all
$i=1,\ldots,n$, which means that the infected curve in each population starts
decreasing, then no $y_{i}$ will have a local minimum point. Indeed, suppose
to multiply both sides of (33) by $b_{i}$ for all $i=1,\ldots,n$ and then sum
over all populations, this implies $\tilde{x}(0)<\gamma$ and therefore
$w_{i}(0)<0$ for all $i=1,\ldots,n$.

##  4 Sufficient conditions for multimodal behaviors

We now consider a particular class of rank-1 interaction matrices in the form

$A=\beta\bm{1}b^{T}\,,$ (41)

with $\beta>0$ and $\bm{1}^{T}b=1$. This is a special case of the one studied
in Section 3 where the vector $a$ has all equal entries and the entries of $b$
sum up to $1$.

######  Remark 4.

This model corresponds to a scenario in which all individuals have the same
susceptibility to the disease but different capabilities of spreading the
disease. For example, individuals wearing medical masks become infected with
the same probability but spread the disease differently. Note that a simple
case of this class of matrices is $A=\bm{1}\bm{1}^{\prime}$, that is studied
in Example 1. The network SIR epidemic model with this interaction matrix is
of interest for control applications, as analyzed in [acemoglu2021optimal].
Indeed, even if the dynamics at the nodes are homogeneous and thus the
infection spreads at the same rate, it may be convenient to divide agents into
multiple groups, for example, to study the effects of differentiated control
policies, especially in cases whereby the cost of applying a control and
epidemic cost for the diffusion of the disease may differ depending on the age
of the agents.

We observe that, for rank-1 interaction matrices in the form (41), the
dynamics become

$\dot{x}_{i}=-\beta x_{i}\bar{y},\qquad\dot{y}_{i}=\beta x_{i}\bar{y}-\gamma
y_{i}\,,$ (42)

for every $i=1,\cdots,n$, and

$\dot{\overline{x}}=-\beta\overline{x}\,\overline{y},\qquad\dot{\overline{y}}=\overline{y}\left(\beta\overline{x}-\gamma\right)\,,$
(43)

since $\overline{x}$ and $\tilde{x}$ differ in a constant term only. The next
result provides sufficient conditions for multimodal behavior of the infection
curve at single node level and encompasses Example 1. We first need to define
auxiliary functions

$g_{i}(\varepsilon)=\frac{1-\varepsilon}{1-b_{i}\varepsilon}\left(1-\frac{\gamma}{\beta}+\frac{\gamma}{\beta}\log\frac{\gamma}{\beta(1-b_{i}\varepsilon)}\right)-\varepsilon\,,$
(44)

Notice that

$g_{i}(0)=1-\frac{\gamma}{\beta}+\frac{\gamma}{\beta}\log\frac{\gamma}{\beta},\quad
g_{i}(1)=-1$

As a consequence, when $\gamma/\beta<1$, $g_{i}$ admits zeroes in $[0,1]$ and
we put

$\overline{\varepsilon}_{i}=\min\left\\{\varepsilon\in[0,1]:\,g_{i}(\varepsilon)=0\right\\}.$
(45)

######  Proposition 5.

Consider the rank $1$ network SIR epidemic model (15a)–(15b) with
$a=\beta\bm{1}$ and $\bm{1}^{T}b=1$. Consider an agent $i\in\\{1,\ldots,n\\}$
and an initial condition $(x(0),y(0))$ that satisfy the following conditions:

$\displaystyle x(0)+y(0)=\bm{1}$ (46) $\displaystyle\beta
x_{i}(0)\bar{y}(0)-\gamma y_{i}(0)<0$ (47)
$\displaystyle\beta\overline{x}(0)>\gamma$ (48) $\displaystyle
0<y_{i}(0)<\overline{\varepsilon}_{i}$ (49)

Then, there exist a local minimum time $\check{t}_{i}$ and a peak time
$\hat{t}_{i}$ such that $0<\check{t}_{i}<\hat{t}_{i}$ and $y_{i}(t)$ is
strictly decreasing on $[0,\check{t}_{i}]$, strictly increasing on
$[\check{t}_{i},\hat{t}_{i}]$, and strictly decreasing on
$[\hat{t}_{i},+\infty)$.

###### Proof.

From (47) and (42) it follows that $\dot{y}_{i}(0)<0$, which implies that
$y_{i}(t)$ is strictly decreasing for sufficiently small $t>0$. On the other
hand, (43) and (48) imply that $\bar{y}(t)$ is strictly increasing for
sufficiently small $t>0$. Since $\overline{x}$ and $\overline{y}$ satisfy the
scalar autonomous SIR epidemic model (43), this implies that $\overline{y}(t)$
has a peak at some time $\hat{t}>0$ and

$\overline{x}(\hat{t})=\frac{\gamma}{\beta}\,.$ (50)

From Proposition 2(i) we obtain that the aggregate peak of infection is

$\displaystyle\overline{y}(\hat{t}\,)$
$\displaystyle=\overline{x}(0)+\overline{y}(0)-\overline{x}(\hat{t})+\frac{\gamma}{\beta}\log\frac{\overline{x}(\hat{t})}{\overline{x}(0)}$
(51)
$\displaystyle=1-\frac{\gamma}{\beta}+\frac{\gamma}{\beta}\log\frac{\gamma}{\beta\overline{x}(0)}\,,$

where the second equality follows from (46) and (50). Moreover, (42) and (43)
imply that

$\frac{x_{i}(t)}{x_{i}(0)}=\frac{\overline{x}(t)}{\overline{x}(0)}\,,$ (52)

for every $i=1,\dots,n,$ and every $t$. Therefore, using (50) we obtain

$x_{i}(\hat{t})=\frac{\gamma}{\beta}\frac{x_{i}(0)}{\overline{x}(0)}\,,$ (53)

We now prove that $y_{i}(t)$ cannot remain decreasing for all $t>0$. Assume by
contradiction that

$\dot{y}_{i}(t)\leq 0\,,\qquad\forall t\in[0,\hat{t}]\,.$ (54)

In particular, this implies that $y_{i}(\hat{t})\leq y_{i}(0)$. This together
with (53) and (51) imply that

$\begin{array}[]{rclcl}0&\geq&\dot{y}_{i}(\hat{t}\,)\\\\[7.0pt]
&=&\displaystyle\beta x_{i}(\hat{t}\,)\bar{y}(\hat{t}\,)-\gamma
y_{i}(\hat{t}\,)\\\\[7.0pt]
&=&\displaystyle\gamma\frac{x_{i}(0)}{\overline{x}(0)}\left(1-\frac{\gamma}{\beta}+\frac{\gamma}{\beta}\log\frac{\gamma}{\beta\overline{x}(0)}\right)-\gamma
y_{i}(\hat{t}\,)\\\\[9.0pt]
&\geq&\displaystyle\gamma\left[\frac{x_{i}(0)}{\overline{x}(0)}\left(1-\frac{\gamma}{\beta}+\frac{\gamma}{\beta}\log\frac{\gamma}{\beta\overline{x}(0)}\right)-y_{i}(0)\right]\end{array}$
(55)

Notice now that, because of the assumptions on $b$, it holds that
$\overline{x}(0)=1-\overline{y}(0)\leq 1-b_{i}y_{i}(0)$. Since the last
expression in (55) is decreasing in $\overline{x}(0)$ and
$x_{i}(0)=1-y_{i}(0)$, we obtain that

$0\geq\displaystyle\gamma g_{i}(y_{i}(0))$ (56)

Notice that by (48), we necessarily have that $\gamma/\beta<1$. This implies
that $g_{i}(0)>0$ and, together with (56), implies that
$y_{i}(0)\geq\overline{\varepsilon}_{i}$, thus violating (45). This
contradiction implies that $y_{i}(t)$ cannot remain decreasing for all $t>0$.
The thesis then follows from Theorem 2. ∎

######  Remark 5.

Observe that the set of model parameters and initial conditions that satisfy
the assumptions of Proposition 5 is nonempty. To prove this, let us consider a
network with $n$ nodes, interaction matrix as defined in (41) with parameters
$\beta>\gamma$ and

$b_{1}<\min\left\\{\frac{\gamma}{\beta},1-\frac{\gamma}{\beta}\right\\}\,,$
(57)

Fix an initial condition such that $y_{1}(0)\in(0,\overline{\epsilon}_{1})$
and $y_{j}(0)=0$ for every $j=2,\cdots,n$. Notice that
$\overline{y}(0)=b_{1}y_{1}(0)$. A straightforward check shows that (47) and
(48) are automatically satisfied putting no further restriction on $y_{1}(0)$.
Therefore, all assumptions of Proposition 5 are satisfied.

######  Remark 6.

Observe that under the assumptions of Proposition 5 we can provide an upper
bound for stationary infection peaks of a node $i$. Let $\hat{t}_{i}$ be the
peak time of node $i$, hence $\dot{y}_{i}(\hat{t}_{i})=0$. This implies by
(42)

$y_{i}(\hat{t}_{i})=\frac{\beta
x_{i}(\hat{t}_{i})}{\gamma}\overline{y}(\hat{t}_{i})\,.$ (58)

Since the aggregate quantity of infected is limited above by its infection
peak, i.e. $\overline{y}(t)\leq\overline{y}(\hat{t}\,)$, for all $t\geq 0$,
and the fraction of susceptibles is monotonically decreasing, then

$\begin{array}[]{rclcl}y_{i}(\hat{t}_{i})&\leq&\displaystyle\frac{\beta
x_{i}(0)}{\gamma}\overline{y}(\hat{t})\\\ &=&\displaystyle\frac{\beta
x_{i}(0)}{\gamma}\left(\overline{x}(0)+\overline{y}(0)-\overline{x}(\hat{t})+\gamma\log\frac{\overline{x}(\hat{t}\,)}{\overline{x}(0)}\right)\\\\[9.0pt]
&=&\displaystyle
x_{i}(0)\left(\frac{\beta}{\gamma}-1+\log\frac{\gamma}{\overline{x}(0)}\right)\,,\\\
\end{array}$

where the first equivalence follows from (51) and the last one from (50) and
(46).

##  5 Numerical Simulations

In this section we provide numerical simulations of the network SIR model. We
start by considering a network with $n=5$ nodes and rank-1 interaction matrix.
We can observe from Figure 3 that nodes $1$ and $3$ exhibit multimodal
behaviors with two changes of monotonicity: for each $i=1,3$, the fraction of
infected $y_{i}(t)$ is strictly decreasing for times $t$ in
$[0,\check{t}_{i}]$, it is strictly increasing in
$[\check{t}_{i},\hat{t}_{i}]$, and it is strictly decreasing for $t$ in
$[\hat{t}_{i},+\infty)$. The bottom plot shows the aggregate variable
$\overline{y}(t)$, which has a unimodal behavior with a peak in $\hat{t}$, as
proved in the Proposition 4. Note also that, as consequence of Theorem 1(ii)
(cfr. Remark 3), $\check{t}_{i}\leq\hat{t}$ for every $i=1,3$, namely the
local minimum of the infected curve in each node cannot occur after the
aggregate infection peak.

Figure 3: Numerical simulations of the network SIR model on a network with
$n=5$ nodes and rank-1 interaction matrix. The interaction matrix is
$A=ab^{T}$ with $a=(0.1\,0.25\,0.6\,1\,0.2)$ and
$b=(0.45\,0.4\,0.6\,0.65\,0.01)$, and the recovery rate is $\gamma=0.6$. The
initial conditions are $x_{1}(0)=0.85$, $x_{2}(0)=0.999$, $x_{3}(0)=0.8$,
$x_{4}(0)=1$, $x_{5}(0)=0.75$, and $y(0)=\bm{1}-x(0)$.

Figure 4 illustrates numerical simulations of a network SIR model with $n=2$
and full rank interaction matrix. This simulation shows that dynamics may
exhibit more complex behaviors than in the rank $1$ case. In particular, the
bottom plot of Figure 4 shows that for suitable initial conditions, the curve
of infected in node 1 exhibits two stationary local maxima (with the second
higher than the first one) with three changes of monotonicity. We remark that
for rank-1 interaction matrices this behavior is ruled out by Theorem 1(iii).

Figure 4: Numerical simulations of the network SIR model on networks with
$n=2$ nodes and full rank interaction matrices. In the top plot,
$\displaystyle\hskip 85.35826ptA=\left(\begin{array}[]{cc}0.8&0.5\\\
0.1&1.5\end{array}\right)\,,$
$x_{1}(0)=0.92=1-y_{1}(0)$, $x_{2}(0)=0.99=1-y_{2}(0)$, and $\gamma=0.9$. In
the bottom plot,
$\displaystyle\hskip 85.35826ptA=\left(\begin{array}[]{cc}0.95&0.2\\\
0.001&1.6\end{array}\right)\,,$
$x_{1}(0)=0.98=1-y_{1}(0)$, $x_{2}(0)=1=1-y_{2}(0)$, and $\gamma=0.9$.

Figure 5:  Numerical simulation of the network SIR epidemic model with $n=4$
nodes, recovery rate $\gamma=0.5$, and interaction matrix
$\displaystyle\hskip
28.45274ptA=\left(\begin{array}[]{cccc}0.05&0.05&0.05&0.05\\\
0.0001&0.8&0.0001&0.0001\\\ 0.0001&0.0001&0.05&0.0001\\\
0.0009&0.0009&0.0009&0.9\end{array}\right)\,.$
Initial conditions are $x_{1}(0)=x_{2}(0)=1$, $y_{1}(0)=y_{2}(0)=0$,
$x_{3}(0)=0.7=1-y_{3}(0)$, $x_{4}(0)=0.99=1-y_{4}(0)$. It can be observed that
the curve of infected individuals in the first node has three peaks.

Figure 6: Numerical simulation of the network SIR epidemic model with $n=4$
nodes, recovery rate $\gamma=0.5$, and interaction matrix
$\displaystyle\hskip
28.45274ptA=\left(\begin{array}[]{cccc}0.05&0.07&0.05&0.05\\\
0.0001&0.8&0.0001&0.0001\\\ 0.0001&0.0001&0.1&0.0001\\\
0.01&0.01&0.01&0.9\end{array}\right)\,.$
Initial conditions are $x_{1}(0)=x_{2}(0)=x_{4}(0)=1$,
$y_{1}(0)=y_{2}(0)=y_{4}(0)=0$, $x_{3}(0)=0.9=1-y_{3}(0)$. It can be observed
that the curve of infected individuals in the first subpopulation has three
peaks, and the second peak is the highest one.

In Figure 5 we illustrate a numerical simulation of a network SIR model with
$n=4$ nodes, where the curves of infected at single node level exhibit three
peaks. Moreover, we observe the presence of delays among the peaks in
different nodes. Such simulations show that the limitation of the number of
peaks to two is a peculiar feature of rank $1$ matrices while for general
network matrices, even in correspondence of a limited number of nodes, such
limitation does no longer hold. It is also important to note that successive
infection peaks can also be greater than the previous one, as observed in
Figure 6. This is an interesting phenomenon especially in the design of a
control, as an epidemic process is simulated in which after one peak of
infection a second larger one may occur.

##  6 Conclusion

In this paper, we studied the network epidemic SIR model and provided
theoretical results on the dynamics for the special case of rank 1 interaction
matrices, which finds applications in epidemics control, as shown in
[acemoglu2021optimal]. We first proved that, in contrast to the scalar SIR
model, in the network SIR model the curve of infected individuals in a single
node may exhibit multimodal behaviors, and established sufficient conditions
for the occurrence of this phenomenon. We then provided a linear combination
of the fraction of infected in each node that exhibits a unimodal behavior,
and characterized all the possible behaviors that the dynamics may exhibit at
single node level, showing that the infection curve in a single node can
undergo two changes of monotonicity at most. We then conducted a numerical
analysis showing that for more general interaction matrices the network SIR
model may exhibit more than two peaks at single node level. Future work aims
to include in the model more complex phenomena, such as waving immunity and
endogenous response of agents in response to the disease, to fully describe
the occurrence of multimodal behaviors in epidemic models.

##  7 Proof of Lemma 1

By taking the derivative of both sides of (13) and substituting (15a), we get

$\dot{\overline{x}}=\sum_{j=1}^{n}b_{j}\dot{x}_{j}=-\bar{y}\sum_{j=1}^{n}b_{j}a_{j}x_{j}=-\overline{y}\tilde{x}\,,$

thus proving (18). Analogously, by taking the derivative of both sides of (14)
and substituting (15b), we get

$\dot{\overline{y}}=\sum_{j=1}^{n}b_{j}\dot{y}_{j}=\bar{y}\sum_{j=1}^{n}b_{j}a_{j}x_{j}-\gamma\overline{y}=\overline{y}\left(\tilde{x}-\gamma\right)\,,$

thus proving (19). On the other hand, taking the derivative of both sides of
(15b) and substituting (15a) and (19) yield

$\begin{array}[]{rcl}\ddot{y}_{i}&=&a_{i}\left(\dot{x}_{i}\bar{y}+x_{i}\dot{\bar{y}}\right)-\gamma\dot{y}_{i}\\\\[3.0pt]
&=&a_{i}x_{i}\overline{y}\left(\tilde{x}-\gamma-
a_{i}\overline{y}\right)-\gamma\dot{y}_{i}\\\\[3.0pt]
&=&a_{i}x_{i}\overline{y}w_{i}-\gamma\dot{y}_{i}\,,\end{array}$

thus proving (20).

##  8 Proof of Lemma 2

(i) From (15a), we have that

$\dot{\tilde{x}}=\sum_{i}a_{i}b_{i}\dot{x}_{i}=-\overline{y}\sum_{i}a_{i}^{2}b_{i}x_{i}\leq
0\,.$ (59)

On the other hand, thanks to Proposition 1(v), the assumption that $y(0)\gneq
0$ implies that

$y(t)>0\,,\qquad\forall t>0\,,$ (60)

so that in particular

$\overline{y}(t)=\sum_{j}b_{j}y_{j}(t)>0\,,\qquad\forall t\geq 0\,.$ (61)

It then follows from (17),(19), (59), and (61) that

$\dot{w}_{i}=\dot{\tilde{x}}-a_{i}\dot{\overline{y}}=\dot{\tilde{x}}-a_{i}\overline{y}(\tilde{x}-\gamma)=\dot{\tilde{x}}-a_{i}\overline{y}w_{i}-a_{i}^{2}\overline{y}^{2}<-a_{i}\overline{y}w_{i}\,,$

thus proving the claim.

(ii) For $0\leq t\leq\overline{t}_{i}$, we have $w_{i}(t)\geq 0$ so that, by
point (i) and (61), we have that

$\dot{w}_{i}(t)<-a_{i}\overline{y}(t)w_{i}(t)\leq 0\,.$

The above implies that $w_{i}(t)$ is strictly decreasing for $0\leq
t\leq\overline{t}_{i}$.

(iii) It follows from point (i) that, for every $t_{i}^{*}\geq 0$ such that
$w_{i}(t_{i}^{*})=0$, we have

$\dot{w}_{i}(t_{i}^{*})<-a_{i}\overline{y}(t_{i}^{*})w_{i}(t_{i}^{*})=0\,,$

so that $w_{i}(t)<0$ for $t>t_{i}^{*}$. If $w_{i}(0)\leq 0$ so that
$\overline{t}_{i}=0$, the above implies that $w_{i}(t)<0$ for every
$t>0=\overline{t}_{i}$. On the other hand, $w_{i}(0)>0$ so that
$\overline{t}_{i}>0$, the above implies that $w_{i}(t)<0$ for every
$t>\overline{t}_{i}$.

(iv) Assume by contradiction that $t\geq\overline{t}_{i}$ is a local minimum
point of $y_{i}(t)$ with $\dot{y}_{i}(t)=0$. By (20), we then have

$\ddot{y}_{i}(t)=a_{i}x_{i}(t)\overline{y}(t)w_{i}(t)-\gamma\dot{y}_{i}(t)\\\
=\gamma y_{i}(t)w_{i}(t)\,,$ (62)

where the last identity holds true since, by (15b), $\dot{y}_{i}(t)=0$ is
equivalent to

$a_{i}x_{i}(t)\bar{y}(t)=\gamma y_{i}(t)\,.$ (63)

Equation (62) and (60) imply that

$\operatorname{sgn}(\ddot{y}_{i}(t))=\operatorname{sgn}(w_{i}(t))\,.$ (64)

Now, recall that by point (iii) we have $w_{i}(t)\leq 0$ for every
$t\geq\overline{t}_{i}$. If $w_{i}(t)<0$, then equation (64) implies that
$\ddot{y}_{i}(t)<0$, so that $t$ cannot be a local minimum point for
$y_{i}(t)$. On the other hand, if $w_{i}(t)=0$, by point (i) and (61) we get
that

$\dot{w}_{i}(t)<-a_{i}\overline{y}(t)w_{i}(t)=0\,,$ (65)

while equation (64) implies that $\ddot{y}_{i}(t)=0$. Hence, by taking the
derivative of both sides of (20), by (63) and $w_{i}(t)=0$, we get

$\begin{array}[]{rcl}\dddot{y}_{i}(t)&=&a_{i}x_{i}(t)\overline{y}(t)\dot{w}_{i}(t)+a_{i}(\dot{(x_{i}(t)\overline{y}(t))})w_{i}(t)-\gamma\ddot{y}_{i}(t)\\\\[5.0pt]
&=&\gamma y_{i}(t)\dot{w}_{i}(t)\\\\[5.0pt] &<&0\,,\end{array}$

where the last inequality following from (60) and (65). Together with
$\dot{y}_{i}(t)=\ddot{y}_{i}(t)=0$, the above implies that $t$ is an
inflection point for $y_{i}(t)$, hence it particular it is not a local minimum
point. We have thus proven that $y_{i}(t)$ cannot have any stationary minimum
points in the interval $[\overline{t}_{i},+\infty)$.

[]Martina Alutto received the B.Sc. in Mathematics for Engineering in 2018
from Politecnico di Torino, Italy, the M.S. in Mathematical Engineering (cum
Laude) in 2021 from Politecnico of Turin, Italy. She is currently a PhD
student in Pure and Applied mathematics at the Department of Mathematical
Sciences, Politecnico di Torino, Italy. Her research interests focuses on
analysis and control of network systems, with application to epidemics and
social networks.

[]Leonardo Cianfanelli received the B.Sc. (cum Laude) in Physics and
Astrophysics in 2014 from Università degli Studi di Firenze, Italy, the M.S.
in Physics of Complex Systems (cum Laude) in 2017 from Università di Torino,
Italy, and the PhD in Pure and Applied Mathematics in 2022 from Politecnico di
Torino. He is currently a Research Assistant at the Department of Mathematical
Sciences, Politecnico di Torino, Italy. He was a Visiting Student at
Laboratory for Information and Decision Systems, Massachusetts Institute of
Technology, in 2018-2020. His research focuses on control in network systems,
with applications to transportation networks and epidemics.

[] Giacomo Como(M’12) is a Professor at the Department of Mathematical
Sciences, Politecnico di Torino, Italy. He is also a Senior Lecturer at the
Automatic Control Department, Lund University, Sweden. He received the B.Sc.,
M.S., and Ph.D. degrees in Applied Mathematics from Politecnico di Torino, in
2002, 2004, and 2008, respectively. He was a Visiting Assistant in Research at
Yale University in 2006–2007 and a Postdoctoral Associate at the Laboratory
for Information and Decision Systems, Massachusetts Institute of Technology in
2008–2011. Prof. Como currently serves as Senior Editor for the IEEE
Transactions on Control of Network Systems, as Associate Editor for
Automatica, and as the chair of the IEEE-CSS Technical Committee on Networks
and Communications. He served as Associate Editor for the IEEE Transactions on
Network Science and Engineering (2015-2021) and for the IEEE Transactions on
Control of Network Systems (2016-2022). He was the IPC chair of the IFAC
Workshop NecSys’15 and a semiplenary speaker at the International Symposium
MTNS’16. He is recipient of the 2015 George S. Axelby Outstanding Paper Award.
His research interests are in dynamics, information, and control in network
systems, with applications to infrastructure, social, and economic networks.

[] Fabio Fagnani received the Laurea degree in Mathematics from the University
of Pisa and the Scuola Normale Superiore, Pisa, Italy, in 1986. He received
the PhD degree in Mathematics from the University of Groningen, Groningen, The
Netherlands, in 1991. From 1991 to 1998, he was an Assistant Professor at the
Scuola Normale Superiore. In 1997, he was a Visiting Professor at the
Massachusetts Institute of Technology. Since 1998, he has been with the
Politecnico of Torino, where since 2002 he has been a Full Professor. From
2006 to 2012, he has acted as Coordinator of the PhD program in Mathematics
for Engineering Sciences at Politecnico di Torino. From 2012 to 2019, he
served as the Head of the Department of Mathematical Sciences, Politecnico di
Torino. His current research focuses on dynamical systems over graphs,
inferential distributed algorithms, and opinion dynamics. He is an Associate
Editor of the IEEE Transactions on Automatic Control and served in the same
role for the IEEE Transactions on Network Science and Engineering and of the
IEEE Transactions on Control of Network Systems.


# Evaluation of cosmogenic production of ${}^{39}Ar$ and ${}^{42}Ar$ for rare-
event physics using underground argon

C. Zhang<EMAIL_ADDRESS>Division of Natural Sciences, Mount Marty
University, Yankton, SD 57078  D.-M. Mei<EMAIL_ADDRESS>Department of
Physics, The University of South Dakota, Vermillion, SD 57069

###### Abstract

Underground argon (UAr) with lower cosmogenic activities of ${}^{39}Ar$ and
${}^{42}Ar$ has been planned as a detector in detecting scintillation light
and charge collection using time projection chambers for dark matter searches
and as a veto detector in suppressing backgrounds for neutrinoless double beta
decay (0$\nu\beta\beta$) experiments. Long-lived radioactive isotopes,
${}^{39}Ar$ and ${}^{42}Ar$, can also be produced on the surface when UAr is
pumped out from a deep well. Understanding the production of long-lived
isotopes in Ar is important for utilizing UAr for dark matter and
0$\nu\beta\beta$ experiments in terms of its production, transportation, and
storage. Ar exposure to cosmic rays at sea-level is simulated using Geant4 for
a given cosmic ray muon, neutron, and proton energy spectrum. We report the
simulated cosmogenic production rates of ${}^{39}Ar$, ${}^{42}Ar$, and other
long-lived isotopes at sea-level from fast neutrons, high energy muons, and
high energy protons. Total production rates of 938.53/kg${}_{Ar}\cdot$day and
5.81$\times$10-3/kg${}_{Ar}\cdot$day for 39Ar and 42Ar are found from our
simulation. Utilizing these production rates, we set a time limit of 954 days
constrained by the production of 39Ar for UAr to be on the surface before it
compromises the sensitivity for a dark matter experiment. Similarly, a time
limit of 1702 days constrained by the production of 42Ar is found for a
0$\nu\beta\beta$ experiment.

###### keywords:

Cosmogenic activation , Underground argon , Geant4 simulation

###### PACS:

13.85.Tp , 23.40-s , 25.40.Sc , 28.41.Qb , 95.35.+d , 29.40.Mc

,

##  1 Introduction

Observational evidence indicates that dark matter accounts for approximately
85% of the matter in our universe[1, 2, 3]. But the nature of dark matter is
still mysterious. Of all the dark matter candidates, the Weakly Interacting
Massive Particle (WIMP) seems promising and is favored by scientific
communities. Assuming the WIMP is a new elementary particle, which is
ubiquitous in galaxies, it could be directly detected through its elastic
recoils on ordinary matter[4, 5]. Many dark matter experiments are actively
seeking WIMP scattering events but none have succeeded so far [6, 7, 8, 9, 10,
11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25]. DAMA/LIBRA has
been claiming for decades that the observation of the annual modulation in the
detection is the recoil signal stroked by WIMPs [26], although their claim is
not confirmed by any other experiments. To further push dark matter detection
sensitivities, the next generation large-scale dark matter experiments are
being built aiming at the WIMP-nucleon cross section down to the level of
$10^{-48}\,\rm{cm}^{2}$ [27, 28, 29]. This level of sensitivity also reaches
the floor of the neutrino background due to coherent neutrino induced nuclear
recoils in the detector [30, 31]. In order to achieve such detection
sensitivity, large-scale detectors with hundreds of ton-years exposure under
an ultra-low background environment are needed. For example, DarkSide-20k,
which will use underground argon (UAr) as the target [29], is a two-phase
liquid Ar detector with a 23 tonnes active volume.

The experiments seeking to measure the half-life of neutrinoless double-beta
decay ($0\nu\beta\beta$) aim to determine the Majorana nature of the neutrino
and help understand the absolute neutrino masses and their mass hierarchy[32,
33, 34, 35]. No discovery has been made by current $0\nu\beta\beta$
experiments with the sensitivity of the decay half-life up to $\sim 10^{26}$
years [36]. This sets an upper limit on the effective Majorana mass of the
electron neutrino less than 0.2 eV[36, 37, 38, 39, 40]. A future large-scale
experiment, LEGEND-1000, with increased sensitivity up to $10^{28}$ years has
been proposed [41]. This aims for a sensitivity of 0.01 event per ton per year
in the region of interest. Aiming at such detection sensitivity, instrumental
backgrounds from radioactive isotopes need to be well suppressed and
understood. LEGEND-1000 plans to use UAr as the active veto detector.

Both dark matter and $0\nu\beta\beta$ experiments require their background
rate at the region of interest to be extremely low. The instrumental
background from long-lived radioactive isotopes in detector components must be
minimized and accurately measured. In addition, cosmogenic activation can add
more radioactivities to the background budget. Liquid Ar, a relatively cheap
noble liquid scintillator, is a widely used medium for the detection of
ionizing radiation. Taking advantage of its high ionization yield and
characteristic scintillation time profile [42], liquid Ar is utilized as an
active target or scintillation veto material in dark matter searches and
$0\nu\beta\beta$ experiments [12, 25, 29, 36, 41]. The problematic issue is
that there are long-lived radioactive isotopes existing in Ar when it is
produced on the surface. Both 39Ar (t1/2=269 years) and 42Ar (t1/2=32.9 years)
are produced through cosmogenic activation [43, 44]. The isotopic abundance of
39Ar in atmospheric argon is at a level of (8.1$\pm$0.6) $\times$10-16g/g,
which results in a decay rate of $\sim$1 Bq/kg in a liquid Ar detector [12,
25]. It is a pure $\beta$ decay to a stable daughter nucleus (39K). This decay
rate and the $\beta$ decay Q value of 565.5 keV produce background and pile-up
concerns in the detectors for dark matter searches. On the other hand,
${}^{42}Ar$ was measured by the Gerda collaboration and the isotopic abundance
is determined at a level of (9.2${}^{+2.2}_{-4.6}$)$\times$10-21 g/g. More
recent measurement from the DEAP collaboration shows the specific activity of
42Ar to be 40.4 $\pm$ 5.9 $\mu$Bq/kg [46] which is lower than the GERDA’s
result of (88${}^{+21}_{-44})$ $\mu$Bq/kg. Its decay product, 42K, a $\beta$
emitter with a decay Q value of 3.52 MeV, is a background to 0$\nu\beta\beta$
decay experiments when using liquid argon as a veto detector [36]. Note that a
fraction of 42Ar could be produced in the Earth’s atmosphere due to nuclear
bomb tests from 1945 to 1962 [43]. The concentration of 42Ar after bomb
testing is estimated to be $\sim$ 10-22 to $\sim$10-23 atoms per natAr atoms
[43]. This is only $\sim$1% to $\sim$10% of the measured concentration of 42Ar
[46, 43]. The current 42Ar in Earth’s atmosphere is dominated by spallation
reactions from the cosmic ray nucleon component through
40Ar($\alpha,2p)^{42}$Ar process.

To obtain Ar with a lower level of 39Ar and 42Ar, it is natural to explore Ar
which has been deep underground where Ar has existed for thousands of years
[45]. Since the production of 39Ar depends on the depth of underground wells
and the surrounding rock condition in terms of its porosity, UAr may not
always have a low level of 39Ar [45]. Because the production of 42Ar is mainly
through 40Ar($\alpha$,2p)42Ar and the interaction threshold energy of $\alpha$
particles is above 12 MeV, it is expected that the production of 42Ar is
strongly suppressed in deep underground since there are no $\alpha$ particles
with energy greater than 12 MeV from natural radioactivity decays. Amazingly,
the UAr with low radioactivity produced by the DarkSide collaboration [47] in
Colorado of the United States has shown a reduction of 1400 in 39Ar when
compared to atmospheric Ar. Note that the UAr produced elsewhere may not have
such a low level of 39Ar [45]. Although there is no measured reduction factor
of 42Ar reported from the UAr produced by DarkSide, it is expected that the
reduction of 42Ar should be at least a factor of 1400 [41]. Therefore, the UAr
produced by DarkSide in Colorado is particularly valuable for DarkSide-20k and
LEGEND-1000. To keep low level of 39Ar and 42Ar for underground experiments,
any potential production processes should be carefully examined and suppressed
since 39Ar and 42Ar are mainly produced through cosmogenic processes when Ar
is on the surface [43, 44], precise calculations of cosmogenic activation are
needed to determine the maximum tolerable surface exposure time of argon
during its production, transport and storage.

In this paper, we evaluate cosmogenic production of radioactive isotopes in
Ar. The production rates of ${}^{39}Ar$, ${}^{42}Ar$, and other long-lived
isotopes from fast neutrons, high energy muons, and high energy protons are
obtained using Geant4-based simulations. The results are also compared with
some experimental data.

##  2 Evaluation of cosmogenic production in Ar

###  2.1 Evaluation tools and input sources

GEANT4 (V10.7p02) with shielding modular physics list [48, 49] is used for
this study. It includes a set of electromagnetic and hadronic physics
processes for high energy or high precision simulation needs.
G4MuonNuclearProcess is activated to simulate the muon-nuclear interactions.

Cosmic ray muon flux at sea level can be parameterized using Gaisser’s
formula[50]. It assumes a flat Earth which is only valid for muons with energy
$E_{\mu}>100GeV$ (see Fig. 1). In our simulation, the input muons are sampled
by using the modified Gaisser’s formula[51]. The total muon flux is normalized
to 0.015 $cm^{-2}s^{-1}$ with energy spanning from 1 GeV to 100 TeV. Stop
muons are not included in the simulation due to their relatively small
contribution to the total cosmogenic activation budget for Ar.

Figure 1:  Vertical muon flux at sea level. The dashed line is given by
Gaisser’s formula[50]. The solid line stands for the modified Gaisser’s
formula[51]. As a comparison, the measured muon distribution at $0^{\circ}$
zenith angle from L3+C[52] and BESS[53] are also presented.

The input fast neutrons with energy greater than 4 MeV are adopted from the
ground-based neutron flux measurements at a reference location at sea level,
New York City (NY data)[54]. The total flux is normalized to be 0.004
$cm^{-2}s^{-1}$. It is worth mentioning that the MeV neutrons may slightly
vary from location to location due to local radioactivities such as neutron
yield from $(\alpha,n)$ process in materials.

Cosmic ray proton induced activation in Ar is also simulated although the
number of protons below a few GeV is much less than the number of neutrons in
the atmosphere (see Fig. 2). The proton spectrum is harder than the neutron
one with energies exceeding a few GeV. The input cosmic ray proton flux is
adopted from the result of MCNPX [55] simulation code. Protons with energy
$E_{p}>10$ MeV are sampled for our use.

Figure 2:  Comparison of neutron and proton energy distribution at sea level.
Black dots are experimental neutron data measured at a reference location at
New York City (NY data)[54]. Blue squares are Monte Carlo generated proton
spectrum using MCNPX code[55] with its flux normalized to measured proton data
(red triangles[56]).

###  2.2 Cosmogenic activation in Ar

Commercial Ar is mostly extracted from the Earth’s atmosphere. It consists
primarily of the stable isotopes ${}^{40}Ar$, ${}^{36}Ar$, and ${}^{38}Ar$.
Although natural argon is stable, it suffers from cosmic-ray bombardment at
the earth’s surface. Ar radionuclides produced by cosmogenic activation, such
as ${}^{37}Ar$, ${}^{39}Ar$ and ${}^{42}Ar$, are irreducible and turn out to
be a significant source of background for low-background Ar-based radiation
detectors. To effectively suppress cosmogenic activation, switching from
atmospheric Ar to deep UAr had been proposed for DarkSide-20k and LEGEND-1000.
Significant production of ultra-low background Ar extracted and purified from
deep $CO_{2}$ wells has been successfully developed. It achieves a factor of
1400 reduction (or 0.73 $\rm{mBq/kg_{Ar}}$) over atmospheric argon in terms of
their ${}^{39}Ar$ concentration[57]. However, the cosmogenic production rates
are still critical and need to be accurately calculated to optimize the
surface exposure time during its production, shipping and storage. Note that
37Ar is not so dangerous for dark matter and 0$\nu\beta\beta$ decay
experiments due to its shorter half-life (35 days).

A comprehensive Geant4 (V10.7p02) simulation has been conducted to evaluate
the cosmogenic activation rate in an Ar target. Geant4 default material
”G4_lAr” (40Ar:99.604%; 36Ar:0.334%; 38Ar:0.063%) is selected and stored in a
thin cylinder tank that is 2.0 m in diameter and 2.0 m in height. The liquid
argon has the density of 1.40 $g/cm^{3}$ which gives the total mass of
$\sim$8.80 ton per cylinder. The cosmic-ray muons, neutrons and protons are
considered as inputs at the top of the Ar target.

The production of 39Ar in 40Ar is through the reaction of 40Ar(n,
2n${}^{{}^{\prime}}$)39Ar when high energy neutrons ($>$10 MeV) bombard the Ar
target. This reaction is expected to be seen in the Geant4 simulation. A main
purpose of this simulation is to determine if 42Ar can be seen in the Geant4
simulation since the production rate of 42Ar is extremely small. The long-
lived argon isotope ${}^{42}Ar$ undergoes a $\beta$ decay with a half-life of
32.9 years. The $\beta$ decay of its daughter isotope, ${}^{42}K$, has a
maximum electron energy of 3.52 MeV which becomes a vital background for
$0\nu\beta\beta$ experiments in the region of interest
($Q_{\beta\beta}=2.039MeV$ for $0\nu\beta\beta$ decay from 76Ge)[41] when
using liquid Ar as a veto detector. It is understood that ${}^{42}Ar$ is
mainly produced through the reaction channel ${}^{40}Ar(\alpha,2p)^{42}Ar$
with a Q-value of -12.77 MeV[43, 44]. This means that the $\alpha$ particles
that can generate 42Ar in 40Ar must have kinetic energies greater than 12.77
MeV.

In order to see the production of 42Ar in our simulation in a reasonable
computing time (one month) in the high precision computing cluster with 72
CPUs at the University of South Dakota, all particles with energy less than 12
MeV are killed (tracks stop and kill). Therefore, our results exclude all
capture processes. For example 36Ar(n,$\gamma$)37Ar is not found in our
simulation. We found that the cosmogenic production rate of ${}^{42}Ar$
results in a rate of $5.81\times 10^{-3}\,\rm{atoms/kg_{Ar}/day}$ in total
which corresponds to the activity of $3.88\times
10^{-12}\,\rm{Bq/kg_{Ar}/day}$. This shows that the Geant4 simulation is able
to predict the production of 42Ar. We summarize the three long-lived argon
isotopes in Table 1 and other long-lived isotopes in Table 2.

Table 1:  Cosmogenic activation rates of three argon isotopes: ${}^{37}Ar$, ${}^{39}Ar$ and ${}^{42}Ar$. The simulation results are also compared with the ones from measurements and estimations[58].  | ${}^{37}Ar$ | ${}^{39}Ar$ | ${}^{42}Ar$
---|---|---|---
| $\rm{atoms/kg_{Ar}/day}$
Neutrons (this work) | 176.01 | 857.73 | $4.60\times 10^{-3}$
Neutrons (measurement[58]) | $51.0\pm 7.4$ | $759\pm 128$ | -
Muons (this work) | 2.40 | 52.27 | $1.57\times 10^{-4}$
Muons (calculation[58]) | - | $172\pm 26$ | -
Protons (this work) | 6.20 | 28.53 | $1.05\times 10^{-3}$
Protons (calculation[58]) | $1.73\pm 0.35$ | $3.6\pm 2.2$ | -
Total (this work) | 184.61 | 938.53 | $5.81\times 10^{-3}$
Total (Ref.[58]) | $52.73\pm 7.75$ | $934.6\pm 156.2$ | -
Table 2:  Cosmogenic activation rates of other long-lived isotopes.  Isotope, Half Life | Neutron | Muon | Proton
---|---|---|---
| $\rm{atoms/kg_{Ar}/day}$
${}^{3}H$, 12.32y | $3.00\times 10^{-4}$ | $6.56\times 10^{-6}$ | $1.05\times 10^{-4}$
${}^{7}Be,1.387\times 10^{6}y$ | $3.43\times 10^{-3}$ | $6.47\times 10^{-3}$ | $9.62\times 10^{-3}$
${}^{10}Be,53.22d$ | $7.05\times 10^{-3}$ | $5.22\times 10^{-3}$ | $1.10\times 10^{-2}$
${}^{14}C,5.703\times 10^{3}y$ | 0.10 | $9.85\times 10^{-3}$ | $4.56\times 10^{-2}$
${}^{22}Na,2.602y$ | 0.37 | $2.35\times 10^{-2}$ | $9.92\times 10^{-2}$
${}^{26}Al,7.17\times 10^{5}y$ | 0.63 | $4.24\times 10^{-2}$ | 0.12
${}^{32}Si,153y$ | 7.00 | 0.14 | 0.47
${}^{32}P,14.268d$ | 15.7 | 0.38 | 1.05
${}^{33}P,25.3d$ | 22.5 | 0.42 | 1.32
${}^{35}S,87.37d$ | 74.5 | 1.66 | 3.18
${}^{36}Cl,3.01\times 10^{5}y$ | 75.5 | 1.23 | 3.32
${}^{40}K,1.248\times 10^{9}y$ | 1.80 | $5.86\times 10^{-2}$ | 0.56
${}^{41}Ca,9.94\times 10^{4}y$ | $1.85\times 10^{-3}$ | $1.02\times 10^{-4}$ | $5.68\times 10^{-4}$

As can be seen from the results shown in Table 1 and 2, the production rates
of long-lived isotopes in Ar is dictated by high energy neutrons. The
production rate of 39Ar by fast neutrons from this work is in a reasonable
agreement with the measurement made by Saldanha et al. [58]. However, the
production rate of 37Ar by neutrons is different from Saldanha et al. by a
factor of $\sim$3\. The discrepancy is caused by the different 40Ar(n,4n)37Ar
production cross section used in the process. This interaction cross section
has not been experimentally measured.

The total production rate of 39Ar measured by Saldanha et al. [58]
(934.6$\pm$156.2 atoms/kg$\cdot$day) is in a good agreement with this work
(938.53 atoms/kg$\cdot$day) using the Geant4 simulation. To validate the
simulated production rate of 42Ar from this work, we compared the cross
section of ${}^{40}\rm{Ar}(\alpha,2p)^{42}\rm{Ar}$ from a measurement made by
Yuki et al. [59], the ALICE code [60], a nuclear reaction simulator - TALYS
[61], and a nuclear database - ENDF [62] as shown in Figure 3. The cross
section extracted from Geant4 increases as the kinetic energy of $\alpha$
particles increases. This tendency is similar to the experimental data [59]
and the ALICE code. However, the value of the cross section extracted from
Geant4 is different by a factor of 2 to 5 (depending on energy) from the
available data provided by Yuki et al. and a factor of $\sim$2 from the ALICE
code. Note that the cross section obtained from the TALYS code is
significantly different from Geant4, the ALICE code, and the available data.
The ENDF-6 database only provides the cross section for energy below 30 MeV.
Thus, we can conclude that the cross section of
${}^{40}\rm{Ar}(\alpha,2p)^{42}\rm{Ar}$ is not well understood and more
experimental data are needed.

Figure 3:  Shown is a comparison for the cross section of
${}^{40}\rm{Ar}(\alpha,2p)^{42}\rm{Ar}$ extracted from Geant4 with the
experimental measurement (Yuki et al. [59]), as well as that of ALICE [60],
TALYS [61] and ENDF-6 [62].

Nevertheless, the Geant4 simulation confirms the production of 42Ar on the
Earth’s surface through the secondary $\alpha$ particles from cosmic ray
muons, neutrons, and protons. Note that the calculated rates of 42Ar
production in the atmospheric argon at sea level cannot be directly compared
to the experimentally observed abundance of 42Ar. This is because 42Ar is
mainly produced at the top of the atmosphere where the primary cosmic ray
fluxes are significantly larger than terrestrial cosmic ray fluxes [63]. In
addition, the energy spectra of the primary cosmic ray fluxes at the top of
the atmosphere are much harder than that of the cosmic rays at sea level.
Therefore, the experimental observed abundance of 42Ar is significant larger
than the calculated rates of 42Ar production at sea level.

###  2.3 Tolerable Exposure Time Limit

One can use the production rates of 39Ar and 42Ar obtained from this work to
set the allowed exposure time for UAr on the surface in terms of production,
transportation, and storage. Since the UAr produced by DarkSide in Colorado is
depleted in 39Ar by a factor of 1400 when compared to atmospheric Ar, we can
also expect that the concentration of 42Ar is depleted by at least a factor of
1400. Therefore, the cosmogenic production of 39Ar and 42Ar in the UAr should
be controlled in a level of less than 10% of the existing concentration of
39Ar and 42Ar in the UAr produced by DarkSide. This level of additional
cosmogenic production in the UAr will not compromise its sensitivity for the
planned DarkSide-20k and LEGEND-1000.

Using a constraint of 10% additional cosmogenic production of 39Ar and 42Ar in
the UAr, we estimate the allowed exposure time using the formula below [64]:

$N(t)=\frac{R}{\lambda}[1-exp(-\lambda t)],$ (1)

where $N(t)$ is the additional number of atoms (10% of the existing 39Ar or
42Ar in the UAr) produced in the UAr through cosmogenic activation when the
UAr is on the surface during the production, transportation, and storage, $R$
is the production rate from the Geant4 simulation, $\lambda$ is the decay
constant, and $t$ is the allowed exposure time. Note the 10% of the existing
39Ar in the UAr is calculated using the atmospheric level of
8.1$\times$10-16g/g divided by 1400 and multiplied by 10%. Similarly, the 10%
of the existing 42Ar in the UAr is calculated using the atmospheric level of
9.2$\times$10-21g/g divided by 1400 (assuming an upper limit for a
conservative consideration) and multiplied by 10%.

Derived from Eq (1), the allowed exposure time can be expressed as:

$t=\frac{ln[1-N(t)\lambda/R]}{-\lambda}.$ (2)

Using Eq (2), the allowed exposure time for the UAr on the surface are 954
days in terms of the production of 39Ar and 1702 days in terms of the
production 42Ar, respectively. This means that the UAr can be on the surface
for more than two years, which is a quite reasonable time to produce enough
UAr before transporting to underground for storage. Note that the cosmic ray
neutron flux on the Earth’s surface varies with altitude [63], which will
affect the production of 39Ar and 42Ar accordingly. The exposure time limit
will change with respect to different altitudes.

One can also show the calculated exposure limits using Eq. 2, for both 42Ar
and 39Ar production at sea level, versus a target activity level since
different experiments will have their own target activity levels. Figure 4
displays the results.

Figure 4:  Exposure time limit versus depleted ${}^{42}Ar/^{39}Ar$
concentration. The blue dots represent the reduction factor of 1400 in terms
of its ${}^{42}Ar$ and ${}^{39}Ar$ concentration in atmospheric argon,
respectively.

##  3 Summary

The cosmogenic activations of Ar for the next generation rare event search
experiments at sea level have been simulated using GEANT4 package. Fast
neutrons, muons and protons at the earth surface are considered individually.
The activation rates of ${}^{37}Ar$, 39Ar and ${}^{42}Ar$ are compared with
the measurement and model predictions[58]. We found that the production rates
of long-lived isotopes through cosmogenic activation on the earth’s surface
are dictated by fast neutrons. Using the production rates of 39Ar and 42Ar
from this work, we set a limit of the exposure time for the UAr on the surface
during its production, transportation, and storage. The limit is 954 days if
the production of 39Ar is a main driver for a dark matter experiment such as
DarkSide-20k. On the other hand, the limit is 1702 days if the production of
42Ar is a main concern for a 0$\nu\beta\beta$ decay experiment such as
LEGEND-100 using UAr as a veto. We conclude that the UAr from Colorado
discovered by DarkSide can be exposed on the earth’s surface in terms of
production, transportation, and storage for more than 900 days without
compromising the sensitivity of planned DarkSide-20k and LEGEND-1000
experiments.

##  4 Acknowledgement

The authors would like to thank Dr. Christina Keller for a careful reading of
this manuscript. This work is supported by NSF OISE-1743790, PHYS-1902577,
OIA-1738695, DOE FG02-10ER46709, the Office of Research at the University of
South Dakota and a research center supported by the State of South Dakota.

## References

  * [1] G. R. Blumenthal, S. Faber, J. R. Primack and M. J. Rees, Nature 311, 517 (1984).
  * [2] M. Davis, G. Efstathiou, C. S. Frenk and S. D. M. White, Astrophys. J. 292, 371 (1985).
  * [3] C. L. Bennett, et al., ApJS., 208, 20B (2013).
  * [4] M. W. Goodman and E. Witten, Phys. Rev. D31, 3059 (1985).
  * [5] J. L. Feng, Ann. Rev. Astron. Astrophys. 48, 495 (2010).
  * [6] CDMS Collaboration, Z. Ahmed et al., Z. Ahmed et al., Science 327 (2010) 1619.
  * [7] CDEX and TEXONO Collaboration, S. Liu et al., Phys. Rev. D90 no. 3, (2014) 032003, arXiv:1403.5421.
  * [8] CoGeNT Collaboration, C. Aalseth et al., 1) Phys. Rev. Lett. 106 (2011) 131301, arXiv:1002.4703.
  * [9] CRESST Collaboration, G. Angloher et al., 1) Phys. Rev. D 100 (2019) 102002.
  * [10] COUPP Collaboration, E. Behnke et al., Phys. Rev. D86 (2012) 052001, arXiv:1204.3094.
  * [11] DAMIC Collaboration, J. Barreto et al., Phys. Lett. B711 (2012) 264, arXiv:1105.5191.
  * [12] DarkSide Collaboration, P. Agnes et al., Phys. Rev. D 98 (2018) 102006.
  * [13] DRIFT Collaboration, J. B. R. Battat et al., Phys. Dark Univ. 9 (2014) 1, arXiv:1410.7821.
  * [14] EDELWEISS Collaboration, E. Armengaud et al., Phys. Rev. D 99 (2019) 082003.
  * [15] KIMS Collaboration, S.-C. Kim, J. Phys. Conf. Ser. 384 (2012) 012020.
  * [16] LUX Collaboration, D. S. Akerib et al., Phys. Rev. Lett. 118 no. 2, (2017) 021303, arXiv:1608.07648.
  * [17] PandaX Collaboration, X. Cui et al., Phys. Rev. Lett. 119 (2017) 181302.
  * [18] PICO Collaboration, C. Amole et al., Phys. Rev. Lett. 114 no. 23, (2015) 231302, arXiv:1503.00008.
  * [19] SuperCDMS Collaboration, I. Alkhatib et al., Phys. Rev. Lett. 127 (2021) 061801.
  * [20] XENON Collaboration, E. Aprile et al., Astropart. Phys. 34 (2011) 679, arXiv:1001.2834.
  * [21] XENON100 Collaboration, E. Aprile et al., Science 349 (2015) 851, arXiv:1507.07747.
  * [22] XENON Collaboration, E. Aprile et al., Phys. Rev. Lett. 123 (2019) 241803.
  * [23] XMASS Collaboration, K. Abe et al., Phys. Lett. B 789 (2019) 45.
  * [24] ZEPLIN-III Collaboration, D. Y. Akimov et al., Phys. Lett. B692 (2010) 180, arXiv:1003.5626.
  * [25] DEAP Collaboration, R. Ajaj et al., Phys. Rev. D 100 (2019) 022004.
  * [26] DAMA/LIBRA Collaboration, R. Bernabei et al., Eur. Phys. J. C67 (2010) 39, arXiv:1002.1028.
  * [27] D. S. Akerib et al. (LZ Collaboration) Nucl.Instrum.Meth.A 953 (2020) 163047, e-Print: 1910.09124.
  * [28] E. Aprile et al. Journal of Cosmology and Astroparticle Physics, Volume 2020, November 2020.
  * [29] C. E. Aalseth et al. Eur. Phys. J. Plus 133 (2018) 131. arXiV: 1707.08145.
  * [30] L. E. Strigari, New J. Phys. 11 (2009) 105011, arXiv:0903.3630.
  * [31] A. Gutlein et al., Astropart. Phys. 34 (2010) 90, arXiv:1003.5530.
  * [32] Steve R. Elliott and Petr Vogel, Annu. Rev. Nucl. Part. Sci. 115 (2002).
  * [33] Steve R. Elliott and Jonathan Engel, J. Phys. G: Nucl. Part. Phys 30, R183 (2004).
  * [34] F.T. Avignone III, G.S. King III and Yu. G. Zdesenko, New Journal of Physcis, 7, 6 (2005).
  * [35] A.S. Barabash, Physics of Atomic Nuclei, 67, No. 3, 438 (2004).
  * [36] M. Agostini et al. (GERDA Collaboration), Phys. Rev. Lett 111, 122503 (2013).
  * [37] J. B. Albert et al. (The EXO-200 Collaboration), Nature 510, 229234 (2014).
  * [38] A. Gando et al. (KamLAND-Zen Collaboration), A. Gando et al., Phys. Rev. Lett. 117 (2016) 082503.
  * [39] SNO+ Collaboration, S. Andringa et al., Adv. High Energy Phys. 2016 (2016) 6194250, arXiv:1508.05759.
  * [40] W. Xu et al. (Majorana Collaboration) Journal of Physics: Conference Series 606, 012004 (2015).
  * [41] N. Abgrall et al. (LEGEND Collaboration), arXiv: 2107.11462.
  * [42] M. G. Boulay and A. Hime, Astropart. Phys. 25 (2006) 179.
  * [43] A. S. Barabash, V. N. Kornoukhov, V. E. Jants, NIM A 385 (1997) 530-534.
  * [44] A. S. Barabash, R. R. Saakyan, V. I. Umatov, Journal of Physics: Conference Series 718 (2016) 062004.
  * [45] D.-M. Mei et al., arXiv:0912.5368, PRC 81 (2010) 055802.
  * [46] R. Ajaj et al. (The DEEP Collaboration), Phys. Rev. D 100 (2019) 072009.
  * [47] Walter M. Bonivento on behalf of the DarkSide-20k Collaboration, Journal of Physics: Conference Series 1468 (2020) 012234.
  * [48] S. Agostinelli, et al., Nucl. Instr. and Meth. A506 (2003) 250-303; K. Amako et al., IEEE Transactions on Nuclear Science 53 (2006) 270-278.
  * [49] http://www.slac.stanford.edu/comp/physics/geant4/slac_physics_lists/shielding/physlistdoc.html.
  * [50] Thomas K. Gaisser, ”Cosmic Rays and Particle Physics”, Cambridge Univeristy Press 1990.
  * [51] M. Guan, et. al., http://escholarship.org/uc/item/6jm8g76d.
  * [52] P. Archard et al. (L3+C Collaboration), Phys. Lett. B598 15 (2004).
  * [53] K. Abe, et. al., Phys. Lett. B564 8 (2003).
  * [54] M. S. Gordon, et. al., IEEE Transactions on Nuclear Science 51, 3427 (2004).
  * [55] C. Hagmann, D. Lange, D. Wright, UCRL-TM-229452, http://nuclear.llnl.gov/simulation/doc_cry_v1.7/cry_physics.pdf.
  * [56] G. Brooke and A. Wolfendale, Proc. Phys. Soc. 83, 843 (1964).
  * [57] P. Agnes, et al. (DarkSide collaboration), Phys. Rev. D 93, 081101(R) (2016).
  * [58] R. Saldanha, et al., Phys. Rev. C 100, 024608 (2019).
  * [59] Yuki H. et. al., CYRIC annual report, 3004, 109 (2004).
  * [60] M. Blann, K. H. Vonach, Phys. Rev. C 28 (1983) 1475.
  * [61] Arjan Koning, ”Status of TENDL: TENDL-2013 and 2014”, May 2014, WPEC meeting, NEA Data Bank.
  * [62] A.J. Koning et. al., Nucl. Data Sheets 155, 1 (2019).
  * [63] J. F. Ziegler, IBM Journal of Research Development, Vol. 40 NO. 1 January 1996.
  * [64] Kenneth S. Krane, Introductory Nuclear Physics, p.169, Published by John Wiley $\&$ Sons, 1987.