#  Soft and hard scales of the transverse momentum distribution in the Color
String Percolation Model

J. R. Alvarado García  Facultad de Ciencias Físico Matemáticas, Benemérita
Universidad Autónoma de Puebla, Apartado Postal 165, 72000 Puebla, Puebla,
México  D. Rosales Herrera  Facultad de Ciencias Físico Matemáticas,
Benemérita Universidad Autónoma de Puebla, Apartado Postal 165, 72000 Puebla,
Puebla, México  P. Fierro  Instituto de Física Luis Rivera Terraza, Benemérita
Universidad Autónoma de Puebla, Apartado Postal 165, 72000 Puebla, Puebla,
México  J. E. Ramírez<EMAIL_ADDRESS>Centro de
Agroecología, Instituto de Ciencias, Benemérita Universidad Autónoma de
Puebla, Apartado Postal 165, 72000 Puebla, Puebla, México  A. Fernández Téllez
Facultad de Ciencias Físico Matemáticas, Benemérita Universidad Autónoma de
Puebla, Apartado Postal 165, 72000 Puebla, Puebla, México  C. Pajares
<EMAIL_ADDRESS>Departamento de Física de Partículas and Instituto
Galego de Física de Altas Enerxías, Universidad de Santiago de Compostela,
E-15782 Santiago de Compostela, España

###### Abstract

In color string models, the transverse momentum distribution (TMD) is obtained
through the convolution of the Schwinger mechanism with the string tension
fluctuations distribution. Considering a $q$-Gaussian distribution for these
fluctuations, the TMD becomes a hypergeometric confluent function that
adequately reproduces the characteristic scales at low and high $p_{T}$
values. In this approach, the hard scale of the TMD is a consequence of
considering a heavy-tailed distribution for the string tension fluctuations
whose width rises as $\sqrt{s}$, multiplicity or centrality increases. In this
paper, we introduce the complete information of the TMD in the color string
percolation model by means of the determination of the color suppression
factor, which now also depends on the parameters of the $q$-Gaussian. To this
end, we analyze the reported data on pp and AA collisions at different center
of mass energies, multiplicities, and centralities. In particular, for minimum
bias pp collisions, we found that the $q$-Gaussian parameters and the
effective temperature are monotonically increasing functions of the center of
mass energy. Similar results are found for AA collisions as a function of the
centrality at fixed $\sqrt{s}$. We summarize these results in a phase diagram
that indicates the $q$-Gaussian parameters region allowing the quark-gluon
plasma formation.

††preprint: APS/123-QED

##  I Introduction

In high energy physics experiments, two projectiles (usually hadrons or
nuclei) traveling at relativistic velocities collide. These projectiles look
like thin disks because of the Lorentz contractionChaudhuri (2014); Busza  _et
al._ (2018). Instants before their interaction, color flux tubes are projected
into the transverse plane from the partons of each hadron or nuclei in the
collision Lappi and McLerran (2006), which are depicted as small disks
carrying color fields, namely color strings, and the fundamental color
interaction is given by the string overlapping Armesto  _et al._ (1996);
Ferreiro and Pajares (2012).

The number of the color strings in the systems grows, in general, with
increasing the center of mass energy, but also with multiplicity for pp
collisions, and size of the colliding objects and centrality for AA collisions
Pajares (2005). In consequence, strings start to overlap forming clusters in a
similar way as in the two-dimensional continuum percolation theory occurs
Nardi and Satz (1998). In this context, the color string percolation model
(CSPM) is a suitable framework for analyzing the properties of the colliding
systems. Here, the multiparticle production is described in terms of the color
strings, which break and create new strings through neutral color objects,
producing latter the observed hadrons. When the string clusters are formed,
they behave as a single color source with a color field corresponding to the
vector sum of the individual ones. The main result due to the random oriented
color summation of overlapped strings is a reduction of the multiplicity and
an increase of the string tension, hence an increase of the mean transverse
momentum Braun  _et al._ (2015a); Bautista  _et al._ (2019).

It is worth mentioning that the CSPM can explain most of the experimental data
on pp, pA, and AA collisions, including but not limited to the azimuthal
distribution of the produced particles, as well as the temperature dependence
of the ratio between the shear and bulk viscosities over the entropy density
Bautista  _et al._ (2012); Braun  _et al._ (2013, 2015b); Dias de Deus  _et
al._ (2016); Andrés  _et al._ (2017); Sahoo  _et al._ (2019). Moreover, it is
possible to introduce a parameter like a temperature in the CSPM (and other
color string models) by assuming the Schwinger mechanism for particle
production. If the string tension fluctuates, the transverse momentum
distribution (TMD) can be estimated by the convolution of the Schwinger
mechanism with the fluctuations distribution of the string tension. In
particular, if these fluctuations are modeled by a normal distribution, the
Schwinger mechanism becomes an exponential decay, the well-known thermal
distribution, which constant decay is related to the inverse of the effective
temperature Bialas (1999); Dias de Deus and Pajares (2006). Even though, this
result has been extensively used, it only adequately describes the TMD at low
transverse momentum values.

Usually, the transverse momentum spectrum is fully described by a fit function
composed by an exponential decay term plus a decreasing power-like
contribution representing the soft non-perturbative and hard perturbative QCD
collisions, respectively Bylinkin and Rostovtsev (2014); Bylinkin  _et al._
(2014). A shortcoming with this approach is the impossibility to identify the
string tension fluctuations that originate this TMD after their convolution
with the Schwinger mechanism.

The importance of the relation between the hard and soft parts of the TMD have
been emphasized recently in connection with the possible fast thermalization
of the quark gluon plasma produced in pp and AA collisions. In the same way
that in conformal field theory the energy cut-off for the ultraviolet modes
sets the effective thermal behavior of the system Calabrese and Cardy (2016);
Berges  _et al._ (2018), a hard parton interaction in a high energy collision
produces a rapid quench of the entangled partonic initial state and thus the
corresponding effective temperature, inferred from the exponential shape of
the TMD. In addition, it was shown that the fluctuations of the hard scale
leads to the effective temperature Baker and Kharzeev (2018); Feal  _et al._
(2019, 2021).

On the other hand, in a recent paper Pajares and Ramírez , C. Pajares and J.
E. Ramírez introduced a TMD obtained by considering that the string tension
fluctuations are a $q$-Gaussian distribution instead of a normal distribution.
In this way, the full TMD is described by a confluent hypergeometric function,
which adequately reproduces the exponential decay at low $p_{T}$ and also has
a power-like behavior at high $p_{T}$, describing the soft and hard scales of
the TMD that only depend on the $q$-Gaussian parameters. Using this approach
it has been possible to fit the TMD of charged hadrons produced in pp
collisions at different center of mass energies. Moreover, the TMD of Higgs
bosons reconstructed from $H\to\gamma\gamma$ and $H\to 4l$ decays in pp
collisions at $\sqrt{s}$=13 TeV are also well described. Notably, this
confluent hypergeometric function has been derived to describe the TMD by
considering complex arguments for string tension fluctuations Feal  _et al._
(2021).

In this paper, we aim to introduce the soft and hard scales of the transverse
momentum distribution into the CSPM. To do this, we explore modifications on
the average of the transverse momentum squared and the color suppression
factor produced by the introduction of the $q$-Gaussian as the distribution
that models the color string fluctuations.

The plan of this paper is as follows. In Sec. II, we present the CSPM and its
fundamental phenomenology. In Sec. III, we discuss the procedure to obtain the
TMD from considering that the string tension fluctuations are modeled by a
$q$-Gaussian distribution and how the hard part of the TMD contributes on
$\langle p_{T}^{2}\rangle$ and the color suppression factor. In Sec. IV we
describe the analysis of the experimental data and we show our main results.
Finally, Sec. V contains our conclusions and perspectives.

##  II Color String Percolation Model

The nature of the color fields and their interactions in particle collisions
are conveniently described by percolation theory of fully penetrable disks, as
we mentioned before. The color strings are represented in the transverse plane
by disks of radius $r\approx$ 0.2-0.3 fm Amelin  _et al._ (1993); Braun and
Pajares (2000); Braun  _et al._ (2015a); Bautista  _et al._ (2019). For
simplicity, all the color strings are considered of the same area $S_{1}$.
Each individual color source has a color charge arbitrarily oriented
$\mathbf{Q}_{1}$ with color field intensity $Q_{1}$ and color density
$\rho_{Q}=Q_{1}/S_{1}$, multiplicity $\mu_{1}$ and average squared transverse
momentum $\langle p_{T}^{2}\rangle_{1}$ Braun and Pajares (2000); Braun  _et
al._ (2015a); Bautista  _et al._ (2019). In this picture, parton interactions
are modeled by the overlapping of objects that carry color charge. Let us
discuss an example of a cluster formed of two overlapping strings as we depict
in Fig. 1. In this scenario, we have three different color sources: green (1),
cyan (2), and blue (3).

Figure 1: Sketch of the interaction between two color sources.

Notice that ${(1)}$ and ${(3)}$ are two independent color sources with equal
area $S^{(3)}=S^{(1)}$, such that $S^{(1)}=S_{1}-S^{(2)}$, while $S^{(2)}$ is
the area of the region (2), where the contributions of color fields sum one
another. We assume that each area individually produces charged particles.
This means that each area has its own color density, multiplicity, and
transverse momentum. Since the overlapping is partial, the color charge on
each region with no overlapping, color sources (1) and (3), is

$Q^{(1)}=\rho_{Q}S^{(1)}=Q_{1}(S^{(1)}/S_{1}).$ (1)

Notice that each color source contributes to the color charge on the
overlapping region (2) as

$\rho_{Q}S^{(2)}=Q_{1}(S^{(2)}/S_{1}).$ (2)

Thus the color charge $Q^{(2)}$ over $S^{(2)}$ is the vector sum considering
that color sources (1) and (3) are of equal color field intensity and in
average $\langle\mathbf{Q}_{1},\mathbf{Q^{\prime}}_{1}\rangle=0$ Biro  _et
al._ (1984); Ferreiro  _et al._ (2004); Ramírez  _et al._ (2021). Therefore

$Q^{(2)}=\sqrt{2}\rho_{Q}S^{(2)}=\sqrt{2}Q_{1}(S^{(2)}/S_{1}).$ (3)

In the same way, since $\mu^{(j)}$ is proportional to the color charge
$Q^{(j)}$ by $\mu^{(j)}=\mu_{1}Q^{(j)}/Q_{1}$ with $j=1,2$, the total
multiplicity can be expressed as Braun and Pajares (2000); Braun  _et al._
(2015a); Bautista  _et al._ (2019)

$\displaystyle\mu$ $\displaystyle=$ $\displaystyle 2\mu^{(1)}+\mu^{(2)}$ (4)
$\displaystyle=$ $\displaystyle
2\mu_{1}\frac{S^{(1)}}{S_{1}}+\sqrt{2}\mu_{1}\frac{S^{(2)}}{S_{1}}.$

Similarly, the average squared transverse momentum is Braun and Pajares
(2000); Braun  _et al._ (2015a); Bautista  _et al._ (2019)

$\displaystyle\langle p_{T}^{2}\rangle$ $\displaystyle=$ $\displaystyle
2\langle p_{T}^{2}\rangle^{(1)}+\langle p_{T}^{2}\rangle^{(2)}$ (5)
$\displaystyle=$ $\displaystyle
2\frac{\mu_{1}}{\mu}\frac{S^{(1)}}{S_{1}}\langle
p_{T}^{2}\rangle_{1}+\sqrt{2}\frac{\mu_{1}}{\mu}\sqrt{2}\frac{S^{(2)}}{S_{1}}\langle
p_{T}^{2}\rangle_{1}.$

By using Eq. (4) and Eq. (5) and taking into account that
$S^{(1)}+S^{(2)}=S_{1}$ we compute Braun and Pajares (2000); Braun  _et al._
(2015a); Bautista  _et al._ (2019)

$\displaystyle\frac{\langle p_{T}^{2}\rangle}{\langle p_{T}^{2}\rangle_{1}}$
$\displaystyle=$
$\displaystyle\frac{2(S^{(1)}/S_{1})+2(S^{(2)}/S_{1})}{2(S^{(1)}/S_{1})+\sqrt{2}(S^{(2)}/S_{1})}$
(6) $\displaystyle=$
$\displaystyle\frac{2}{2(S^{(1)}/S_{1})+\sqrt{2}(S^{(2)}/S_{1})}.$

Generalizing for a cluster of $N$ strings, $\mu$ and $\langle
p_{T}^{2}\rangle$ are given by

$\displaystyle\mu$ $\displaystyle=$
$\displaystyle\sum^{M}_{i}\mu^{(i)}=\sum^{M}_{i}\frac{\sqrt{n_{i}}{S}^{(i)}}{S_{1}}\mu_{1},$
(7) $\displaystyle\langle p_{T}^{2}\rangle$ $\displaystyle=$
$\displaystyle\sum^{M}_{i}\langle
p_{T}^{2}\rangle^{(i)}=\frac{{\sum_{i}}^{M}\left(\frac{n_{i}{S}^{(i)}}{S_{1}}\right)}{{\sum_{i}}^{M}\left(\frac{\sqrt{n_{i}}{S}^{(i)}}{S_{1}}\right)}\langle
p^{2}_{T}\rangle_{1},$ (8)

with $M$ being the total number of generated color sources in the cluster
Braun and Pajares (2000); Braun  _et al._ (2015a). It counts the number of the
regions $S^{(i)}$ with $n_{i}$ overlapped strings. Note that in the case
$n_{i}=1$ corresponds to the remnant fraction surfaces with no overlapping
(see Fig. 1). We use relations (8) and $\sum_{i}^{M}n_{i}S^{(i)}=NS_{1}$ to
establish

$N=\frac{\mu}{\mu_{1}}\frac{\langle p_{T}^{2}\rangle}{\langle
p_{T}^{2}\rangle_{1}},$ (9)

which suggests a conservation law for the transverse momentum of the produced
color sources Braun and Pajares (2000); Braun  _et al._ (2015a); Bautista  _et
al._ (2019).

Notice that Eq. (7) can be split as a sum expressed in terms of the total area
$S^{\text{tot}}_{n}$ defined as the sum of all regions that has $n=n_{i}$
overlapped disks. Thus,

$\displaystyle\frac{\mu}{\mu_{1}}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{1}}{S_{1}}\sum_{i,n_{i}=1}{S}^{(i)}_{n_{i}=1}+\frac{\sqrt{2}}{S_{1}}\sum_{i,n_{i}=2}{S}^{(i)}_{n_{i}=2}+\cdots+\frac{\sqrt{N}}{S_{1}}\sum_{i,n_{i}=N}{S}^{(i)}_{n_{i}=N}$
(10) $\displaystyle=$
$\displaystyle\sum_{n=1}^{N}\frac{\sqrt{n}}{S_{1}}\sum_{i,n_{i}=n}{S}^{(i)}_{n_{i}=n}=\sum_{n=1}^{N}\frac{\sqrt{n}}{S_{1}}S_{n}^{\text{tot}},$

and from Eq.(9) Braun and Pajares (2000); Braun  _et al._ (2015a); Bautista
_et al._ (2019)

$\frac{\langle p^{2}_{T}\rangle}{\langle
p^{2}_{T}\rangle_{1}}=N\frac{\mu_{1}}{\mu}=\frac{N}{\sum_{n=1}^{N}\frac{\sqrt{n}}{S_{1}}S_{n}^{\text{tot}}}.$
(11)

Assuming uniformly distributed strings, Eq. (10) becomes

$\frac{\mu}{\mu_{1}}=\frac{\langle\sqrt{n}\text{
}\rangle}{S_{1}}S=N\frac{\langle\sqrt{n}\text{ }\rangle}{\xi},$ (12)

where $\xi=NS_{1}/S$ is the string density or filling factor, which describes
the occupation of $N$ strings of area $S_{1}$ over the total interaction area
$S$ Braun and Pajares (2000); Braun  _et al._ (2015a); Bautista  _et al._
(2019). On the other hand, from Eq. (12) we can see a damping on multiplicity
given by

$F(\xi)=\frac{\mu}{N\mu_{1}}=\frac{\langle\sqrt{n}\text{ }\rangle}{\xi},$ (13)

where the average on the number of strings is taken from a Poisson
distribution and represents a reduction in multiplicity on a certain number of
strings Braun and Pajares (2000); Braun  _et al._ (2015a); Bautista  _et al._
(2019). We call this term the color suppression factor $F(\xi)$, which emerges
naturally from cluster formation. Notice that $F\leq 1$, because the number of
strings times the multiplicity of a single string is higher or equal to the
total multiplicity $\mu$. In the thermodynamic limit $F(\xi)$ is given by
Braun and Pajares (2000); Braun  _et al._ (2015a); Bautista  _et al._ (2019)

$\frac{1}{F^{2}(\xi)}=\frac{\xi}{1-e^{-\xi}}.$ (14)

Finally, from Eqs. (12) and (9) we can express multiplicity and squared
average transverse momentum as functions of the color suppression factor Dias
de Deus and Pajares (2006)

$\displaystyle\mu$ $\displaystyle=$ $\displaystyle\mu_{1}NF(\xi),$ (15)
$\displaystyle\langle p_{T}^{2}\rangle$ $\displaystyle=$ $\displaystyle\langle
p_{T}^{2}\rangle_{1}/F(\xi).$ (16)

In Sec. III we discuss the method to compute the color suppression factor by
using the information of the transverse momentum distribution. In a general
way, to estimate the filling factor $\xi$, it is necessary to solve the
equation $F(\xi)=b$, or equivalently

$\frac{1-\exp(-\xi)}{\xi}=b^{2},$ (17)

which solution is given by

$\xi=\frac{1}{b^{2}}+W\left(-\frac{\exp(-1/b^{2})}{b^{2}}\right),$ (18)

where $W$ is the Lambert function.

##  III Soft and hard scales of the TMD

In some color string models, the transverse momentum distribution of the
particles produced is described by the Schwinger mechanism Schwinger (1962);
Wong (1994):

$\frac{dN}{dp_{T}^{2}}\sim\exp\left(-\frac{\pi p_{T}^{2}}{x^{2}}\right),$ (19)

where $x$ is the string tension, which can fluctuate. If $P(x)$ is the
distribution describing the string tension fluctuations, then the TMD is
computed as the convolution of the Schwinger mechanism in Eq. (19) with the
string tension fluctuations as follows Bialas (1999)

$\frac{dN}{dp_{T}^{2}}\sim\int_{0}^{\infty}\exp\left(-\frac{\pi
p_{T}^{2}}{x^{2}}\right)P(x)dx.$ (20)

In particular, if $P(x)$ is a Gaussian distribution centered at 0 with
variance $\sigma^{2}=\langle x^{2}\rangle$, the Schwinger mechanism (19)
becomes a thermal distribution Bialas (1999); Dias de Deus and Pajares (2006)

$\frac{dN}{dp_{T}^{2}}\sim\exp\left(-\beta p_{T}\right),$ (21)

where $\beta=(2\pi/\sigma^{2})^{1/2}$ can be interpreted as the inverse
temperature of the system since Eq. (21) is similar to the Boltzmann
distribution Dias de Deus and Pajares (2006). The average of the transverse
momentum squared is

$\langle p_{T}^{2}\rangle=\frac{\int_{0}^{\infty}p_{T}^{2}\exp(-\beta
p_{T})dp_{T}}{\int_{0}^{\infty}\exp(-\beta
p_{T})dp_{T}}=\frac{\sigma^{2}}{\pi}.$ (22)

We can define a temperature for the color string percolation model by
comparing $\langle p_{T}^{2}\rangle$ computed from the Schwinger mechanism
(22) and the one deduced from the overlapping of color strings (16), which
implies $\sigma^{2}=\pi\langle p_{T}^{2}\rangle_{1}/F(\xi)$ Dias de Deus and
Pajares (2006). Thus

$T(\xi)=\sqrt{\frac{\langle p_{T}^{2}\rangle_{1}}{2F(\xi)}}.$ (23)

Even though Eq. (21) reproduces the soft part of the transverse momentum
distribution, it does not describe the hard part.

One natural and economical way (by introducing only one new parameter) to
describe both scales of the TMD is considering that the string tension
fluctuations are described by a $q$-Gaussian distribution Budini (2015)

$P(x)=\mathcal{N}\left(1+\frac{(q-1)x^{2}}{2\sigma^{2}}\right)^{\frac{1}{1-q}}$
(24)

centered at 0 with width of the distribution $\sqrt{\sigma^{2}}$. We also
assume that $1<q<3$. Here, $\mathcal{N}$ is the constant that guarantees the
normalization of the distribution on the domain of the string tension
fluctuation. Under these considerations and introducing the variable

$\tau=\frac{2\sigma^{2}}{(q-1)x^{2}},$ (25)

the TMD is computed as follows

$\frac{dN}{dp_{T}^{2}}\sim\frac{\Gamma\left(\frac{1}{q-1}\right)}{\sqrt{\pi}\Gamma\left(\frac{q-3}{2(1-q)}\right)}\int_{0}^{\infty}\exp\left(-\frac{\pi
p_{T}^{2}(q-1)}{2\sigma^{2}}\tau\right)\tau^{\frac{1}{q-1}-\frac{3}{2}}(1+\tau)^{\frac{1}{1-q}}d\tau.$
(26)

On the other hand, the confluent hypergeometric function is defined as Arfken
(1995)

$U(a,b,z)=\frac{1}{\Gamma(a)}\int_{0}^{\infty}\exp(-zt)t^{a-1}(1+t)^{b-a-1}dt.$
(27)

Comparing (26) and (27) we identify

$\displaystyle a=\frac{1}{q-1}-\frac{1}{2}$ and $\displaystyle b=\frac{1}{2}.$
(28)

Thus, we can write the TMD as

$\frac{dN}{dp_{T}^{2}}\sim\frac{1}{\sqrt{\pi}}\Gamma\left(\frac{1}{q-1}\right)U\left(\frac{1}{q-1}-\frac{1}{2},\frac{1}{2},\pi
p_{T}^{2}\frac{q-1}{2\sigma^{2}}\right),$ (29)

where $U$ is the confluent hypergeometric function, which has two well-known
asymptotic behaviors Slater (1972). At low $p_{T}$,

$\frac{dN}{dp_{T}^{2}}\sim\exp\left(-\frac{\sqrt{2\pi(q-1)}\Gamma\left(\frac{1}{q-1}\right)p_{T}}{\Gamma\left(\frac{1}{q-1}-\frac{1}{2}\right)\sigma}\right).$
(30)

Similarly to Eq. (21), we define the soft scale as

$T_{th}=\sigma\frac{\Gamma\left(\frac{1}{q-1}-\frac{1}{2}\right)}{\sqrt{2\pi(q-1)}\Gamma\left(\frac{1}{q-1}\right)}.$
(31)

At high $p_{T}$, the TMD (29) behaves as

$\frac{dN}{dp_{T}^{2}}\sim\frac{\Gamma\left(\frac{1}{q-1}\right)}{\sqrt{\pi}}\left(\frac{\pi
p_{T}^{2}(q-1)}{2\sigma^{2}}\right)^{\frac{1}{2}-\frac{1}{q-1}},$ (32)

and thus we define the hard scale as

$T_{H}=\sigma\sqrt{\frac{2}{\pi(q-1)}}\left(\frac{\sqrt{\pi}}{\Gamma\left(\frac{1}{q-1}\right)}\right)^{\frac{q-1}{q-3}}.$
(33)

Notice that both the soft and hard scales depend only on the parameters $q$
and $\sigma$ of the $q$-Gaussian distribution. Moreover, the ratio between the
soft and hard scales

$\frac{T_{H}}{T_{th}}=2\left(\frac{\sqrt{\pi}}{\Gamma\left(\frac{1}{q-1}\right)}\right)^{\frac{q-1}{q-3}}\frac{\Gamma\left(\frac{1}{q-1}\right)}{\Gamma\left(\frac{1}{q-1}-\frac{1}{2}\right)}$
(34)

only depends on $q$. In Fig. 2 we depict this ratio.

Figure 2: Ratio $T_{H}/T_{th}$ in Eq. (34) as a function of the parameter $q$.

Notice that $q-1$ marks the departure from the thermal behavior and determines
the power like behavior of the TMD at high $p_{T}$. When $q$ tends to 1 there
is not hard part. Later we will show that the normalized fluctuations of the
TMD depends only on $q-1$, increasing as $q-1$ does.

To connect this approach to the color string percolation model, we compare the
average of the transverse momentum squared computed in both cases. For the TMD
obtained by considering the string tension fluctuations as the $q$-Gaussian
distribution, we found

$\displaystyle\langle p_{T}^{2}\rangle$ $\displaystyle=$
$\displaystyle\frac{\int_{0}^{\infty}p_{T}^{2}U\left(\frac{1}{q-1}-\frac{1}{2},\frac{1}{2},\pi
p_{T}^{2}\frac{q-1}{2\sigma^{2}}\right)dp_{T}}{\int_{0}^{\infty}U\left(\frac{1}{q-1}-\frac{1}{2},\frac{1}{2},\pi
p_{T}^{2}\frac{q-1}{2\sigma^{2}}\right)dp_{T}}$ (35) $\displaystyle=$
$\displaystyle\frac{\sigma^{2}}{\pi(3-2q)},$

for $1<q<3/2$, otherwise $\langle p_{T}^{2}\rangle$ diverges. The upper bound
3/2 for $q$ has also been reported as the limit for the Bose-Einstein
condensation in nonextensive thermodynamics approaches Megías  _et al._
(2022). Comparing Eq. (35) with Eq. (16), the width of the $q$-Gaussian
distribution can be related to the color string percolation model through

$\sigma^{2}=\frac{\pi\langle p_{T}^{2}\rangle_{1}^{qG}(3-2q)}{F(\xi)}.$ (36)

Here, we have denoted $\langle p_{T}^{2}\rangle_{1}^{qG}$ with the superscript
$qG$ to emphasize that the mean of the transverse momentum squared of a single
string may depend on $q$. Notice that $\langle p_{T}^{2}\rangle$ rises as $q$
increases (see Eq. (35)), which can be done by increasing the string density
or enhancing $\langle p_{T}^{2}\rangle_{1}^{qG}$. The former case leads us to
an incompatible description of the CSPM with the experimental data. In what
follows, we focus our efforts on the description of the case where increasing
the value of $\langle p_{T}^{2}\rangle_{1}$ correctly reproduces $\langle
p_{T}^{2}\rangle$ in Eq. (35).

Figure 3: Ratio between $\sqrt{\langle p_{T}^{2}\rangle_{1}^{qG}}$ and
$\sqrt{\langle p_{T}^{2}\rangle_{1}}$ as a function of $q$.

We expect the analysis of the TMD at low $p_{T}$ give us a $T_{th}$ value that
matches the corresponding $T$ value for the CSPM. The condition $T_{th}=T$
implies that the color suppression factor is now computed as $F(\xi)=\langle
p_{T}^{2}\rangle_{1}/2T_{th}^{2}$, and also depends on the $q$-Gaussian
parameters. Moreover, we found that the modified average of the squared
transverse momentum of a single string is given by

$\sqrt{\langle p_{T}^{2}\rangle_{1}^{qG}}=\sqrt{\langle
p_{T}^{2}\rangle_{1}}\sqrt{\frac{q-1}{3-2q}}\frac{\Gamma\left(\frac{1}{q-1}\right)}{\Gamma\left(\frac{1}{q-1}-\frac{1}{2}\right)},$
(37)

where $\sqrt{\langle p_{T}^{2}\rangle_{1}}$ is estimated by direct comparison
with the value of the critical temperature computed in the lattice QCD
framework Dias de Deus and Pajares (2006) or the chemical freeze-out
temperature (experimentally determined) Scharenberg  _et al._ (2011), taking
values around 200 MeV. In Fig. 3 we show the behavior of the quotient
$\sqrt{\langle p_{T}^{2}\rangle_{1}^{qG}}/\sqrt{\langle
p_{T}^{2}\rangle_{1}}$. Notice this case means that the string density
coincides with those described by the CSPM. However, the average of the
transverse momentum squared of a single string should increase in order to
reproduce the hard scale of the TMD. One implication of this approach is the
modification of the average of the transverse momentum squared of the CSPM in
Eq. (16), which now reads

$\langle p_{T}^{2}\rangle=\frac{\langle
p_{T}^{2}\rangle_{1}^{qG}}{F(\xi)}=\frac{\langle
p_{T}^{2}\rangle_{1}}{F(\xi)}\frac{q-1}{3-2q}\left(\frac{\Gamma\left(\frac{1}{q-1}\right)}{\Gamma\left(\frac{1}{q-1}-\frac{1}{2}\right)}\right)^{2}.$
(38)

Clearly this new $\langle p_{T}^{2}\rangle$ is enhanced by the contribution of
the hard part of the TMD and it is larger than the predicted value by the CSPM
($\langle p_{T}^{2}\rangle_{1}/F(\xi)$).

Moreover, the color suppression factor is well determined by the parameters of
the $q$-Gaussian distribution, but we still require the value of $\langle
p_{T}^{2}\rangle_{1}$.

It is interesting to evaluate the variance normalized by the average squared
of the string tension fluctuations. To this end, we compute the moments
$\langle x^{n}\rangle=I_{n}/I_{0}$, with $I_{n}$ being the integral

$\displaystyle I_{n}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}x^{n}\left[1+\frac{q-1}{2\sigma^{2}}x^{2}\right]^{1/(1-q)}dx$
(39) $\displaystyle=$
$\displaystyle\frac{1}{2}\left(\frac{2\sigma^{2}}{q-1}\right)^{(n+1)/2}B\left(\frac{1}{q-1}-1-\frac{1}{2}(n-1),1+\frac{1}{2}(n-1)\right),$

where $B(k_{1},k_{2})=\Gamma(k_{1})\Gamma(k_{2})/\Gamma(k_{1}+k_{2})$ is the
beta function. The integrals $I_{n}$ are well defined if $q<(n+3)/(n+1)$. In
particular, the first two moments of the string tension fluctuations are

$\displaystyle\langle x\rangle$ $\displaystyle=$
$\displaystyle\sigma\sqrt{\frac{2}{\pi(q-1)}}\frac{\Gamma\left(\frac{1}{q-1}-1\right)}{\Gamma\left(\frac{1}{q-1}-\frac{1}{2}\right)},$
(40) $\displaystyle\langle x^{2}\rangle$ $\displaystyle=$
$\displaystyle\frac{2\sigma^{2}}{5-3q}.$ (41)

Thus, the variance normalized by the squared average is given by

$\frac{\langle x^{2}\rangle-\langle x\rangle^{2}}{\langle
x\rangle^{2}}=\frac{\pi(q-1)}{5-3q}\left(\frac{\Gamma\left(\frac{1}{q-1}-\frac{1}{2}\right)}{\Gamma\left(\frac{1}{q-1}-1\right)}\right)^{2}-1$
(42)

which only depends on $q$. This observable is a bounded monotonic increasing
function for $1<q<3/2$, taking values between $\pi/2-1$ and $\pi^{2}/4-1$. On
the other hand, given the value of the color suppression factor for a
particular process, the string density can be computed by using Eq. (18).
However, this estimation is adequate for large collisions systems like those
produced in heavy ion collisions but may fail to describe pp collisions. In
fact, recent studies have pointed out the relevance of finite size effects in
all the observables of the CSPM, including the color suppression factor, which
significantly deviates from the estimation of the thermodynamic limit Ramírez
_et al._ (2017); Ramírez and Pajares (2019); Texca García  _et al._ (2022).

##  IV Results

Figure 4: Samples of TMD (figures) and their fits (solid lines) for (a)
minimum bias pp collisions at different center of mass energies, (b) pp
collisions at 13 TeV with multiplicity classification, and (c) central heavy
ion collisions at different center of mass energies.

In this section we discuss our analysis on the experimental data of the
production of charged hadrons on pp and AA collisions under different
conditions of center of mass energy, centrality, and multiplicity. To this
end, we analyse the reported data on Refs. STAR Collaboration (2003); ALICE
Collaboration (2018); ALICE Collaboration (2013a); ALICE Collaboration (2019,
2016); ATLAS Collaboration (2016) with the ROOT 6 software. In all cases, the
confluent hypergeometric function in Eq. (29) is fitted to the TMD data. In
this way, we found the values of the $q$-Gaussian parameters $q$ and $\sigma$,
and then both the soft and hard scales are computed together with the color
suppression factor of the corresponding processes. In Fig. 4 we show a
selected sample of TMD and their corresponding fits. For the sake of notation,
we denote as $dN/dp_{T}^{2}$ the invariant yield of charged particles as a
function of $p_{T}$ normalized by the number of events, the pseudorapidity
interval, and $2\pi$ corresponding to the azimuthal angle. To improve
visualization, the TMD is scaled by a factor of $10^{L}$ in some instances. It
is worth mentioning that the approach discussed in this manuscript is also
valid for fitting the normalized differential primary charged particle cross
sections measured at 0.9, 2.76, and 7 TeV, which is proportional to the
invariant yield of produced charged particles by the measured luminosity ALICE
Collaboration (2013b) and only differ from the TMD by a multiplicative
constant. Moreover, some experiments report a TMD without the $p_{T}$
normalization. In those cases, the fitting function would be $p_{T}U$ instead
of $U$.

###  IV.1 Minimum bias pp collisions

In Fig. 4 (a), we show the minimum bias data of pp collisions at different
center of mass energies together with their fits according to Eq. (29).

For this case, we found that $q$ and the center of mass energy are related as
follows

$q(\sqrt{s}\text{ })=a_{q}\left(\ln\left(\sqrt{\frac{s}{s_{0}}}\text{
}\right)-\ln(0.2)\right)^{0.75}+c_{q}$ (43)

with $a_{q}$=0.0247(5), $c_{q}$=1.218(8), and $\sqrt{s_{0}}=$1 TeV. Figure 5
(a) shows the $q$ fitted data. Notice the agreement between the data and the
proposed function $q(\sqrt{s}\text{ })$.

Figure 5: (a) $q$ fitted values (squares) as a function of $\sqrt{s}$ and
their relation through Eq. (43) (dashed line). (b) $\sigma$ fitted values
(squares) as a function of $\sqrt{s}$ and their relation through Eq. (45)
(dashed line). (c) $T_{th}$ estimated (squares) by plugging the $q$ and
$\sigma$ fitted values on Eq. (31) as a function of $\sqrt{s}$ and their
relation through Eq. (44) (dashed line). (d) Ratio $T_{H}/T_{th}$ estimated
(squares) by plugging the $q$ and $\sigma$ fitted values on Eq. (34) as a
function of $\sqrt{s}$. Red solid lines on (c) and (d) are estimations of
$T_{th}$ and $T_{H}/T_{th}$ by using relations $q(\sqrt{s}\text{ })$ and
$\sigma(\sqrt{s}\text{ })$ on (43) and (45), respectively.

On the other hand, for $\sigma$ and $T_{th}$, we propose as fitting functions

$\displaystyle{T_{th}}(\sqrt{s}\text{ })$ $\displaystyle=$ $\displaystyle
a_{{T_{th}}}\left(\sqrt{\frac{s}{s_{0}}}\text{ }\right)^{c_{{T_{th}}}},$ (44)
$\displaystyle\sigma(\sqrt{s}\text{ })$ $\displaystyle=$ $\displaystyle
a_{\sigma}\left(\sqrt{\frac{s}{s_{0}}}\text{ }\right)^{c_{\sigma}}.$ (45)

We found $a_{T_{th}}=$0.1750(5) GeV and $c_{T_{th}}=$0.046(2), and
$a_{\sigma}=$0.392(1) GeV and $c_{\sigma}=$0.043(1). In Figs. 5 (b) and (c) we
show the values obtained for $\sigma$ and the value estimated of $T_{th}$,
respectively. In particular, Eq. (44) has a power law trend as the previous
relation reported in Ref. Bylinkin  _et al._ (2014) for $T_{th}$; however, we
estimated a slightly lesser value for the exponent. Moreover, in Fig. 5 (c) we
also depict the soft scale estimation by plugging into Eq. (31) the functions
$q(\sqrt{s}\text{ })$ and $\sigma(\sqrt{s}\text{ })$ from Eqs. (43) and (45),
finding an excellent match with the function (44).

Additionally, $T_{H}$ is estimated by plugging the fits for $q$ and $T_{th}$
on Eq. (34) and using the functions $q(\sqrt{s}\text{ })$ and
$\sigma(\sqrt{s}\text{ })$. These results are shown in Fig. 5 (d).

###  IV.2 Dependence on the multiplicity of pp collisions

Figure 6: Fitted values of (a) $q$ and (b) $\sigma$ for pp collisions at
$\sqrt{s}$= 5.02 (purple) and 13 TeV (pink) for the V0M multiplicity
classification, and (c) the estimation of the soft scale $T_{th}$, and (d) the
ratio $T_{H}/T_{th}$.

Other TMDs that we can study are those corresponding to the pp collisions with
multiplicity classification (see Fig. 4 (b)). In particular, we analyze the
data reported for pp collisions at $\sqrt{s}$=5.02 and 13 TeV by the VZERO
ALICE detector, which are classified by the V0M nomenclature. In Figs. 6 (a)
and (b) we show the behavior of $q$ and $\sigma$ as a function of
multiplicity, respectively. We observe some key features. Both $q$ and
$\sigma$ increase as multiplicity and center of mass energy does. Then, the
soft scale $T_{th}$ also increases with multiplicity for fixed $\sqrt{s}$ (see
Fig. 6 (c)). On the other hand, the hard scale $T_{H}$ can be computed by
multiplying $T_{th}$ by the quotient $T_{H}/T_{th}$ in Eq. (34), which only
depends on $q$ and it is a decreasing function.

###  IV.3 Heavy ions collisions

Figure 7: Fitted values of (a) $q$ and (b) $\sigma$ for AuAu collisions at
$\sqrt{s}$= 0.2 TeV (red), PbPb collisions at 2.76 (blue) and 5.02 (orange)
TeV, (c) estimation of $T_{th}$ and (d) ratio $T_{H}/T_{th}$, as functions of
centrality.

In principle, the TMD reported for heavy ion collision can also be described
by using the confluent hypergeometric function in Eq. (29). This approach gets
to describe the soft and hard parts of the TMD but deviates from the spectra
at mesoscale values of $p_{T}$ (see Fig. 4 (c)). It happens because the final
state interactions with a high density medium, the high $p_{T}$ particles are
suppressed and induce the jet quenching Pajares (2005). Then we fit lower
values of $q$ for the most central heavy ion collisions, as it is plotted in
Fig. 7 (a). We observe that those deviations are less important as the
centrality decreases and the momentum spectra look similar to the TMD of pp
collisions. Nevertheless, the soft scale $T_{th}$ increases as centrality and
center of mass energy does, as expected. Our results for $\sigma$ and $T_{th}$
have the same behavior as those for pp collisions. The dependence of $\sigma$
on the energy, the centrality, or the multiplicity is now codified into
$F(\xi)$ via Eq. (36). As energy or centrality (multiplicity) increases,
$F(\xi)$ decreases and $\sigma$ increases. Similar behavior is held for
$T_{th}$ because its value is mainly determined by $\sigma$, see Eq. (31). In
pp collisions as the energy or the multiplicity increases, there are more hard
collisions and therefore $q$ increases. This is not the case of heavy ion
collisions where is well known the suppression of the production of high
$p_{T}$ particles and jet quenching due to the interaction with the quark
gluon plasma produced, thus, $q$ decreases.

Above the percolation threshold of the string density a cluster of the
produced strings is formed covering most of the surface of the collision. The
diversity of string tensions would be small because there is only one large
cluster and a few others with small number of strings. In particular, we
observe that the normalized variance decreases with centrality. This is
clearly inferred from Fig. 7 (a). On the other hand, for the case of pp
collisions, the large number of hard collisions produces an increment on $q$,
and then an increment of the normalized variance. This different behavior
between pp and AA collisions can be traced back to the absence of jet
quenching in pp high multiplicity collisions, in contrast to the heavy ion
case.

###  IV.4 Phase diagram for QGP formation

As we stated in Sec. III, the computation of the color suppression factor is
now straightforward. For the value of $\sqrt{\langle p_{T}^{2}\rangle_{1}}$,
we consider the estimation using the chemical-freeze out temperature, taking
the value 0.207 GeV Scharenberg  _et al._ (2011). Thus, the color suppression
factor is calculated as

$F(\xi)=\frac{0.0214245\text{ GeV}^{2}}{T_{th}^{2}},$ (46)

where $T_{th}$ is the estimated soft scale from the $q$-Gaussian parameters.
In Fig. 8, we show our evaluations of $F(\xi)$. In all cases, the higher
values of $F(\xi)$, the lower string density the systems will have.
Particularly, the string density rises with the increasing center of mass
energy, multiplicity, or centrality, as expected.

Figure 8: Color suppression factor for the cases analyzed: (a) minimum bias pp
collisions, (b) pp collisions at $\sqrt{s}$= 5.02 and 13 TeV with multiplicity
classification, and (c) heavy ion collisions. Figures, lines, and colors are
the same as in Figs. 5, 6, and 7, respectively. Figure 9: Variance normalized
by the average square of the string tension fluctuations as a function of the
filling factor of the system for the TMD analyzed.

In Fig. 9, we show the normalized variance as a function of the string
density. For all cases, the variance increases as the filling factor grows. In
particular, the normalized variance for pp collisions has higher values for
the higher center of mass energies or multiplicities; meanwhile, it decreases
with growing centrality for heavy ion collisions.

Figure 10: Diagram for QGP formation in the $q$ \- $\sigma$ plane. Dotted and
dashed lines are the critical pairs $q$ \- $\sigma$ allowing the QGP formation
estimated in the thermodynamic limit and considering finite size effects for
pp collisions, respectively.

To close this section, let us extend the discussion on the implication of the
relation between the color suppression factor and the soft scale on Eq. (46).
In this way, $F(\xi)$ is also related to the $q$-Gaussian parameters as
mentioned before. Then, the analyzed LHC and RHIC data (at mid-rapidity
region) can be placed in the $q$ \- $\sigma$ plane. We can infer some relevant
implications. When no strings are overlapping, the color field is not
suppressed, and $F(\xi)=1$, which occurs for very low dense systems. Under
these conditions, the minimal soft temperature physically acceptable is
$T_{th,\text{min}}$=0.146 GeV. This leads to an exclusion region of forbidden
$q$ \- $\sigma$ pairs, depicted as a grey shaded region in Fig. 10. Another
noteworthy situation happens when the system reaches the percolation
threshold. At this point, the spanning cluster of color strings emerges, which
marks the onset of the QGP in high energy collisions. The best estimation for
the critical value of the color suppression factor in the thermodynamic limit
is $F_{c}=$0.77430816(4) Texca García  _et al._ (2022). Then, the $q$ \-
$\sigma$ pairs that give soft scales greater than the critical temperature
$T_{th,c}=$0.16634… GeV assures the formation of the QGP for heavy ion
collisions. The critical curve $T_{th,c}=T_{th}(q,\sigma)$ is represented as a
dotted line in Fig. 10. On the other hand, for small systems, the finite-size
effects on $F_{c}$ have been found to be $F_{c}-F_{cL}\propto L^{-1.3}$,
giving the value $F_{cL}\approx$0.7066 for pp collisions Texca García  _et
al._ (2022). The latter leads to a shift in the critical temperature for the
QGP formation. Under these conditions, $T_{cL}\approx 0.174...$ GeV for pp
collisions. Thus, small systems require more energetic collisions to form the
QGP than heavy ion collisions. The critical curve
$T_{th,c,L}=T_{th}(q,\sigma)$ for pp collisions is represented as a dashed
line in Fig. 10.

##  V Conclusions

We presented a natural extension of the color string percolation model to
include all the information of the transverse momentum distribution. To this
end, we have explored the possibility of having string tension fluctuations
distributed according to a $q$-Gaussian distribution instead of a normal one.
After the convolution of the Schwinger mechanism with the string tension
fluctuations we found that the TMD becomes a confluent hypergeometric
function, which can adequately fit the TMD experimental data reported for the
production of charged hadrons in pp and heavy ion collisions under a wide
range of conditions of center of mass energy, multiplicity classification or
(when applicable) centrality. It is worth mentioning that the confluent
hypergeometric function adequately reproduces the expected characteristics of
the TMD: exponential decay and power-like behaviors at low and high $p_{T}$
values, respectively. Then, the soft and hard scales of the TMD have natural
definitions.

To connect with the color string percolation model, we compared the prediction
of $\langle p_{T}^{2}\rangle$ by using both approaches. We observe that the
incorporation of the hard scale of the TMD rises the value of $\langle
p_{T}^{2}\rangle$, which can be possible by increasing the value of $\langle
p_{T}^{2}\rangle_{1}$. It means that the collision system has the same string
density as these modeled by the CSPM, but $\langle p_{T}^{2}\rangle_{1}$ grows
as $q$ takes higher values. In particular, notice that the ratio
$\sqrt{\langle p_{T}^{2}\rangle_{1}^{qG}}/\sqrt{\langle p_{T}^{2}\rangle_{1}}$
takes values between 1.19 and 1.4 for the typical values of $q$ and $\sigma$
obtained by fitting the TMD of the processes discussed above.

As we pointed out in Sec. IV.4, now the color suppression factor is
straightforwardly computed by its relation with the soft scale, given by Eq.
(46), which also depends on the parameters of the $q$-Gaussian. In this way,
we have introduced the soft and hard scales of the TDM into the color string
percolation model.

The introduction of the $q$-Gaussian fluctuations allow us to relate directly
the power like hard TMD with its normalized fluctuations, both depend only on
$q$. This fact is in line with the possibility that a hard collision produces
a rapid quench in the entangled initial partonic state, giving rise to an
exponential behavior at low $p_{T}$ and subsequently to the thermal
temperature.

We condensed all the values obtained for $q$ and $\sigma$ from fits in the $q$
\- $\sigma$ plane. By analyzing the characteristic values of the color
suppression factor, we found that the low-density limit produces an excluded
region of forbidden $q$ \- $\sigma$ pairs (shaded region in Fig. 10). Another
notable value is the critical color suppression factor. The picture of the
CSPM in the thermodynamic limit is adequate to describe heavy ion collision.
In the percolation threshold, the $F_{c}$ dictates the values of the $q$ \-
$\sigma$ pairs that mark the departure for QGP formation (see Fig. 10).
Particularly, our results for AuAu collisions at RHIC energies agree with the
claim of the QGP observed in those experiments. On the other hand, for small
systems, it is necessary to take into account the finite-size effects, leading
to an increment of the critical soft temperature for pp collisions which is
consistent with previous estimations. Notice that minimum bias pp collisions
require center of mass energies above 2 TeV to expect the QGP formation. On
the contrary, the QGP formation is not expected for pp collisions with a small
production of charged particles, as we showed in Fig. 10.

This work can be extended in several ways. For instance, it is possible to
analyze the TMD of pp collisions in producing a particular charged hadron.
This result could give insight into the hardness of the collisions required
for producing such charged hadron. It will be of interest if different $q$
values are observed for different charged hadrons. Moreover, since the $q$ and
$\sigma$ parameters are involved in determining the color suppression factor,
modifications on all observables of the CSPM are expected. Finally, we must
emphasize that the results presented here are also applicable to study TMDs
far from mid-rapidity region, which corresponds to non vanishing baryon-
chemical potential.

###### Acknowledgements.

C. P. thanks the grant Maria de Maeztu Unit of Excellence under the project
MDM-2016 0682 of Ministry of Science and Innovation of Spain. This work has
been funded by the projects PID2020-119632GB-100 of the Spanish Research
Agency, Centro Singular de Galicia 2019-2022 of Xunta de Galicia and the ERDF
of the European Union. This work was funded by Consejo Nacional de Ciencia y
Tecnología (CONACyT-México) under the project CF-2019/2042, graduated
fellowships grant numbers 645654, 1140160, and 848955, and postdoctoral
fellowship grant number 289198.

## References

  * Chaudhuri (2014) A. K. Chaudhuri,  _A short course on Relativistic Heavy Ion Collisions_ (IOPP, 2014) arXiv:1207.7028 [nucl-th] .
  * Busza  _et al._ (2018) W. Busza, K. Rajagopal, and W. van der Schee, Ann. Rev. Nucl. Part. Sci. 68, 339 (2018), arXiv:1802.04801 [hep-ph] .
  * Lappi and McLerran (2006) T. Lappi and L. McLerran, Nucl. Phys. A 772, 200 (2006), arXiv:hep-ph/0602189 .
  * Armesto  _et al._ (1996) N. Armesto, M. A. Braun, E. G. Ferreiro, and C. Pajares, Phys. Rev. Lett. 77, 3736 (1996), arXiv:hep-ph/9607239 .
  * Ferreiro and Pajares (2012) E. G. Ferreiro and C. Pajares, Phys. Rev. C 86, 034903 (2012), arXiv:1203.5936 [hep-ph] .
  * Pajares (2005) C. Pajares, Eur. Phys. J. C 43, 9 (2005), arXiv:hep-ph/0501125 .
  * Nardi and Satz (1998) M. Nardi and H. Satz, Phys. Lett. B 442, 14 (1998), arXiv:hep-ph/9805247 .
  * Braun  _et al._ (2015a) M. A. Braun, J. Dias de Deus, A. S. Hirsch, C. Pajares, R. P. Scharenberg, and B. K. Srivastava, Phys. Rept. 599, 1 (2015a), arXiv:1501.01524 [nucl-th] .
  * Bautista  _et al._ (2019) I. Bautista, C. Pajares, and J. E. Ramírez, Rev. Mex. Fis. 65, 197 (2019).
  * Bautista  _et al._ (2012) I. Bautista, J. Dias de Deus, and C. Pajares, Eur. Phys. J. C 72, 2038 (2012), arXiv:1102.3837 [hep-ph] .
  * Braun  _et al._ (2013) M. Braun, C. Pajares, and V. Vechernin, Nucl. Phys. A 906, 14 (2013), arXiv:1204.5829 [hep-ph] .
  * Braun  _et al._ (2015b) M. A. Braun, C. Pajares, and V. V. Vechernin, Eur. Phys. J. A 51, 44 (2015b), arXiv:1407.4590 [hep-ph] .
  * Dias de Deus  _et al._ (2016) J. Dias de Deus, A. S. Hirsch, C. Pajares, R. P. Scharenberg, and B. K. Srivastava, Phys. Rev. C 93, 024915 (2016), arXiv:1602.03437 [nucl-th] .
  * Andrés  _et al._ (2017) C. Andrés, M. Braun, and C. Pajares, Eur. Phys. J. A 53, 41 (2017), arXiv:1609.03927 [hep-ph] .
  * Sahoo  _et al._ (2019) P. Sahoo, S. K. Tiwari, S. De, R. Sahoo, R. P. Scharenberg, and B. K. Srivastava, Mod. Phys. Lett. A 34, 1950034 (2019), arXiv:1708.06689 [hep-ph] .
  * Bialas (1999) A. Bialas, Phys. Lett. B 466, 301 (1999), arXiv:hep-ph/9909417 .
  * Dias de Deus and Pajares (2006) J. Dias de Deus and C. Pajares, Phys. Lett. B 642, 455 (2006), arXiv:hep-ph/0607101 .
  * Bylinkin and Rostovtsev (2014) A. A. Bylinkin and A. A. Rostovtsev, Nucl. Phys. B 888, 65 (2014), arXiv:1404.7302 [hep-ph] .
  * Bylinkin  _et al._ (2014) A. A. Bylinkin, D. E. Kharzeev, and A. A. Rostovtsev, Int. J. Mod. Phys. E 23, 1450083 (2014), arXiv:1407.4087 [hep-ph] .
  * Calabrese and Cardy (2016) P. Calabrese and J. Cardy, J. Stat. Mech. 1606, 064003 (2016), arXiv:1603.02889 [cond-mat.stat-mech] .
  * Berges  _et al._ (2018) J. Berges, S. Floerchinger, and R. Venugopalan, Phys. Lett. B 778, 442 (2018), arXiv:1707.05338 [hep-ph] .
  * Baker and Kharzeev (2018) O. K. Baker and D. E. Kharzeev, Phys. Rev. D 98, 054007 (2018), arXiv:1712.04558 [hep-ph] .
  * Feal  _et al._ (2019) X. Feal, C. Pajares, and R. A. Vazquez, Phys. Rev. C 99, 015205 (2019), arXiv:1805.12444 [hep-ph] .
  * Feal  _et al._ (2021) X. Feal, C. Pajares, and R. A. Vazquez, Phys. Rev. C 104, 044904 (2021), arXiv:2012.02894 [hep-ph] .
  * (25) C. Pajares and J. E. Ramírez, “On the relation between the soft and hard parts of the transverse momentum distribution,” arXiv:2212.08887 .
  * Amelin  _et al._ (1993) N. S. Amelin, M. A. Braun, and C. Pajares, Phys. Lett. B 306, 312 (1993).
  * Braun and Pajares (2000) M. A. Braun and C. Pajares, Eur. Phys. J. C 16, 349 (2000), arXiv:hep-ph/9907332 .
  * Biro  _et al._ (1984) T. S. Biro, H. B. Nielsen, and J. Knoll, Nucl. Phys. B 245, 449 (1984).
  * Ferreiro  _et al._ (2004) E. G. Ferreiro, F. del Moral, and C. Pajares, Phys. Rev. C 69, 034901 (2004), arXiv:hep-ph/0303137 .
  * Ramírez  _et al._ (2021) J. E. Ramírez, B. Díaz, and C. Pajares, Phys. Rev. D 103, 094029 (2021), arXiv:2012.07920 [hep-ph] .
  * Schwinger (1962) J. Schwinger, Phys. Rev. 128, 2425 (1962).
  * Wong (1994) C.-Y. Wong,  _Introduction to high-energy heavy-ion collisions_ (World scientific, 1994).
  * Budini (2015) A. A. Budini, Phys. Rev. E 91, 052113 (2015), arXiv:1501.01037 [cond-mat.stat-mech] .
  * Arfken (1995) G. B. Arfken,  _Mathematical Methods for Physicists_ (Academic Press, California, 1995) Chap. 13, pp. 802–803.
  * Slater (1972) L. J. Slater, in  _Handbook of Mathematical Functions With Formulas, Graphs, and Mathematical Tables_ , edited by M. Abramowitz and I. A. Stegun (National Bureau of Standards, Washington, DC, 1972) Chap. 13, pp. 503–536.
  * Megías  _et al._ (2022) E. Megías, V. Timóteo, A. Gammal, and A. Deppman, Physica A 585, 126440 (2022).
  * Scharenberg  _et al._ (2011) R. P. Scharenberg, B. K. Srivastava, and A. S. Hirsch, Eur. Phys. J. C 71, 1510 (2011), arXiv:1006.3260 [nucl-ex] .
  * Ramírez  _et al._ (2017) J. E. Ramírez, A. Fernández Téllez, and I. Bautista, Physica A 488, 8 (2017), arXiv:1707.06395 [cond-mat.stat-mech] .
  * Ramírez and Pajares (2019) J. E. Ramírez and C. Pajares, Phys. Rev. E 100, 022123 (2019), arXiv:2004.10067 [cond-mat.stat-mech] .
  * Texca García  _et al._ (2022) J. C. Texca García, D. Rosales Herrera, J. E. Ramírez, A. Fernández Téllez, and C. Pajares, Phys. Rev. D 106, L031503 (2022), arXiv:2208.03769 [hep-ph] .
  * STAR Collaboration (2003) STAR Collaboration, Phys. Rev. Lett. 91, 172302 (2003), arXiv:nucl-ex/0305015 .
  * ALICE Collaboration (2018) ALICE Collaboration, JHEP 11, 013 (2018), arXiv:1802.09145 [nucl-ex] .
  * ALICE Collaboration (2013a) ALICE Collaboration, Eur. Phys. J. C 73, 2662 (2013a), arXiv:1307.1093 [nucl-ex] .
  * ALICE Collaboration (2019) ALICE Collaboration, Eur. Phys. J. C 79, 857 (2019), arXiv:1905.07208 [nucl-ex] .
  * ALICE Collaboration (2016) ALICE Collaboration, Phys. Lett. B 753, 319 (2016), arXiv:1509.08734 [nucl-ex] .
  * ATLAS Collaboration (2016) ATLAS Collaboration, Phys. Lett. B 758, 67 (2016), arXiv:1602.01633 [hep-ex] .
  * ALICE Collaboration (2013b) ALICE Collaboration, Eur. Phys. J. C 73, 2456 (2013b), arXiv:1208.4968 [hep-ex] .


exampleExample remarkRemark SIR–models for householdsP. Dönges, T. Götz, T.
Krüger, e.a.

# SIR–Model for Households

Philipp Dönges Max Planck Institute for dynamics and self–organization,
Göttingen, Germany<EMAIL_ADDRESS>Thomas Götz Mathematical
Institute, University Koblenz, Germany<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>Tyll Krüger
Polytechnic University Wroclaw, Poland<EMAIL_ADDRESS>Karol
Niedzielewski University Warsaw, Poland<EMAIL_ADDRESS>Viola
Priesemann11footnotemark: 1 Moritz Schäfer22footnotemark: 2

###### Abstract

Households play an important role in disease dynamics. Many infections
happening there due to the close contact, while mitigation measures mainly
target the transmission between households. Therefore, one can see households
as boosting the transmission depending on household size. To study the effect
of household size and size distribution, we differentiated the within and
between household reproduction rate. There are basically no preventive
measures, and thus the close contacts can boost the spread. We explicitly
incorporated that typically only a fraction of all household members are
infected. Thus, viewing the infection of a household of a given size as a
splitting process generating a new, small fully infected sub–household and a
remaining still susceptible sub–household we derive a compartmental ODE–model
for the dynamics of the sub–households. In this setting, the basic
reproduction number as well as prevalence and the peak of an infection wave in
a population with given households size distribution can be computed
analytically. We compare numerical simulation results of this novel
household–ODE model with results from an agent–based model using data for
realistic household size distributions of different countries. We find good
agreement of both models showing the catalytic effect of large households on
the overall disease dynamics.

###### keywords:

COVID–19, Epidemiology, Disease dynamics, SIR–model, Social Structure
MSC2020:

92D30, 93-10

##  1 Introduction

The spread of an infectious diseases strongly depends on the interaction of
the considered individuals. Traditional SIR–type models assume a homogeneous
mixing of the population and typically neglect the increased transmission
within closed subcommunities like households or school–classes. However,
literature indicates, that household transmission plays an important role [4,
3, 9].

In case of the COVID–pandemic, several studies have quantified the secondary
attack rate within households, i.e. the probability household-members get
infected, given that one household member is infected [10, 12, 7, 3, 9, 11].
The secondary attack rate depends on the virus variant, immunity and
vaccination status, cultural differences and mitigation measures within a
household (like e.g. early quarantine). In December 2021, when both the Delta
and Omicron variants were spreading, it ranged between about 19 % (Delta
variant in Norway) [9] to 39 % (Omicron in Spain) [3]. Hence, interestingly,
the secondary attack rate within a household is far from 100 % despite the
close contacts.

There are several models reported that try to include the contribution of
in–household transmission to the overall disease dynamics. The Reed–Frost
model describing in–household infections as a Bernoulli–process has be used by
[5, 4]. An average model assuming households always get completely infected
and ignoring the temporal dynamics has been proposed in [2]. Their findings
for the effective reproduction number agree with our results, see (5).
Extensions of differential equation based SIR–models in case of a uniform
household distribution have be proposed by [6, 8]. However, in real
populations, the household distribution is far from uniform, see [14]. Ball
and co–authors provided in [1] an overview of challenges posed by integrating
household effects into epidemiological models, in particular in the context of
compartmental differential equations models.

In this paper we develop an extended ODE SIR–model for disease transmission
within and between households of different sizes. We will treat the infection
of a given household as a splitting process generating two new sub–households
or household fragments of smaller size - representing the susceptible ($S$) or
fully infected ($I$) members. Each of the susceptible household fragments can
get infected and split later in time, while infected household segments
recover with a certain rate and are then immune (recovered or removed, $R$).
The dynamics of the susceptible, infected and recovered household segments is
modeled in Section 2. In Section 3, we compute the basic reproduction number
for our household model combining the attack rate inside a single households
and the transmission rate between individual households. The prevalence, as
the limit of the recovered part of the population can be computed analytically
— at least in case of small maximal household size, see Section 4. The peak of
a single epidemic wave in a population with given household distribution is
considered in Section 5. Again, for small maximal household size, we will be
able to compute analytically the maximal number of infected. Numerical
simulations based on realistic household size distributions for different
countries and a comparison with an agent–based model demonstrate the
applicability of the presented model. As outlook with focus on
non–pharmaceutical interventions we consider an extension of our model
quarantining infected sub–households based on a certain detection rate for a
single infected.

##  2 Household Model

Instead of real households we consider effective sub–households or household
fragments, still called households throughout this paper, being either fully
susceptible, infected or recovered. Let $S_{j},I_{j}$ and $R_{j}$ denote the
number of the respective households consisting of $j$ persons, where $1\leq
j\leq K$ and $K$ denoting the maximal household size. Then
$H_{j}=S_{j}+I_{j}+R_{j}$ equals to the total number of current sub–households
of size $j$. Furthermore, we introduce $H=\sum_{j=1}^{K}H_{j}$ as the total
number of households and $h_{j}=H_{j}/H$. The total population is given by
$N=\sum_{j=1}^{K}jH_{j}$. For further reference we also introduce the first
two moments of the household size distribution

$\displaystyle\mu_{1}:=\sum_{j=1}^{K}jh_{j}\quad\text{and}\quad\mu_{2}:=\sum_{j=1}^{K}j^{2}h_{j}\;.$

If an infection is brought into a susceptible sub–household of size $j$,
secondary infections will occur inside the household. We assume that each of
the remaining $j-1$ household members can get infected with probability $a$,
called the _in–household attack rate_. On average, we expect

$E_{j}=a(j-1)+1$

infections (including the primary one) inside a household of size $j$.
Existing field studies, see [12, 7] indicate an in–household attack rate in
the range of $16\%$–$30\%$ depending on the overall epidemiological situation,
household size and vaccinations. For the sake of tractability and simplicity,
our model assumes a constant attack rate $a$ independent of the household
size.

Let $b_{j,k}$ denote the probability, that a primary infection in a household
of size $j$ generates in total $k$ infections inside this household, where
$1\leq k\leq j$. The secondary infections give rise to a splitting of the
initial household of size $j$ into a new, fully infected sub–household of size
$k$ and another still susceptible sub–household of size $j-k$.

An infected household of size $k$ recovers with a rate $\gamma_{k}$ and
contributes to the overall force of infection by an out–household infection
rate $\beta_{k}$. We assume that the out–household reproduction number is
independent of the household size, i.e.

(1) $\mathcal{R}^{\ast}=\frac{\beta_{k}}{\gamma_{k}}=\text{constant
independent of $k$}\;.$

Now, the dynamical system governing the dynamics of the susceptible, infected
and recovered households of size $k$ reads as

(2a) $\displaystyle S_{k}^{\prime}$
$\displaystyle=Y\left[-kS_{k}+\sum_{j=k+1}^{K}\\!\\!jS_{j}\cdot
b_{j,j-k}\right]\;,$ (2b) $\displaystyle I_{k}^{\prime}$
$\displaystyle=-\gamma_{k}I_{k}+Y\sum_{j=k}^{K}jS_{j}\cdot b_{j,k}\;,$ (2c)
$\displaystyle R_{k}^{\prime}$ $\displaystyle=\gamma_{k}I_{k}\;,$ where
$\displaystyle Y$ $\displaystyle:=\frac{1}{N}\sum_{k=1}^{K}\beta_{k}\cdot
kI_{k}$

denotes the total force of infection.

For the recovery rate of an infected household of size $k$ we can consider the
two extremal cases and an intermediate case.

  1. 1.

_Simultaneous_ infections: All members of the infected household get infected
at the same time and recover at the same time, hence the recovery rate
$\gamma_{k}=\gamma_{1}$ is independent of the household size. Assumption (1)
leads to a constant out–household infection rate $\beta_{k}=\beta$.

  2. 2.

_Sequential_ infections: All members of the household get infected one after
another and the total recovery time for the entire household equals to $k$
times the individual recovery time. Hence the recovery rate
$\gamma_{k}=\gamma_{1}/k$ and by (1) we get $\beta_{k}=\beta_{1}/k$.

  3. 3.

_Parallel_ infections: The recovery times $T_{i}$, $i=1,\dots,k$ for each of
the $k$ infected individuals are modeled as independent exponentially
distributed random variables. Hence the entire household is fully recovered at
time $\max(T_{1},\dots,T_{k})$. The recovery rate $\gamma_{k}$ equals to the
inverse of the expected recovery time, i.e.

$\gamma_{k}=\frac{1}{E[\max(T_{1},\dots,T_{k})]}=\frac{1}{E[T_{1}]\cdot\theta_{k}}=\frac{\gamma_{1}}{\theta_{k}}\;,$

where $\theta_{k}=1+\tfrac{1}{2}+\dots+\tfrac{1}{k}\sim\log k+g$ denotes the
$k$–th harmonic number and $g\approx 0.5772\dots$ denotes the
_Euler–Mascheroni_ constant.

All these cases can be subsumed to

$\gamma_{k}=\frac{\gamma_{1}}{\eta_{k}},\quad\beta_{k}=\frac{\beta_{1}}{\eta_{k}}\;,$

where $\eta_{k}$ models the details of the temporal dynamics inside an
infected household.

In Figure 1 we compare these three cases in the scenario of a population with
maximal household size $K=6$. Each household size represents $1/6$ of the
entire population. Initially, $1\permil$ of the population is infected. The
out–household reproduction number is assumed to be $\mathcal{R}^{\ast}=1.33$,
in–household attack rate equals $a=0.2$ and the recovery rate $\gamma_{1}$ for
an infected individual equals $0.1$. Shown are the incidences, i.e. the daily
new infections over time. The three cases differ in the timing and the height
of the peak of the infection; for the simultaneous infections ($\gamma_{k}$
constant) the disease spreads fastest and for the sequential infections
($\gamma_{k}=\gamma_{1}/k$) the spread is significantly delayed. The cases of
parallel infections, i.e. $\gamma_{k}\simeq\gamma_{1}/(\log k+g)$, lies
between theses two extremes.

Figure 1:  Simulation of one epidemic wave in case of the three different
scalings of the recovery rates for households of size $k$. Shown is the
incidence, i.e. daily new infections over time. The green curve corresponds
equal recovery rates for all households, i.e. $\gamma_{k}=\gamma_{1}$. The red
curve corresponds to $\gamma_{k}=\gamma_{1}/k$ (sequential infections). The
blue curve corresponds to recovery rates obtaind from the maximum of
exponential distributions, i.e. $\gamma_{k}\simeq\gamma_{1}/(\log k)$.

######  Example 2.1.

Modeling the infections inside the household by a Bernoulli–process with
in–household attack rate (infection probability) $a\in[0,1]$, the total number
of infected persons inside the household follows a binomial distribution

(3) $b_{j,k}:=\binom{j-1}{k-1}a^{k-1}(1-a)^{j-k}\;.$

For the binomial in–household infection (3) it holds, that $E_{j}=a(j-1)+1$.

######  Theorem 2.2.

In model (2) the total population $N=\sum_{k=1}^{K}k\cdot(S_{k}+I_{k}+R_{k})$
is conserved.

######  Proof 2.3.

We have

$N^{\prime}=Y\left[-\sum_{k=1}^{K}k^{2}S_{k}+\sum_{k=1}^{K}k\sum_{j=k+1}^{K}jS_{j}b_{j,j-k}+\sum_{k=1}^{K}k\sum_{j=k}^{K}jS_{j}b_{j,k}\right]$

We consider the last two summands separately and reverse the order of
summation. Hence

$\displaystyle\sum_{k=1}^{K}k\sum_{j=k}^{K}jS_{j}b_{j,k}$
$\displaystyle=\sum_{j=1}^{K}jS_{j}\sum_{k=1}^{j}kb_{j,k}=\sum_{j=1}^{K}jE_{j}S_{j}$
and analogously $\displaystyle\sum_{k=1}^{K}k\sum_{j=k+1}^{K}jS_{j}b_{j,j-k}$
$\displaystyle=\sum_{j=2}^{K}jS_{j}\sum_{k=1}^{j-1}kb_{j,j-k}$
$\displaystyle=\sum_{j=2}^{K}jS_{j}\sum_{l=1}^{j-1}(j-l)b_{j,l}=\sum_{j=2}^{K}jS_{j}\sum_{l=1}^{j}(j-l)b_{j,l}$
$\displaystyle=\sum_{j=2}^{K}jS_{j}\left(j\sum_{l=1}^{j}b_{j,l}-\sum_{l=1}^{j}lb_{j,l}\right)=\sum_{j=2}^{K}j^{2}S_{j}-jE_{j}S_{j}$
$\displaystyle=\sum_{j=1}^{K}j^{2}S_{j}-jE_{j}S_{j}-\left(S_{1}-E_{1}S_{1}\right)\;.$

Due to $E_{1}=1$, the last term vanishes. Summing all the contributions, we
finally get

$N^{\prime}=Y\left[\sum_{j=1}^{K}-(j^{2}S_{j})+jE_{j}S_{j}+(j^{2}S_{j}-jE_{j}S_{j})\right]=0\;.$

##  3 Basic Reproduction Number

To compute the basic reproduction number for the household–model (2), we
follow the next generation matrix approach by Watmough and van den Driessche,
see [15]. We split the state variable
$x=(S_{1},\dots,R_{K})\in\mathbb{R}^{3K}$ into the infected compartment
$\xi=(I_{1},\dots I_{K})\in\mathbb{R}^{K}$ and the remaining components
$\chi\in\mathbb{R}^{2K}$. The infected compartment satisfies the differential
equation

$\xi_{k}^{\prime}=\mathcal{F}_{k}(\xi,\chi)-\mathcal{V}_{k}(\xi,\chi)$

where $\mathcal{F}_{k}=Y\sum_{j=k}^{K}jS_{j}\,b_{j,k}$,
$\mathcal{V}_{k}=\gamma_{k}\xi_{k}$ and
$Y=\tfrac{1}{N}\sum_{j=1}^{K}\beta_{j}j\xi_{j}$. Now, the Jacobians at the
disease free equilibrium $(0,\chi^{\ast})$ are given by

$\displaystyle F_{kj}$
$\displaystyle=\frac{\operatorname{\partial}\mathcal{F}_{k}}{\operatorname{\partial}\xi_{j}}(0,\chi^{\ast})=\frac{1}{N}j\beta_{j}\sum_{m=k}^{K}mH_{m}\,b_{m,k}$
$\displaystyle V_{kk}$
$\displaystyle=\frac{\operatorname{\partial}\mathcal{V}_{k}}{\operatorname{\partial}\xi_{k}}(0,\chi^{\ast})=\gamma_{k}\quad\text{and
$V_{kj}=0$ for $k\neq j$}$

The next generation matrix $G\in\mathbb{R}^{K\times K}$ is given by

$\displaystyle G$
$\displaystyle=FV^{-1}=\frac{1}{N}\left(\frac{j\beta_{j}}{\gamma_{j}}\sum_{m=k}^{K}mH_{m}\,b_{m,k}\right)_{k,j}$
Using the relation (1) we obtain $\displaystyle G$
$\displaystyle=\mathcal{R}^{\ast}\left(j\sum_{m=k}^{K}\frac{mH_{m}}{N}\,b_{m,k}\right)_{k,j}=\mathcal{R}^{\ast}\begin{pmatrix}\sum_{m=1}^{K}\tfrac{mH_{m}}{N}\,b_{m,1}\\\\[2.15277pt]
\sum_{m=2}^{K}\tfrac{mH_{m}}{N}\,b_{m,2}\\\ \vdots\\\
\tfrac{KH_{K}}{N}\,b_{KK}\end{pmatrix}\cdot\begin{pmatrix}1,2,\dots,K\end{pmatrix}$

where the last term represents the next generation matrix as a dyadic product
$G=\mathcal{R}^{\ast}\cdot dc^{T}$.

The basic reproduction number $\mathcal{R}=\rho(G)$ is defined as the spectral
radius $\rho(G)$ of the next generation matrix $G\in\mathbb{R}^{K\times K}$.
As a dyadic product, $G$ has rank one and hence there is only one non–zero
eigenvalue.

######  Lemma 3.1.

Let $c,d\in\mathbb{R}^{n}$ be two vectors with $c^{T}d\neq 0$. The non–zero
eigenvalue of the dyadic product $A=dc^{T}\in\mathbb{R}^{n\times n}$ is given
by $c^{T}d=\operatorname{tr}(A)$.

######  Proof 3.2.

Let $v$ be an eigenvector of $A$ to the eigenvalue $\lambda$. Then $\lambda
v=Av=(dc^{T})v=(c^{T}v)d$, where $c^{T}v\in\mathbb{R}$. If $\lambda\neq 0$,
then $v=\frac{c^{T}v}{\lambda}d$. Set
$\mu=\frac{c^{T}v}{\lambda}\in\mathbb{R}$. Now, $Av=A(\mu
d)=\mu(dc^{T})d=(c^{T}d)\mu d=(c^{T}d)v$; hence $\lambda=c^{T}d$.

As an immediate consequence we obtain

######  Theorem 3.3.

The basic reproduction number of system (2) is given by

(4)
$\mathcal{R}=\mathcal{R}^{\ast}\sum_{i=1}^{K}i\sum_{m=i}^{K}\frac{mH_{m}}{N}b_{m,i}=\mathcal{R}^{\ast}\sum_{m=1}^{K}\frac{mH_{m}}{N}E_{m}\;.$

######  Corollary 3.4.

In the particular situation, when the expected number of infections inside a
household of size $m$ is given by $E_{m}=a(m-1)+1$, the basic reproduction
number equals

(5)
$\mathcal{R}=\mathcal{R}^{\ast}\left[1+a\left(\frac{\mu_{2}}{\mu_{1}}-1\right)\right]\;.$

Assuming, that the infection of any member of a household results in the
infection of the entire household, i.e. in–household attack rate $a=1$, our
result $\mathcal{R}=\frac{\mu_{2}}{\mu_{1}}\mathcal{R}^{\ast}$ agrees with the
result obtained by Becker and Dietz in [2, Sect. 3.2].

##  4 Computing the prevalence

Let $z=\frac{1}{N}\sum_{k=1}^{K}kR_{k}$ denote the fraction of recovered
individuals. Then $z$ satisfies the ODE

$\displaystyle z^{\prime}$ $\displaystyle=\frac{Y}{\mathcal{R}^{\ast}}\;.$

In the sequel we will derive an implicit equation for the prevalence
$\lim_{t\to\infty}z(t)$ in two special cases:

  1. 1.

for maximal household size $K=3$. The procedure used here allows for immediate
generalization but the resulting expression gets lengthy and provide only
minor insight into the result.

  2. 2.

for in–household attack rate $a=1$ and arbitrary household sizes. In this
setting the equations (2a) for the susceptible households decouple and allow
and complete computation of the equation for the prevalence.

We will start with the second case. Let us consider $a=1$, i.e. inside
households infections are for sure and $E_{k}=k$. Then the ODE for the
susceptible households reads as $S_{k}^{\prime}=-YkS_{k}$ and we can insert
the recovered $z$

$\displaystyle S_{k}^{\prime}$ $\displaystyle=-k\mathcal{R}^{\ast}S_{k}z\;.$
After integration with respect to $t$ from $0$ to $\infty$, we get
$\displaystyle\ln\frac{S_{k}(\infty)}{S_{k}(0)}$
$\displaystyle=k\mathcal{R}^{\ast}\left(z(0)-z(\infty)\right)\;.$

Assuming initially no recovered individuals, i.e. $z(0)=0$ and considering the
total population at the the end of time, i.e.
$N=Nz(\infty)+\sum_{k}kS_{k}(\infty)$, we arrive at the system

$\displaystyle N$
$\displaystyle=Nz(\infty)+\sum_{k=1}^{K}kS_{k}(0)e^{-k\mathcal{R}^{\ast}z(\infty)}\;.$
Scaling with $N$, i.e. introducing $s_{k,0}=S_{k}(0)/N$, we get (6)
$\displaystyle 1$
$\displaystyle=z+\sum_{k=1}^{K}k\,s_{k,0}\,e^{-k\mathcal{R}^{\ast}z}\;.$ In
the limit $\mathcal{R}^{\ast}z\ll 1$ we obtain the approximation
$\displaystyle z$
$\displaystyle\sim\frac{1-\sum_{k}ks_{k,0}}{1-\mathcal{R}^{\ast}\sum_{k}k^{2}s_{k,0}}\;.$

Figure 2:  Prevalence in case of maximal household size $K=2$ vs. reproduction
number. If double households dominate (green), the prevalence is larger than
in the case of mostly single households (red), since households speed up the
infection dynamics.

Figure 2 shows the numerical solution of the prevalence equation (6) in case
of $K=2$ for reproduction numbers $\mathcal{R}^{\ast}\in[0,3]$ and initially
$2\%$ of the entire population being infected. The three different graphs
correspond to different initial values for the susceptible households: both
single and double households contain $49\%$ of the population (blue), $89\%$
of susceptibles live in single households and only $9\%$ in double households
(red) or just $9\%$ in single households and $89\%$ in double households
(green).

In case of arbitrary in–household attack rates $a\in[0,1]$, we consider the
situation for $K=3$. Setting $z=Z/N$, the relevant equations read as

$\displaystyle S_{3}^{\prime}$
$\displaystyle=-3\mathcal{R}^{\ast}z^{\prime}\,S_{3}$ $\displaystyle
S_{2}^{\prime}$
$\displaystyle=\mathcal{R}^{\ast}z^{\prime}\left(-2S_{2}+3b_{3,1}S_{3}\right)$
$\displaystyle S_{1}^{\prime}$
$\displaystyle=\mathcal{R}^{\ast}z^{\prime}\left(-S_{1}+2b_{2,1}S_{2}+3b_{3,2}S_{3}\right)\;.$

Solving the equations successively starting with
$S_{3}(t)=S_{3}(0)e^{-3\mathcal{R}^{\ast}z(t)}$ and using variation of
constants, we arrive at

$\displaystyle S_{2}(t)$
$\displaystyle=S_{2}(0)e^{-2\mathcal{R}^{\ast}z(t)}+3b_{3,1}S_{3}(0)\left(1-e^{-\mathcal{R}^{\ast}z(t)}\right)e^{-2\mathcal{R}^{\ast}z(t)}$
$\displaystyle S_{1}(t)$
$\displaystyle=S_{1}(0)e^{-\mathcal{R}^{\ast}z(t)}+2b_{2,1}\left(S_{2}(0)+3b_{3,1}S_{3}(0)\right)\left(1-e^{-\mathcal{R}^{\ast}z(t)}\right)e^{-\mathcal{R}^{\ast}z(t)}$
$\displaystyle\qquad+\tfrac{3}{2}\left(b_{3,2}-2b_{2,1}b_{3,1}\right)S_{3}(0)\left(1-e^{-2\mathcal{R}^{\ast}z(t)}\right)e^{-\mathcal{R}^{\ast}z(t)}\;.$

For the prevalence $z=\lim_{t\to\infty}z(t)$ we obtain the implicit equation

(7)
$1=z+\left[s_{1,0}+2b_{2,1}s_{2,0}+3(b_{2,1}b_{3,1}+\tfrac{1}{2}b_{3,2})s_{3,0}\right]e^{-\mathcal{R}^{\ast}z}\\\
\qquad+2(1-b_{2,1})\left[s_{2,0}+3b_{3,1}s_{3,0}\right]e^{-2\mathcal{R}^{\ast}z}\\\
+3\left[1+(b_{2,1}-2)b_{3,1}-\tfrac{1}{2}b_{3,2}\right]s_{3,0}e^{-3\mathcal{R}^{\ast}z}$

or in short

$1=z+c_{1}e^{-\mathcal{R}^{\ast}z}+c_{2}e^{-2\mathcal{R}^{\ast}z}+c_{3}e^{-3\mathcal{R}^{\ast}z}\;,$

where the coefficients $c_{1},c_{2}$ and $c_{3}$ are the above, rather lengthy
expressions involving the initial conditions and the in–household infection
probabilities $b_{j,k}$. In case of $K=2$, i.e. setting $s_{3,0}=0$, we get

(8)
$1=z+\left(s_{1,0}+2b_{2,1}s_{2,0}\right)e^{-\mathcal{R}^{\ast}z}+2\left(1-b_{2,1}\right)s_{2,0}e^{-2\mathcal{R}^{\ast}z}\;.$

In case of arbitrary household size $K>3$, the resulting equation for the
prevalence will have the same structure.

For the binomial infection distribution with $a=1$ this reduces to

$\displaystyle 1$
$\displaystyle=z+s_{1,0}e^{-\mathcal{R}^{\ast}z}+2s_{2,0}e^{-2\mathcal{R}^{\ast}z}+3s_{3,0}e^{-3\mathcal{R}^{\ast}z}$
and in case of $a=0$ we arrive at $\displaystyle 1$
$\displaystyle=z+\left(s_{1,0}+2s_{2,0}+3s_{3,0}\right)e^{-\mathcal{R}^{\ast}z}\;.$

In case of small initial infections, i.e.
$s_{1,0}+2s_{2,0}+3s_{3,0}=c_{1}+c_{2}+c_{3}\approx 1$, Eqn. (7) allows for
the trivial, disease free solution $z=0$. However, above the threshold
$\mathcal{R}_{c}=\frac{1}{c_{1}+2c_{2}+3c_{3}}$ the non–trivial endemic
solution shows up. Expanding the exponential for $\mathcal{R}^{\ast}z\ll 1$,
we arrive at

$\displaystyle 1\approx
z+\left(c_{1}+c_{2}+c_{3}\right)-\mathcal{R}^{\ast}z\left(c_{1}+2c_{2}+3c_{3}\right)+\frac{{\mathcal{R}^{\ast}}^{2}z^{2}}{2}\left(c_{1}+4c_{2}+9c_{3}\right)$
and hence $\displaystyle
z\left(\mathcal{R}^{\ast}(c_{1}+2c_{2}+3c_{3})-1\right)\approx\frac{{\mathcal{R}^{\ast}}^{2}z^{2}}{2}\left(c_{1}+4c_{2}+9c_{3}\right)\;.$

Besides the trivial solution $z=0$, this approximation has the second solution

(9)
$z\simeq\frac{2\left(\mathcal{R}^{\ast}(c_{1}+2c_{2}+3c_{3})-1\right)}{{\mathcal{R}^{\ast}}^{2}\left(c_{1}+4c_{2}+9c_{3}\right)}\;.$

If $\mathcal{R}^{\ast}>\mathcal{R}_{c}=(c_{1}+2c_{2}+3c_{3})^{-1}$, this
second root is positive.

The following Figure 3(left) shows the numerical solution of the prevalence
equation (8) in case of $K=2$ for reproduction numbers
$\mathcal{R}^{\ast}\in[0,3]$ and initially $2\%$ of the entire population
being infected while both single and double households contain $49\%$ of the
susceptible population. The different curves show the results depending on the
in–household attack rate $a$. The graph on the right depict the case $K=3$ and
$s_{1,0}=s_{2,0}=s_{3,0}=\frac{1}{6}$ for three different in–household attack
rates $a$. The approximation (9) is shown by the dashed curves. For
$\mathcal{R}^{\ast}<\mathcal{R}_{c}$, the trivial disease free solution $z=0$
is the only solution.

Figure 3:  Visualization of the prevalence vs. out–household reproduction
number $\mathcal{R}^{\ast}$ and for different in–household attack rate $a$.
(Left:) Maximal household size $K=2$. (Right:) Maximal household size $K=3$.
Also shown are the asymptotic approximation (9) and the critical value
$\mathcal{R}_{c}$.

##  5 Computing the peak of the infection

Let $J=\sum_{k=1}^{K}kI_{k}$ denote the total number of infected. Then $J$
satisfies

$\displaystyle J^{\prime}$
$\displaystyle=\sum_{k=1}^{K}kI_{k}^{\prime}=Y\left[-\frac{N}{\mathcal{R}^{\ast}}+\sum_{k=1}^{K}k\sum_{j=k}^{K}jS_{j}b_{j,k}\right]=Y\left[-\frac{N}{\mathcal{R}^{\ast}}+\sum_{j=1}^{K}jS_{j}\sum_{k=1}^{j}kb_{j,k}\right]$
$\displaystyle=Y\left[-\frac{N}{\mathcal{R}^{\ast}}+\sum_{j=1}^{K}jS_{j}E_{j}\right]\;.$

For $K=2$ we have to consider the problem

$\displaystyle S_{1}^{\prime}$
$\displaystyle=Y\left[-S_{1}+2S_{2}b_{2,1}\right]$ $\displaystyle
S_{2}^{\prime}$ $\displaystyle=Y\left[-2S_{2}\right]$ $\displaystyle
J^{\prime}$
$\displaystyle=Y\left[-\frac{N}{\mathcal{R}^{\ast}}+S_{1}+2E_{2}S_{2}\right]\;.$

Note, that $E_{2}=b_{2,1}+2b_{2,2}$ and $b_{2,1}+b_{2,2}=1$, hence
$b_{2,1}=2-E_{2}$. For sake of shorter notation and easier interpretation, we
introduce the scaled compartments $s_{1}=S_{1}/N$, $s_{2}=S_{2}/N$ and
$j=J/N$. Writing $s_{1}$ as a function of $s_{2}$, we get

$\displaystyle\frac{ds_{1}}{ds_{2}}$
$\displaystyle=E_{2}-2+\frac{s_{1}}{2s_{2}}$ with the solution $\displaystyle
s_{1}$ $\displaystyle=c\sqrt{s_{2}}-2(2-E_{2})s_{2}\;,$

where $c=2(2-E_{2})\sqrt{s_{2,0}}+s_{1,0}/\sqrt{s_{2,0}}$. For $j$ we have the
equation

$\displaystyle\frac{dj}{ds_{2}}$
$\displaystyle=\frac{-1/\mathcal{R}^{\ast}+s_{1}+2E_{2}s_{2}}{-2s_{2}}=\frac{1}{2\mathcal{R}^{\ast}s_{2}}-\frac{ds_{1}}{ds_{2}}-E_{2}$
with the solution given by $\displaystyle j$
$\displaystyle=\frac{1}{2\mathcal{R}^{\ast}}\ln\frac{s_{2}}{s_{2,0}}-s_{1}(s_{2})-E_{2}(s_{2}-s_{2,0})\;.$

The maximum of $j$ is either attained in case of an under–critical epidemic at
initial time at the necessary condition $j^{\prime}=0$ has to hold. Inserting
$s_{1}$ as a function of $s_{2}$ we arrive at the quadratic equation

$\displaystyle 4(E_{2}-1)x^{2}+cx-\frac{1}{\mathcal{R}^{\ast}}=0$ for
$x=\sqrt{s_{2}}$ with the positive root $\displaystyle
x=\sqrt{s_{2}}=\frac{c}{8(E_{2}-1)}\left(\sqrt{1+\frac{16}{\mathcal{R}^{\ast}}\frac{E_{2}-1}{c^{2}}}-1\right)\;.$

This solution is only meaningful, if $s_{1}+2s_{2}\leq 1$, i.e.

(10) $\frac{1}{\mathcal{R}^{\ast}}\leq 1+2(E_{2}-1)x^{2}$

which is an implicit equation for a threshold value of $\mathcal{R}^{\ast}$.

To obtain the value of $j$ at the maximum, we plug this root into the above
solution for $j$ and arrive at

$\displaystyle j_{max}$
$\displaystyle=1-\frac{1}{2\mathcal{R}^{\ast}}\left(1+\ln\frac{s_{2,0}}{x^{2}}\right)-\frac{cx}{2}\;.$

Hence, the peak of infection $j_{max}$ is a function of the out–household
reproduction number $\mathcal{R}^{\ast}=\beta_{k}/\gamma_{k}$, the expected
infections in double households $E_{2}$ and the initial conditions $s_{2,0}$,
$s_{1,0}$. In particular, it is independent of the scaling of the recovery
periods $\gamma_{k}$, as can be seen also in the more complex setting
presented in Figure 1.

The following Figure 4 shows the peak of infection $j_{max}$ versus the
out–household reproduction number $\mathcal{R}$. Solutions are only plotted,
if the condition (10) is satisfied. The left figure shows the situation for
$E_{2}=1.5$ and different initial conditions. The right figure shows the
variation with respect to $E_{2}=1+a$ keeping the initial conditions fixed as
$s_{1,0}=0.5$ and $s_{2,0}=0.25$.

Figure 4:  Peak height $j_{\max}$ of a single epidemic wave vs. out–household
reproduction number $\mathcal{R}^{\ast}$. (Left:) Variation with respect to
different initial household size distributions. (Right:) Variation with
respect to different in–household attack rates.

##  6 Simulations and comparison with agent–based model

The analysis of our household model presented in the previous sections, shows
that larger households have a strong effect on the spread of the epidemic. To
illustrate this, we simulate a single epidemic wave in populations with close
to realistic household distribution. For comparison we choose the household
size distributions for Bangladesh (BGD), Germany (GER) and Poland (POL) as
published by the UN statistics division in 2011, see [14]. Using a sample of
almost $1\,000\,000$ individuals, we obtain the distribution shown in Table 1.

|  |  Number $H_{k}$ of households of size
---|---|---
| Total | $k=1$ | $k=2$ | $k=3$ | $k=4$ | $k=5$ | $k\geq 6$
Bangladesh | 999 995 | 7 366 | 24 351 | 44 022 | 55 989 | 42 037 |  53 960 ($k=7$)
Germany | 999 999 | 173 640 | 154 920 | 67 846 | 48 585 | 15 201 |  7 106 ($k=6$)
Poland | 999 996 | 85 032 | 91 220 | 71 221 | 57 605 | 26 156 |  22 523 ($k=7$)
Table 1:  Distribution of household sizes in Bangladesh, Germany and Poland
for the year 2011, see [14]. The sample is based on a scaled population of
$1\,000\,000$. Differences are due to round–off effects.

######  Remark 6.1.

Note, that in Table 1 we have redistributed for Poland and Bangladesh the
fraction of population living in households of size $6$ or bigger to household
of size equal $7$ to match the total population. How one treats the population
represented by the tail in the household distribution can make a substantial
difference in the simulation outcomes. To visualize that, we show in Figure 5
simulations for out–household reproduction number $\mathcal{R}^{\ast}=0.9$,
in–household attack rate $a=0.2$ and three different versions of how to treat
the tail in the Polish household size distribution. The blue curve represents
the population distribution given in Table 1. The green curve uses the
extended census data including households up to size $10+$. All households of
size $10+$ are treated as households of size _equal_ to $10$. The red curve
show a simplified treatment of the tail, where all households of size $6+$ are
treated as being of size _equal_ to $6$.

The simulation redistributing the households of size $6+$ to size $7$ (blue)
agrees within the simulation accuracy with the detailed distribution (green).
The simplified treatment (red) shows already a significant difference. This
difference will be even more pronounced for smaller reproduction numbers
$\mathcal{R}^{\ast}$ since the smaller households first get subcritical with
decreasing $\mathcal{R}^{\ast}$. This result clearly visualizes the need for
careful treatment of the tail in the household size distribution, since big
households have an overproportional impact on the dynamics. Unfortunately,
most publicly available data truncates the household size distribution at size
$6$. From a modeling and simulation point of view, there is a need for more
detailed data on the household size distribution including information for
households of size bigger than $6$.

Figure 5:  Simulation of one infection wave using different treatment of the
tail in the household size distribution. (Blue:) Population distribution given
in Table 1 where households of size $6+$ are redistributed to size $7$.
(Green:) Polish census data including household sizes up to $10+$. (Red:)
Households of size $6+$ treated as being of size equal to $6$.

As initial condition for our simulation we assume $100$ infected single
households. The recovery rates inside the households are described using the
model of parallel infections, i.e. $\gamma_{k}=\gamma_{1}/(1+1/2+\cdots+1/k)$.
For the out–household reproduction number $\mathcal{R}^{\ast}$ we consider the
range $0.6\leq\mathcal{R}^{\ast}\leq 2.3$ and the in–household attack is
chosen as $a=0.25$. The ODE–system (2) is solved using a standard
RK4(5)–method.

Figures 6 and 7 show a comparison of the three countries. The prevalance and
the maximum number of infected shown in Fig. 6 clearly visualize that large
households are drivers of the infection. For moderate out–household
reproduction number $\mathcal{R}^{\ast}\simeq 1$, Bangladesh, with an average
household size $\mu_{1}=4.4$, faces a relative prevalence being approx. $50\%$
higher than Germany with an average household size of $\mu_{1}=2.1$ persons.
Also the peak number of infected is almost twice as high in Bangladesh
compared to Germany for $\mathcal{R}^{\ast}\in[1,1.25]$. Fig. 7 shows the
simulation of a single infection wave for moderate out–household reproduction
number $\mathcal{R}^{\ast}=1.1$ and in-household attack rate $a=0.25$. The
graph illustrates the faster and more severe progression of the epidemics in
case of larger households. The peak of the wave occurs in Bangladesh $60$ days
earlier and affects almost double the number of individuals compared to
Germany.

Figure 6:  Simulation for countries with different household size
distributions according to Table 1. The $x$–axis shows the out–household
reproduction number $\mathcal{R}^{\ast}$. The in–household attack rate is
fixed as $a=0.25$. Left: Relative Prevalence $z$ after $T=2\,000$ days. Right:
Absolute height of peak number of infected $J_{\max}$ (in thousands).

Figure 7:  Simulation of one infection wave for countries with different
household size distributions according to Table 1. Out–household reproduction
number $\mathcal{R}^{\ast}=1.1$ and in–household attack rate $a=0.25$. Plotted
is the total number of infected $J$ versus time.

To validate our model, we compared it to a stochastic, microscopic agent–based
model developed at the Interdisciplinary Centre for Mathematical and
Computational Modelling at the University of Warsaw. Complete details of this
model are given in [13]. In this model agents have certain states
(susceptible, infected, recovered, hospitalized, deceased, etc.) and infection
events occur in certain context, i.e. in households, on the streets or
workplaces and several more. Besides the household–context, the street–context
is used to capture infection events outside households. Since the ODE–model
(2) is a variant of an SIR–model, the agent–based model also just uses the
SIR–states and ignores all other states. The agent–based model uses for each
infected individual a recovery time that is sampled from an exponential
distribution with mean $10$ days. Based on the household distribution for
Poland, Figure 8 provides a comparison of the computed prevalence vs. the
out–household reproduction number for our model 2 and the agent–based model
(blue squares). The solid lines show the results of the ODE–model for
different in–household attack rates. The relative prevalence plotted in the
figure is defined as the total number of recovered individuals for time
$t\to\infty$; here we use $T=2\,000$ days for practical reasons.

Figure 8:  Simulation of one infection wave, plotted is the relative
prevalence $z$ versus the out–household reproduction number
$\mathcal{R}^{\ast}$. Initially $100$ infected single households in a
population according to Table 1. The solid lines show the results of the ODE
model for different in–household attack rates. The results of the agent–based
model are shown by the blue squares (exponentially distributed recovery time).

Figure 9 shows the results for the peak of the infection wave. The graph on
the left shows the peak of the incidences, i.e. the maximum number of daily
new infections and the graph on the right shows the time, when this peak
occurs.

Figure 9:  Simulation of one infection wave, initially $100$ infected single
households in a population according to Table 1. Left: Peak number of daily
new infections versus the out–household reproduction number
$\mathcal{R}^{\ast}$. Right: Time, when the peak occurs vs.
$\mathcal{R}^{\ast}$. The solid lines show the results of the ODE model for
different in–household attack rates. The results of the agent–based model are
shown by the blue squares.

For all three presented criteria: prevalence, peak height and peak time, our
ODE–household model (2) matches quite well with the agent–based simulation.
Best agreement can be seen for in–household attack rate $a$ in the range
between $0.16$ and $0.20$.

##  7 Effect of household quarantine

In this section we will extend our basic model (2) to households put under
quarantine after an infection is detected. Therefore, we introduce the
additional compartments $Q_{k}$ of quarantined households of size $k$. The
extended household model reads as

(11a) $\displaystyle S_{k}^{\prime}$
$\displaystyle=Y\left[-kS_{k}+\sum_{j=k+1}^{K}\\!\\!jS_{j}\cdot
b_{j,j-k}\right]$ (11b) $\displaystyle I_{k}^{\prime}$
$\displaystyle=-\gamma_{k}I_{k}-q_{k}I_{k}+Y\sum_{j=k}^{K}jS_{j}\cdot b_{j,k}$
(11c) $\displaystyle Q_{k}^{\prime}$ $\displaystyle=q_{k}I_{k}$ (11d)
$\displaystyle R_{k}^{\prime}$ $\displaystyle=\gamma_{k}I_{k}$

Here $q_{k}>0$ denotes the detection and quarantine rate for a household of
size $k$.

Let $q>0$ be the probability of infected individual to get detected. The
probability that a household consisting of $k$ infected gets detected equals
$d_{k}(q)=1-(1-q)^{k}$. Given recovery and detection rates $\gamma_{k}$ and
$q_{k}$ for a household of size $k$, the probability, that the household gets
detected before recovery equals $q_{k}/(\gamma_{k}+q_{k})$, which equals
$d_{k}(q)$. Therefore

$\frac{q_{k}}{\gamma_{k}}=\frac{d_{k}(q)}{1-d_{k}(q)}\sim kq$

for $q\ll 1$.

To compute the basic reproduction number, we repeat the computations for the
next generation matrix. The vector $\xi$ of infected compartments remains
unaltered and the new quarantine compartments $Q_{1},\dots,Q_{K}$ are included
in the vector $\chi$. The gain term $\mathcal{F}_{k}$ for the infected
compartments remains unaltered as well, but the loss term reads as
$\mathcal{V}_{k}=(\gamma_{k}+q_{k})\xi_{k}$. Accordingly, its Jacobian
modifies to $V_{kk}=\gamma_{k}+q_{k}$ and finally the next generation matrix
$G$ reads as

$\displaystyle G$
$\displaystyle=\frac{\beta}{N}\left(\frac{j}{\gamma_{j}+q_{j}}\sum_{m=i}^{K}mH_{m}\,b_{m,i}\right)_{i,j}$
$\displaystyle=\mathcal{R}^{\ast}\begin{pmatrix}\sum_{m=1}^{K}\tfrac{mH_{m}}{N}\,b_{m,1}\\\\[2.15277pt]
\sum_{m=2}^{K}\tfrac{mH_{m}}{N}\,b_{m,2}\\\ \vdots\\\
\tfrac{KH_{K}}{N}\,b_{KK}\end{pmatrix}\cdot\begin{pmatrix}\frac{1}{1+q_{k}/\gamma_{k}},\frac{2}{1+q_{2}/\gamma_{2}},\dots,\frac{K}{1+q_{K}/\gamma_{K}}\end{pmatrix}$

Its non–zero eigenvalue equals

$\mathcal{R}_{q}:=\mathcal{R}^{\ast}\sum_{k=1}^{K}\frac{k}{1+q_{k}/\gamma_{k}}\sum_{m=k}^{K}\frac{mH_{m}}{N}b_{m,k}=\mathcal{R}^{\ast}\sum_{m=1}^{K}\frac{mH_{m}}{N}\sum_{k=1}^{m}\frac{kb_{m,k}}{1+q_{k}/\gamma_{k}}\;.$

Due to denominator $1+q_{k}/\gamma_{k}$ we cannot interpret the last sum as
the expected infected cases in a household of size $m$. In the asymptotic case
of small quarantine rates $q_{k}\ll 1$, we can use the series expansion
$(1+q_{k}/\gamma_{k})^{-1}\sim 1-q_{k}/\gamma_{k}\sim 1-kq$ and get

$\displaystyle\mathcal{R}_{q}\sim\mathcal{R}^{\ast}\sum_{m=1}^{K}\frac{mH_{m}}{N}\left(E_{m}-q\sum_{k=1}^{m}k^{2}b_{m,k}\right)\;.$

So, the second moment of the in–household infection distribution comes into
play.

######  Remark 7.1.

If the expectation $E_{m}$ and the the second moment
$\sum_{k=1}^{m}k^{2}b_{m,k}$ of the in–household infections scale linear with
the household size $m$, then the above reproduction number depends on also on
the third moment of the household size distribution. This leads to a
straight–forward extension of the result (5).

In case of a binomial distribution of the in–household infections, it holds
that $E_{m}=a(m-1)+1$ and
$\sum_{k=1}^{m}k^{2}b_{m,k}=a(m-1)\left(3+a(m-2)\right)+1$. For the
reproduction number we arrive at the expression

$\mathcal{R}_{q}\sim\mathcal{R}^{\ast}\left[1-\frac{q}{\gamma}+a\left(1-\frac{3q}{\gamma}(1-a)\right)\left(\frac{\mu_{2}}{\mu_{1}}-1\right)-\frac{q}{\gamma}a^{2}\left(\frac{\mu_{3}}{\mu_{1}}-1\right)\right]\;,$

where $\mu_{3}:=\sum_{k=1}^{K}k^{3}h_{k}$ denotes the third moment of the
household size distribution.

The effectiveness of quarantine measures can be assessed by relating the
reduction in prevalence to the quarantined individuals. For a given detection
probability $q>0$, let $Z(q)$ denote the prevalence observed under these
quarantine measures and let $Q(q)$ denote the total number of individuals
quarantined. We define the effectiveness of the quarantine measures as

$\eta_{q}:=\frac{\left|Z(q)-Z(0)\right|}{Q(q)}\;.$

The numerator equals to the reduction in prevalence and hence $\eta$ can be
interpreted as infected cases saved per quarantined case.

In Figure 10 we compare the effectiveness in the setting of household
distributions resembling Germany (predominantly small households) and
Bangladesh (large households are dominant) for out–household reproduction
numbers in the range between $0.9$ and $1.2$. The results indicate, that
quarantine measures are less effective in societies with larger households
since infection chains inside households are not affected by the quarantine.
On the other hand, quarantine measures seem most effective in presence of
small households and for low reproduction numbers.

Figure 10:  Effectiveness $\eta$ of quarantine measures for the different
household distributions in Germany (GER) and Bangladesh (BGD), see Table 1.

##  8 Conclusion and Outlook

The findings of the compartmental household model (2) are in rather good
agreement with microscopic agent–based simulations. Despite its several
simplifying assumptions, the prevalence and the peak of the infection wave are
reproduced quite well.

However, field data suggests that the in–household attack rate $a$ is
significantly higher for $2$–person households than for larger households, see
[7, 12]. Analyzing the data given in [7] the attack rate in two–person
households is by a factor of $1.25$–$1.5$ larger than in households with three
or more members. According to [12] this may be caused by spouse relationship
to the index case reflecting intimacy, sleeping in the same room and hence
longer or more direct exposure to the index case. An extension of the model
(3) to attack rates $a_{k}$ depending on the household size $k$ is
straightforward. Carrying over the analytical computation of the reproduction
number (5) would lead to a variant of (5) including weighted moments of the
household size distribution.

When considering quarantine for entire households, the model shows some
interesting outcomes. Defining the ratio of reduction in prevalence and number
of quarantined persons as the effectivity of a quarantine, the model suggests
that quarantine is more effective in populations where small households
dominate. Whether this finding is supported by field data is a topic for
future literature research.

Vaccinations are not yet included in the model. From the point of view of
public health, the question arises whether the limited resources of vaccines
should be distributed with preference to larger households. Given the
catalytic effect of large households on the disease dynamics, one could
suggest to vaccinate large households first. However, a detailed analysis of
this question will be subject of follow–up research.

When extending our ODE SIR–model to an SIRS–model including possible
reinfection a recombination of the recovered sub–households is required.
Again, this extension will be left for future work. A relaxation of the
distribution of the recovered sub–household distribution to the overall
household distribution seems a tractable approach to tackle this issue.

## References

  * [1] F. Ball, T. Britton, T. House, V. Isham, D. Mollison, L. Pellis, and G. S. Tomba, Seven challenges for metapopulation models of epidemics, including households models, Epidemics, 10 (2015), pp. 63–67.
  * [2] N. G. Becker and K. Dietz, The effect of household distribution on transmission and control of highly infectious diseases, Mathematical biosciences, 127 (1995), pp. 207–219.
  * [3] J. Del Águila-Mejía, R. Wallmann, J. Calvo-Montes, J. Rodríguez-Lozano, T. Valle-Madrazo, and A. Aginagalde-Llorente, Secondary attack rate, transmission and incubation periods, and serial interval of sars-cov-2 omicron variant, spain, Emerging Infectious Diseases, 28 (2022), p. 1224.
  * [4] C. Fraser, D. A. Cummings, D. Klinkenberg, D. S. Burke, and N. M. Ferguson, Influenza transmission in households during the 1918 pandemic, American journal of epidemiology, 174 (2011), pp. 505–514.
  * [5] K. Glass, J. McCaw, and J. McVernon, Incorporating population dynamics into household models of infectious disease transmission, Epidemics, 3 (2011), pp. 152–158.
  * [6] T. House and M. J. Keeling, Deterministic epidemic models with explicit household structure, Mathematical biosciences, 213 (2008), pp. 29–39.
  * [7] T. House, H. Riley, L. Pellis, K. B. Pouwels, S. Bacon, A. Eidukas, K. Jahanshahi, R. M. Eggo, and A. Sarah Walker, Inferring risks of coronavirus transmission from community household data, Statistical Methods in Medical Research, 31 (2022), pp. 1738–1756.
  * [8] G. Huber, M. Kamb, K. Kawagoe, L. M. Li, B. Veytsman, D. Yllanes, and D. Zigmond, A minimal model for household effects in epidemics, Physical Biology, 17 (2020), p. 065010.
  * [9] S. B. Jørgensen, K. Nygård, O. Kacelnik, and K. Telle, Secondary attack rates for omicron and delta variants of sars-cov-2 in norwegian households, JAMA, 327 (2022), pp. 1610–1611.
  * [10] W. Li, B. Zhang, J. Lu, S. Liu, Z. Chang, C. Peng, X. Liu, P. Zhang, Y. Ling, K. Tao, et al., Characteristics of household transmission of covid-19, Clinical Infectious Diseases, 71 (2020), pp. 1943–1946.
  * [11] F. P. Lyngse, L. H. Mortensen, M. J. Denwood, L. E. Christiansen, C. H. Møller, R. L. Skov, K. Spiess, A. Fomsgaard, R. Lassaunière, M. Rasmussen, et al., Household transmission of the sars-cov-2 omicron variant in denmark, Nature communications, 13 (2022), pp. 1–7.
  * [12] Z. J. Madewell, Y. Yang, I. M. Longini, M. E. Halloran, and N. E. Dean, Household transmission of sars-cov-2: a systematic review and meta-analysis, JAMA network open, 3 (2020), pp. e2031756–e2031756.
  * [13] K. Niedzielewski, J. M. Nowosielski, R. P. Bartczuk, F. Dreger, Ł. Górski, M. Gruziel-Słomka, A. Kaczorek, J. Kisielewski, B. Krupa, A. Moszyński, et al., The overview, design concepts and details protocol of icm epidemiological model (pdyn 1.5), (2022).
  * [14] UN statistics division, Demographic and social statistics.  https://unstats.un.org/unsd/demographic-social/products/dyb/documents/household/4.pdf, 2015\.  Accessed: 2023-01-09.
  * [15] P. Van den Driessche and J. Watmough, Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, Mathematical biosciences, 180 (2002), pp. 29–48.