India Based Neutrino Observatory,
Physics Reach and Status Report

D. Indumathi1,

1 The Institute of Mathematical Sciences, Chennai and Homi Bhabha National
Institute, Mumbai

<EMAIL_ADDRESS>
21st International Symposium on Very High Energy Cosmic Ray Interactions
(ISVHE- CRI 2022)
Online, 23-28 May 2022 10.21468/SciPostPhysProc.?

## Abstract

The India-based Neutrino Observatory (INO) is a proposed underground facility
located in India that will primarily house the magnetised Iron CALorimeter
(ICAL) detector to study atmospheric neutrinos produced by interactions of
cosmic rays with Earth’s atmosphere. The physics goal is to to make precision
measurements of the neutrino mixing and oscillation parameters through such a
study. We present here the results from detailed simulations studies, as well
as a status report on the project. In particular, we highlight the sensitivity
of ICAL to the open issue of the neutrino mass ordering, which can be
determined independent of the CP phase at ICAL.

###### Contents

  1. 1 Introduction
  2. 2 The INO project
    1. 2.1 The ICAL detector
  3. 3 Atmospheric Neutrinos and ICAL
    1. 3.1 ICAL Detector Response
    2. 3.2 Precision measurement of the 2–3 parameters
    3. 3.3 Other measurements at ICAL
  4. 4 Prototypes and the mini-ICAL detector
  5. 5 Conclusion

##  1 Introduction

The possibility of neutrino oscillations was originally suggested by
Pontecorvo (actually, he examined the possibility of neutrino anti-neutrino
oscillations) [1]. The theory of neutrino mixing in three-flavours was given
by Maki, Nakagawa and Sakata [2] and the effect of matter on the propagation
of neutrinos was discussed by Wolfenstein [3], Mikheev and Smirnov [4]. It is
now well-established that neutrino flavours mix, and hence neutrino
oscillations have been observed in diverse instances such as in solar
neutrinos, atmospheric neutrinos, reactor neutrinos, etc.; see Ref. [5] for a
historical review. Neutrino mixing in three flavours is parametrised by the
PMNS matrix, $U_{PMNS}$, given by

$\displaystyle
U_{PMNS}=\left(\begin{array}[]{ccc}c_{12}c_{13}&s_{12}c_{13}&s_{13}e^{-i\delta_{CP}}\\\
-c_{23}s_{12}-s_{23}s_{13}c_{12}e^{i\delta_{CP}}&c_{23}c_{12}-s_{23}s_{13}s_{12}e^{i\delta_{CP}}&s_{23}c_{13}\\\
s_{23}s_{12}-c_{23}s_{13}c_{12}e^{i\delta_{CP}}&-s_{23}c_{12}-c_{23}s_{13}s_{12}e^{i\delta_{CP}}&c_{23}c_{13}\end{array}\right).$
(4)

Here $c_{12}=\cos\theta_{12}$, $s_{12}=\sin\theta_{12}$ etc., and
$\delta_{CP}$ denotes the CP violating (Dirac) phase. By definition, the 3
$\times$ 3 neutrino mass matrix $M_{\nu}$ is diagonalised in the charged-
lepton mass basis. The parameters involved are the mixing angles $\theta_{ij}$
and the mass-squared differences $\delta_{ij}\equiv m_{i}^{2}-m_{j}^{2}$. A
combined analysis [6] of the data from various experiments indicates that the
across-generation mixing angle $\theta_{13}\sim 8.5^{\circ}$ is smaller than
$\theta_{12}\sim 34^{\circ}$, $\theta_{23}\sim 45^{\circ}$, while the CP phase
is still not well-established. There are two independent mass squared
differences, $|\delta_{31}|\gg\delta_{21}$. While $\delta_{21}\sim 7.6\times
10^{-5}$ eV2 is known to be positive, the sign of the other mass squared
difference is unknown, with $|\delta_{31}|\sim 2.5\times 10^{-3}$ eV2, with
the sum of all the neutrino masses constrained by cosmological considerations
to be less than about 2 eV. One of the important open issues to be settled by
future experiments is the so-called neutrino mass ordering, viz., whether the
third mass eigenstate is heavier or lighter than the others. We will discuss
the proposed INO in this context.

##  2 The INO project

The INO is a proposed mega-science project that is jointly funded by the
Department of Atomic Energy and the Department of Science and Technology,
Government of India. The immediate goal is the creation of an underground (at
least 1 km deep) cavern for scientific research purposes, with the main ICAL
detector to study atmospheric neutrinos and cosmic muons. It will incorporate
a centre for particle physics and detector technology at Madurai, South India,
which currently houses the mini-ICAL prototype. The INO graduate programme
trains students on both the theoretical and experimental aspects of
neutrino/particle physics.

###  2.1 The ICAL detector

The proposed [7] 51 kton ICAL detector will consist of 151 layers of 56 mm
thick soft iron that can be magnetised to about 15 kGauss. It will consist of
three identical modules with the layers of iron plates separated by 40 cm air
gaps in which the roughly 30,000 active Resistive Plate Chambers (RPCs) of
dimensions $2\times 2$ m2 will be housed; see Fig. 1. The high DC voltage
across the 3 mm glass plates separated by 2 mm of gas gap creates a discharge
when a charged particle passes through them. These are digitised and stored as
hits using the nearly 4 million channels of readout electronics, which are
later analysed to reconstruct the kinematics of each event. While the minimum
ionising muons leave clean long tracks in the detector, the hadrons shower.
Hence the momentum, including direction, of the muons is well-reconstructed
compared to that of the hadrons; in addition, the sign of the muon’s charge
can also be reconstructed to better than 98% for the few GeV events of
interest.

Figure 1: Schematic of ICAL detector and RPCs.

##  3 Atmospheric Neutrinos and ICAL

Cosmic rays interact with the Earth’s atmosphere to produce muons, which
further decay to neutrinos, producing muon and electron neutrinos (and anti-
neutrinos) in roughly 2:1 ratio. These neutrinos will mainly interact with the
iron in the detector to produce charged particles in charged-current (CC)
interactions:

$\displaystyle\nu_{\mu}+N$ $\displaystyle\to$ $\displaystyle\mu^{-}+X~{},$
$\displaystyle\overline{\nu}_{\mu}+N$ $\displaystyle\to$
$\displaystyle\mu^{+}+X~{},$ (5)

where $X$ is any hadronic debris. The magnetic field bends the differently
charged muons in different directions and hence differentiates between
neutrino and anti-neutrino induced events. This separation allows for a
sensitive study of matter effects which are different for neutrinos and anti-
neutrinos and hence can pin down the neutrino mass ordering.

###  3.1 ICAL Detector Response

The physics reach of ICAL using detailed GEANT4 based simulations [8] is shown
in this and the following sections. Unless specifically mentioned, all results
are from the INO white paper [7]; we give the main highlights here.

The magnetic field in ICAL allows for a separation of $\mu^{-}$ and $\mu^{+}$
(neutrino or anti-neutrino induced) events as well as helps to determine the
magnitude and direction of the muon’s momentum. A typical atmospheric neutrino
event is shown in Fig. 2. The lomg tracks of the mimimum ionising muons curve
in the presence of the magnetic field and their momentum is reconstructed
using a Kalman filter algorithm. The hadrons shower and appear as hits
clustered close to the vertex. Typically 10,000 events with fixed momentum and
direction were propagated in the simulated ICAL detector with GEANT, and their
behaviour analysed.

Figure 2: Typical event in ICAL simulations.

The muon momentum resolution (defined as $\sigma/P_{in}$, obtained by
reconstructing fixed momentum muons in ICAL and fitting the resultant
distribution to a gaussian distribution) is shown in Fig. 3 for different
values of $(\cos\theta,\phi)$, with the upward direction corresponding to
$\cos\theta=1$. The azimuthal dependence occurs due to the presence of the
magnetic field which breaks the azimuthal symmetry. (In addition, the angular
resolution (for the zenith angle) is better than $1^{\circ}$ for muons with
more than $\approx 2$ GeV energy.)

Figure 3: Muon resolution as a function of input momentum and $\cos\theta$,
for different regions of azimuthal angles, $|\phi|\leq\pi/4$ (top left),
$\pi/4\leq|\phi|\leq\pi/2$ (top right), $\pi/2\leq|\phi|\leq 3\pi/4$ (bottom
left), $3\pi/4\leq|\phi|\leq\pi$ (bottom right).

Fig. 4 shows the reconstruction $\epsilon_{rec}$ and relative charge-id
efficiency $\epsilon_{cid}$. The relative charge-id efficiency is the fraction
of reconstructed events that were reconstructed with the correct charge sign.
We see that $\epsilon_{cid}\gtrsim 98$% over a large range of momentum
magnitude and direction.

Figure 4: Reconstruction efficiency $\epsilon_{rec}$ (left) and charge-id
efficiency $\epsilon_{cid}$ (right) as a function of the input momentum for
different values of the zenith angle.

Finally, the hadron energy resolution for fixed energy charged pion events is
shown in Fig. 5. The final hadronic component in neutrino-iron interactions of
interest at ICAL can contain multiple hadrons; however, pions dominate the
sample; for details, please see Ref. [7].

Figure 5: Hadron energy resolution as a function of pion energy. Here the
reconstructed hadron energy was fitted to a Vavilov distribution and its
resolution $\sigma/E$ computed. See Ref. [7] for details.

###  3.2 Precision measurement of the 2–3 parameters

Honda atmospheric neutrino fluxes [9] were used to generate 1000 years of
events using the ICAL geometry and the NUANCE neutrino generator [10] and then
oscillated and scaled appropriately for the analysis. Typical events as a
function of the muon zenith angle $\cos\theta$ are shown in Fig. 6 for events
in the muon energy bin 2–3 GeV with an exposure time of 10 years. Finite
detector resolutions smear these numbers somewhat. Notice the poor
reconstruction in the horizontal direction.

Figure 6: Zenith angle distributions of oscillated atmospheric muon $\mu^{-}$
events with energies from 2 to 3 GeV for an exposure time of 10 years.
Statistical errors are shown. Note that there are about 2.5 times fewer
$\mu^{+}$ events due to the smaller cross sections although their fluxes are
approximately the same.

Fig. 7 shows the precision that can be obtained in 10 years in the 2–3
parameters, $\sin^{2}\theta_{23}$ and $\delta_{32}$ for input values of
$(0.5,2.4\times 10^{-3}$ eV${}^{2})$ in comparison with other experiments.
Note that ICAL is yet to be built! The results are marginalised over the
allowed $3\sigma$ range of $\sin^{2}\theta_{13}$ and include the reconstructed
muon momentum and direction information, as well as hadron energy information.
Results without hadron energy information are labelled 2D and those with
hadron energy included are labelled 3D in the figure.

Figure 7: Left: Precision reach in 10 years of ICAL for $\sin^{2}\theta_{23}$
and $\delta_{32}$, including hadron energy information. Right: Sensitivity to
the mass ordering for assumed input normal ordering as a function of exposure
time; similar results are obtained for inverted ordering. The sensitivity
increases with inclusion of larger muon energy range and inclusion of hadron
energy information.

Fig. 7 also shows the sensitivity to the mass ordering: plotted is the minimum
$\chi^{2}$ obtained when the data is fitted to the wrong (inverted) mass
ordering as a function of the exposure time when the true ordering is normal.
The result is marginalised over the magnitude of $\delta_{32}$ as well as
$\theta_{23}$ and $\theta_{13}$. The inclusion of hadron information improves
the result considerably. Similar results are obtained when the true hierarchy
is inverted.

The neutrino mass ordering is the centrepiece of ICAL physics. Although many
other experiments such as DUNE and JUNO will also measure this parameter, ICAL
is complementary to these other experiments because of the separation of the
neutrino and anti-neutrino events using the magnetic field. As a consequence,
the sensitivity is independent of the CP phase $\delta_{CP}$, unlike in other
experiments. Since this phase is currently unknown, the effects of this
parameter must be disentangled from those due to the mass ordering in the
other experiments. This can only be achieved for a fraction of the total range
of $\delta_{CP}$. This feature of ICAL can also be used to obtain synergies
with other experiments. For instance, the NO$\nu$A or ESS$\nu$SB sensitivity
to the mass ordering depends significantly on $\delta_{CP}$; this sensitivity
will be greatly improved on combining with ICAL data; see Fig. 8. The solid
line shows the sensitivity of NO$\nu$A to the mass ordering: the “data" is
generated with the normal ordering and then fitted assuming the inverted
ordering, to obtain the minimum value of $\chi^{2}$. The results are plotted
as a function of the true value of $\delta_{CP}$, and marginalised over all
oscillation parameters, including $\delta_{CP}$. It is seen that NO$\nu$A has
very little sensitivity to the mass ordering for values of
$\delta_{C}<180^{\circ}$. Hence the reach of NO$\nu$A alone for the neutrino
mass ordering is very sensitive to the actual value of $\delta_{CP}$. The
sensitivity in this small-$\delta_{CP}$ region improves somewhat (dotted
lines) on including reactor (especially T2K) data; and this improves
dramatically (dashed lines) on including the ICAL information.

Fig. 8 also shows that inclusion of ICAL data will increase the sensitivity to
$\delta_{CP}$ in regions of the $\delta_{CP}$ plane which earlier had poor
sensitivity to this parameter.

Figure 8: Left: Synergies with other experiments that improve the
determination of the mass ordering over the entire $\delta_{CP}$ range. Here
the full proposed runs of the long baseline and reactor experiments are taken,
with 10 years’ running of ICAL. Right: Synergies that improve the
determination of $\delta_{CP}$.

###  3.3 Other measurements at ICAL

Electrons from CC interactions of electron neutrinos (and anti-neutrinos) will
also produce showers in ICAL. It is not possible to separate the trackless
events arising from neutral current interactions of all flavours of neutrinos
in ICAL from these electron induced events. Low energy CC muon events where
the muon track is not reconstructed are an additional contamination. It turns
out that the trackless events also yield information [11] on the 2–3 neutrino
mixing parameters, see Fig. 9. Here, the solid line indicates the sensitivity
to $\sin^{2}\theta_{23}$ from just the trackless events at ICAL, when all
other oscillation parameters are kept fixed. The inclusion of systematic
uncertainties (dashed line) and further marginalising over the remaining
oscillation parameters (dotted line) reduces the sensitivity, which remains
significant.

Figure 9: Sensitivity to $\sin^{2}\theta_{23}$ from a simulations analysis of
“trackless" events. See text for details.

Finally, charged tau are produced in CC interactions of tau neutrinos with the
detector. Since atmospheric tau neutrinos are suppressed by a factor of about
$10^{5}$ compared to atmospheric muon neutrinos, significant yield of tau
neutrino events are a direct signature of neutrino oscillations, dominantly of
the atmospheric muon neutrinos. A simulations study of the tau neutrino-
induced CC events at ICAL, with the resulting taus decaying hadronically [12],
are also sensitive to neutrino oscillations. Since the threshold for CC tau
production is 3.5 GeV, such events are visible as an excess over the neutral
current (NC) background although they cannot be classified event-by-event. An
analysis of the combined NC and CC$\tau$ events shows sensitivity to both the
2–3 oscillation parameters $\theta_{23}$ and (the magnitude of) $\Delta
m_{31}^{2}$. Since CC muon events arise from the same atmospheric neutrino
fluxes, several systematic uncertainties are the same, and in fact, simulation
studies show that a combined study of muon and tau CC events improves the
sensitivity to the 2–3 oscillation parameters, in particular, the unknown
octant of the mixing angle $\theta_{23}$, as seen in Fig. 10.

Figure 10: Sensitivity to the octant of $\sin^{2}\theta_{23}$ from a combined
simulations analysis of CC muon and tau events, for an input value of
$\sin^{2}\theta_{23}=40^{\circ}$ (left) and $50^{\circ}$ (right).

There are many other studies of the sensitivity of ICAL to cosmic ray muons,
to BSM physics such as neutrino decays and sterile neutrinos, non-standard
interactions and Lorentz invariance violations, and the ability to search for
dark matter (WIMPS) and exotics such as monopoles.

##  4 Prototypes and the mini-ICAL detector

Cosmic muons have been observed at several RPC stacks that have been built for
calibration and testing. One such stack at Madurai in South India has measured
[13] the cosmic muon flux as a function of both the zenith angle and the
azimuthal angle and has also measured the east-west asymmetry of cosmic ray
muons; see Fig. 11.

Figure 11: Left: Measured east-west asymmetry of cosmic ray muon flux at the
Madurai RPC stack. Right: cosmic muon spectrum from the mini-ICAL detector at
Madurai.

The mini-ICAL detector; see Fig. 12, is a $4\times 4$ m2 85 ton scaled model
with 11 iron layers that has been constructed two years ago. After testing and
calibration, the cosmic muon spectrum has been measured with this detector and
is consistent with theoretical expectations; in addition, the $\mu^{-}$ and
$\mu^{+}$ events have been separately observed [14] and are shown as the
spectra with positive and negative momenta respectively in Fig. 11.

Figure 12: The prototype mini-ICAL detector at Madurai, South India.

##  5 Conclusion

Simulations studies show that ICAL has excellent physics potential, especially
to precisely measure the 2–3 neutrino mixing parameters as well as the
neutrino mass ordering. Detailed studies of charged current muon events
arising from interactions of atmospheric muon neutrinos have been performed;
studies of CC electron neutrino and tau neutrino events are on-going.
Separately, simulations studies have also examined the sensitivity of ICAL to
beyond-the-standard-model (BSM) physics, exotic particles such as WIMPS and
monopoles, etc. ICAL will be complementary in its reach to other detectors
around the world such as DUNE and JUNO. Finally, ICAL is a completely
indigenous detector with excellent R&D on RPC design and construction
including closed-loop gas system and automation by industry. In addition, the
detector electronics (analog and digital front end), including chip design and
development is also complete. A 700 ton large scale $8\times 8\times 2$ m3
engineering module is also to be constructed shortly.

In parallel, research is going on into neutrinoless double beta decay
experiments using tin bolometry, and dark matter experiments, which are also
proposed to be housed at INO. While it is expected that INO, when built, will
galvanise students as well as local industry due to cutting-edge technology
transfer, there have been delays in obtaining some technical clearances for
the underground lab construction. The INO collaboration is looking forward to
going ahead soon with this ambitious project.

## Acknowledgements

I thank the conference organisers and the members of the INO collaboration for
this opportunity to present the physics reach and status report of INO.

## References

  * [1] B. Pontecorvo,  _Mesonium and anti-mesonium_ ,  Sov. Phys. JETP 6, 429 (1957).
  * [2] Z. Maki, M. Nakagawa and S. Sakata,  _Remarks on the unified model of elementary particles_ ,  Prog. Theor. Phys. 28, 870 (1962),  10.1143/PTP.28.870.
  * [3] L. Wolfenstein,  _Neutrino Oscillations in Matter_ ,  Phys. Rev. D 17, 2369 (1978),  10.1103/PhysRevD.17.2369.
  * [4] S. P. Mikheev and A. Y. Smirnov,  _Resonant amplification of neutrino oscillations in matter and solar neutrino spectroscopy_ ,  Nuovo Cim. C 9, 17 (1986),  10.1007/BF02508049.
  * [5] S. Bilenky,  _Neutrino oscillations: From a historical perspective to the present status_ ,  Nucl. Phys. B 908, 2 (2016),  10.1016/j.nuclphysb.2016.01.025,  1602.00170.
  * [6] I. Esteban, M. C. Gonzalez-Garcia, M. Maltoni, T. Schwetz and A. Zhou,  _The fate of hints: updated global analysis of three-flavor neutrino oscillations_ ,  JHEP 09, 178 (2020),  10.1007/JHEP09(2020)178,  2007.14792.
  * [7] S. Ahmed _et al._ ,  _Physics Potential of the ICAL detector at the India-based Neutrino Observatory (INO)_ ,  Pramana 88(5), 79 (2017),  10.1007/s12043-017-1373-4,  1505.07380.
  * [8] S. Agostinelli _et al._ ,  _GEANT4–a simulation toolkit_ ,  Nucl. Instrum. Meth. A 506, 250 (2003),  10.1016/S0168-9002(03)01368-8.
  * [9] M. Honda, T. Kajita, K. Kasahara and S. Midorikawa,  _Improvement of low energy atmospheric neutrino flux calculation using the JAM nuclear interaction model_ ,  Phys. Rev. D 83, 123001 (2011),  10.1103/PhysRevD.83.123001,  1102.2688.
  * [10] D. Casper,  _The Nuance neutrino physics simulation, and the future_ ,  Nucl. Phys. B Proc. Suppl. 112, 161 (2002),  10.1016/S0920-5632(02)01756-5,  hep-ph/0208030.
  * [11] A. Chacko, D. Indumathi, J. F. Libby and P. K. Behera,  _First simulation study of trackless events in the INO-ICAL detector to probe the sensitivity to atmospheric neutrinos oscillation parameters_ ,  Phys. Rev. D 102(3), 032005 (2020),  10.1103/PhysRevD.102.032005,  1912.07898.
  * [12] R. T. Senthil, D. Indumathi and P. Shukla,  _A simulation study of tau neutrino events at the ICAL detector in INO_ (2022),  2203.09863.
  * [13] S. Pethuraj, G. Majumder, V. M. Datar, N. K. Mondal, K. C. Ravindran and B. Satyanarayana,  _Measurement of azimuthal dependent muon flux by 2 m x 2 m RPC stack at IICHEP-Madurai_ ,  Exper. Astron. 49(3), 141 (2020),  10.1007/s10686-020-09655-y,  1905.00739.
  * [14] A. Bhatt,  _Measurement of Atmospheric Muons at IICHEP in Madurai, for better estimation of Neutrino Fluxes at INO Site in Theni_ (2019).


# A non-absorbing SIR stochastic lattice gas model on hybrid lattices

Carlos Handrey Araujo Ferraz<EMAIL_ADDRESS>Exact and Natural Sciences
Center, Universidade Federal Rural do Semi-Árido-UFERSA, PO Box 0137, CEP
59625-900, Mossoró, RN, Brazil

###### Abstract

In this paper, we perform Monte Carlo calculations to study the critical
behavior of the spread of infectious diseases through a novel approach to the
SIR epidemiological model. A stochastic lattice gas version of the model was
applied on hybrid lattices which, in turn, are generated from typical square
lattices when inserting a connection probability $p$ that a given lattice site
has both first- and second-nearest-neighbor interactions. By combining
percolation theory and finite-size scaling analysis, we estimate both the
critical threshold and leading critical exponent ratios of the non-absorbing
SIR model in different cases of hybrid lattices. An analysis of the average
size of the percolating cluster of recovered individuals, as well as the size
distribution of non-percolating clusters, was carried out to determine the
universality class of the model.

###### keywords:

non-absorbing SIR model , hybrid lattices , percolating cluster , critical
exponent ratios , Monte Carlo simulation

##  1 Introduction

In the late 1920s, Kermack and McKendrick [1] proposed a system of ordinary
differential equations to determine the temporal evolution of the population
of individuals that interacted with each other when exposed to an infectious
disease. These individuals were initially divided into three classes:
susceptible (S), infected (I), and recovered (R) individuals. As an acronym
for these classes, the model became known as SIR model. Variations of this
model have been proposed over the years, such as the SIRS, SIS, SEIR models,
considering different classes of individuals that make up a given population.
The relevance of these approximations can not be underestimated, as such
approximations have been successful in modeling the most varied forms of
epidemics: cholera [2, 3], measles[4, 5], rubella [6, 7], hepatitis [8, 9],
influenza [10], AIDS [11, 12], COVID-19 [13, 14, 15], and many others.
Furthermore, epidemiological models are important tools to obtain information
on the speed of disease spreading, as well as to test protocols adopted by
public bodies to contain or mitigate such diseases, as they encompass both
features of collective dynamics [16, 17] and complex systems [18].

The SIR model is one of the simplest epidemiological models, where infected
individuals can be removed from the dynamic, either by permanent immunity or
by death. This peculiarity makes the model suitable for simulating epidemic
outbreaks such as influenza, SARS, AIDS, etc. It was shown that the SIR model
is in the same universality class as that of dynamic percolation [19, 20, 21].
In this paper, we consider a stochastic lattice-gas version of the model SIR
with asynchronous site updates. Similar models have been applied to other
population dynamics [22, 23]. For both synchronous and asynchronous versions,
it has been observed that a phase transition takes place when varying the
control parameters of the model. This transition is found to be of second
order between two distinct regimes: one in which the population remains
susceptible (inactive or endemic regime), and another in which the disease
spreads throughout the network (active or epidemic regime), where a
significant portion of the population becomes infected or eventually recovered
(immune or dead). At the transition point, the system becomes critical and
corresponds to an epidemic outbreak threshold. It is important to emphasize
that the connection between the stochastic and deterministic descriptions (via
coupled differential equations) can be achieved through a mean-field
approximation [24, 25] that results in the so-called Langevin equations
associated with the Fokker-Planck equation, which describes the temporal
evolution of the probability density of the states of the system.

The lattice gas method [26] used in the present study has been widely used to
address problems involving fluid dynamics [27, 28], general diffusion
processes [29], damage spreading [30, 31], and transport phenomena [32]; it
basically consists of a discretization both in time and space, where sites in
the network are seen as the most likely locations of the particle. In
molecular dynamics, for example, we only have a temporal discretization in the
description of the particle state. The lattice-gas method is especially useful
when one wants to describe the dynamics of a particles’ system without caring
about the detailed microscopic aspects of that system. These aspects are often
irrelevant to the description of the general behavior of the system. Here we
will employ this method as a way of describing the spread of the cloud of
recovered individuals over time as the epidemic evolves.

In particular, an analysis of the average size of the percolating cluster of
recovered individuals, as well as the size distribution of non-percolating
clusters, will be carried out through the use of the called Newman-Ziff
algorithm [33] in order to determine the critical exponent ratios of the
model. The SIR stochastic lattice gas model is said to be absorbing-type
because its active phase is characterized by an infinite number of absorbing
configurations, where the final state of the system is composed only by
recovered individuals. However, in the present work, the simulations are
halted as long as the existence of the percolating cluster (spanning cluster)
in the system is verified, thus setting the non-absorbing state of the model.
This procedure speeds up the analysis of simulation data without compromising
the guarantee that the asymptotic limit of the system has been reached. It is
known that only a single percolation cluster is generated in the critical
regime of the system. By exploiting the analysis of these clusters, we will be
able to determine the universality class of the model.

Since the network topology represents an important feature of the system and
directly influences the dynamics of the processes involved [30, 34, 31], we
will study here a hybrid network topology. Hybrid lattices are formed from
regular square lattices ($N=L\times L$) by inserting a probability $p$ that a
given site has both first- and second-nearest-neighbor interactions. Such
networks would simulate a more realistic population formed mainly by two types
of individuals: type I, with low connectivity, and type II, with high
connectivity. The special cases with $p=0$ and $p=1$ are also treated here and
correspond to pure lattices having all nodes with a connectivity of first
neighbors and first and second neighbors, respectively.

The main goal of this study is to understand how extended connectivity effects
can influence the critical behavior of the non-absorbing SIR model. We can
accomplish that by estimating both the epidemic thresholds and the leading
critical exponent ratios of the model. In the present paper, we deal with the
non-absorbing SIR model coupled to square lattices with first-neighbor
interactions (case $p=0$), first- and second-neighbor interactions (case
$p=1$), and to hybrid lattices (case $p=1/2$).

The contents of the article are organized as follow. In section 2, we outline
the SIR stochastic lattice gas model. In section 3, we describe details of the
Monte Carlo (MC) simulation background and lattice generation. In section 4,
we present and discuss the results. Finally, in section 5, we make the
conclusions.

##  2 SIR stochastic lattice gas model

(a) Case $p=0$

(b) Case $p=1/2$

Figure 1: Typical spanning clusters along with a few smaller clusters formed
close to the epidemic threshold $\lambda_{c}$ on a $100\times 100$ lattice for
two different $p$ cases. (a) Case $p=0$ (at $\lambda=0.16$) and (b) Case
$p=1/2$ (at $\lambda=0.21$). Susceptible, infected, and recovered individuals
are represented by blue, red, and green vertexes, respectively.

The SIR stochastic lattice gas model is defined on a given lattice of N sites
in which each site can be occupied by just one individual that can be either a
susceptible (state $S$), an infected (state $I$), or a recovered
(immunized/dead) individual (state $R$). The dynamics consists of two
subprocesses, namely, an auto-catalytic one, $I+S\rightarrow I+I$; and a
spontaneous one, $I\rightarrow R$. At each time step a site is randomly chosen
and then a set of dynamic rules are taken in the following way:
i) If the site is in the state $S$ and there is at least one neighboring site
in the sate $I$ then the site becomes $I$ with probability proportional to a
parameter $\mu$ and the number $z$ of neighboring sites, i.e., $\mu{\kern
1.0pt}m/z$, where $m$ is the number of neighboring $I$ sites.
ii) If the site is in state $I$ it becomes $R$ spontaneously with probability
$\lambda$.
iii) If the site is $R$ it remains unchanged.
At each site $i$ of a 2D lattice we assign a stochastic variable $\sigma_{i}$
that takes the values 0, 1 or 2, according to whether the site is in the state
$S$, $I$, or $R$, respectively. Since the transitions between states in this
model are nonequilibrium ones, the allowed transitions of the state $i$ of a
site are cyclic, this is, $0\rightarrow 1\rightarrow 2$. The corresponding
transition rate is represented by $w_{i}(\sigma)$ and describes the transition
$\sigma\rightarrow\sigma^{\prime}$ in which the whole microscopic
configuration (microstate)
$\sigma^{\prime}\equiv(\sigma_{1},\ldots,\sigma^{\prime}_{i},\ldots,\sigma_{N})$
differs from $\sigma$ only by the state of $i$-th site. It is given by

$w_{i}(\sigma)=\frac{\mu}{z}\delta(\sigma_{i},0)\sum\limits_{j}{\delta(\sigma_{j},1)}+\lambda{\kern
1.0pt}\delta(\sigma_{i},1),$ (1)

where the summation runs over the nearest neighbors of site $i$ and
$\delta(x,y)$ denotes the Kronecker delta. The parameters $\mu$ and $\lambda$
are related to the subprocesses above described, and are chosen such that
$\mu+\lambda=1$.

The system evolves in time according to a master equation for the probability
distribution $P(\sigma,t)$ described by

$\frac{d}{{dt}}P(\sigma,t)=\sum\limits_{i}{\\{w_{i}(\bar{\sigma})}P(\bar{\sigma},t)-w_{i}(\sigma)P(\sigma,t)\\},$
(2)

where the microstate $\bar{\sigma}$ is obtained from $\sigma$ by an anticyclic
permutation of the state of site $i$ ($2\rightarrow 1\rightarrow 0$).

##  3  Monte Carlo Simulation and Lattice Generation

Table 1:  Estimates of the epidemic threshold $\lambda_{c}$ and critical exponents ratios for each $p$ cases. Lattice | Epidemic threshold | $1/\nu$ | $\beta/\nu$ | $\gamma/\nu$
---|---|---|---|---
Square with $p=0$ | $\lambda_{c}=0.176(6)$ | $0.779\pm 0.049$ | $0.105\pm 0.002$ | $1.765\pm 0.006$
Hybrid with $p=1/2$ | $\lambda_{c}=0.228(4)$ | $0.663\pm 0.029$ | $0.097\pm 0.003$ | $1.770\pm 0.006$
Square with $p=1$ | $\lambda_{c}=0.275(0)$ | $0.683\pm 0.046$ | $0.084\pm 0.004$ | $1.809\pm 0.022$
Exact values | $-$ | $3/4$ | $5/48$ | $43/24$

We can implement an asynchronous, non-absorbing SIR model on computer by
following the kinetic Monte Carlo rules below:

  1. 1.

First, we start with a single central infected site (seed) and the remaining
ones being all susceptible on a two-dimensional lattice in which each
individual of the population is attached to its respective lattice site. In
order to speed up the simulation we create two lists that are updated at each
algorithm step: a list of infected individuals (infected list) and a list of
recovered individuals (recovered list), which begins empty.

  2. 2.

Next, we update the system state by randomly choosing an available infected
site from the infected list and proceed as follows:

    1. (a)

Generate a random number $x$ in the interval $(0,1)$. If $x\leq\lambda$, the
infected site is removed from the infected list and placed in the recovered
list;

    2. (b)

Otherwise (if $x>\lambda$), pick randomly one nearest neighbor of the infected
site and make it also infected provided that it is susceptible, adding it to
infected list.

  3. 3.

Repeat asynchronously the step (2) several times until either there is no
infected sites (endemic phase) or there is a percolating cluster of recovered
sites (non-absorbing epidemic phase).

One could determine the MC time $t$ by incrementing $t$ by $\delta t=1/n_{I}$,
where $n_{I}$ is the current number of infected sites, each time an infected
site is pick from the list. However, we do not keep track of time here since
we are more interested in static quantities such as the fraction and the mean
cluster size of $R$ sites.

Remarkably, it was shown that the SIR model on square lattices belongs to the
same universality class then dynamic percolation (DP). This allows us to
investigate the phase transition which takes place in the present non-
absorbing SIR model by making use of the percolation theory. Thus we can
define the epidemic phase of the model when it is observed a formed
percolating cluster of recovery sites in the system and the endemic phase when
it is not. Such graphical analogy has been used for other compartmental models
as well. Following the classical percolation theory it is important first to
determine the cluster distribution of recovery sites, i.e., the number of
clusters with $s$ recovery sites $n_{p}(s)$. That can be accomplished by using
the Newmann-Ziff algorithm, which possess also the built-in feature of
identifying whether a percolating cluster was formed or not, depending on the
$\lambda$ value taken for the model. Notice that for systems with non-periodic
boundaries like the ones concerned here, the percolating cluster is actually a
spanning cluster [33, 35].

From the cluster size distribution, we have the fraction of recovery sites in
the finite (non-percolating) cluster with $s$ size

$P_{s}=s\frac{{n_{p}(s)}}{{n_{R}}},$ (3)

where $n_{R}$ is the total number of recovery sites and $n_{p}(s)$ is the
number of clusters with $s$ recovery sites. Furthermore, the fraction of
recovery sites in the percolating cluster $P_{\infty}$ can be obtained by

$P_{\infty}=1-\frac{1}{{n_{R}}}\sum\limits_{s}{s{\kern 1.0pt}n_{p}(s)},$ (4)

such that the above summation excludes the percolation cluster. Now we can
define the order parameter from Eq. (4) as

$P=<P_{\infty}>,$ (5)

where $<x>$ means an average taken over different dynamic realizations. The
epidemic phase of the model is reached when $P\neq 0$, that is, when the
percolating cluster density is non-zero; while the endemic phase happens when
$P=0$. Other important quantities are the mean cluster size

$S=\frac{1}{{n_{R}}}\sum\limits_{s}{s^{2}{\kern 1.0pt}n_{p}(s)},$ (6)

which plays the rule of the susceptibility in classical percolation theory
[36, 37, 38] when taking the average over different runs, i.e.,

$\chi=<S>,$ (7)

the overall mean cluster size

$S^{\prime}=\frac{1}{{n_{R}}}\sideset{}{{}^{\prime}}{\sum}_{s}{s^{2}{\kern
1.0pt}n_{p}(s)},$ (8)

and the mean quadratic cluster size

$M^{\prime}=\frac{1}{{n_{R}}}\sideset{}{{}^{\prime}}{\sum}_{s}{s^{3}{\kern
1.0pt}n_{p}(s)},$ (9)

where the primed summations above also include the percolating cluster. It is
worth remarking that the last two quantities $S^{\prime}$ and $M^{\prime}$
only make sense for finite lattice as in the asymptotic limit
$(N\rightarrow\infty)$, the percolating cluster size diverges.

At criticality, the cluster size distribution should obey a power-law scaling
[39, 37, 20] as

$n_{p}(s)=s^{-\tau}F[s^{\alpha}\;(\lambda-\lambda_{c})],$ (10)

where $\lambda_{c}$ is the epidemic threshold and $F$ is a scaling function.
Similarly, the remaining quantities also obey scaling relations in accordance
to classical percolation theory given by

$\displaystyle P$ $\displaystyle=L^{-\beta/\nu\;}\tilde{P}(L^{1/\nu}{\kern
1.0pt}|\lambda-\lambda_{c}|),$ (11) $\displaystyle\chi$
$\displaystyle=L^{\gamma/\nu\;}\tilde{\chi}(L^{1/\nu}{\kern
1.0pt}|\lambda-\lambda_{c}|),$ (12) $\displaystyle<S^{\prime}>$
$\displaystyle=L^{\gamma/\nu\;}\tilde{S}^{\prime}(L^{1/\nu}{\kern
1.0pt}|\lambda-\lambda_{c}|),$ (13) $\displaystyle<M^{\prime}>$
$\displaystyle=L^{(\beta+2\gamma)/\nu\;}\tilde{M}^{\prime}(L^{1/\nu}{\kern
1.0pt}|\lambda-\lambda_{c}|).$ (14)

In addition, we can define a universal quantity $U$ in which the scaling
dependencies cancel out by combining Eqs. (11), (13), and (14) in the
following way [20]

$U=P\frac{{<M^{\prime}>}}{{<S^{\prime}>^{2}}},$ (15)

being analogous to the Binder cumulant for ferromagnetic spin model [40, 41],
and obeying also a scaling relation

$U=\tilde{U}\;(L^{1/\nu}{\kern 1.0pt}|\lambda-\lambda_{c}|).$ (16)

The crossing point of the $U$ curves for different lattice sizes allow us to
estimate the epidemic threshold $\lambda_{c}$, whereas a finite-size scaling
analysis of the observables $P$ and $\chi$ by using Eqs. (11) and (12) yields
the according critical exponent ratios $\beta/\nu$and $\gamma/\nu$. The
leading critical exponents $\beta$, $\gamma$ and $\nu$ define the universality
class of the system.

The hybrid lattices used in the present study were constructed starting from
regular square lattices with free boundary conditions. First, we begin from a
regular square lattice consisting of nodes linked to their four first nearest
neighbors by both outgoing and incoming links. Then, with probability $p$, we
connect a chosen site also to their second nearest neighbors. After repeating
this procedure for every site, a new lattice is constructed with a density $p$
of nodes with both first- and second-nearest-neighbor connections. Such
networks would mimic a hypothetical population formed basically by two types
of individuals: type I, with low connectivity, and type II, with high
connectivity. The special cases with $p=0$ and $p=1$ correspond to pure
lattices having all nodes with connectivity of first neighbors and first and
second neighbors, respectively. Figs. 1(a) and 1(b) display typical spanning
clusters along with a few smaller clusters formed close to $\lambda_{c}$ for
the cases with $p=0$ and $p=1/2$, respectively. These clusters arise from the
SIR model dynamics. The spanning cluster is formed when the cloud of recovered
individuals reaches any two opposing edges of the lattice, in other words,
when such a cloud spans the lattice from one side to the other. We grew more
than $10^{5}$ spanning clusters for every considered $\lambda$ value to take
reliable averages of the above quantities. On average it took about 26 ms to
grow a single spanning cluster like those shown in Fig. 1 on an 3.70 GHz Intel
Xeon workstation.

##  4  Results and Discussion

Figure 2: Static quantities for the non-absorbing SIR stochastic lattice gas
model on square lattices with first-neighbor interactions (case $p=0$). Panels
(a), (b), and (c) display the order parameter $P$, Binder cumulant $U$, and
susceptibility $\chi$ as a function of the recovery rate $\lambda$,
respectively. From the Binder cumulant crossing, we can estimate the epidemic
threshold at $\lambda_{c}=0.176(6)$. The inset in panel (b) is a refinement of
the calculations for $U$ inside the critical region. Panels (d), (e), and (f)
display the log-log plot of the quantities $\phi$, $P$, and $\chi$ calculated
at $\lambda_{c}$ as a function of the linear size of the system $L$,
respectively. Red straight lines are the best linear fit to the corresponding
data.

Figure 3: Static quantities for the non-absorbing SIR stochastic lattice gas
model on hybrid lattices (case $p=1/2$). Panels (a), (b), and (c) display the
order parameter $P$, Binder cumulant $U$, and susceptibility $\chi$ as a
function of the recovery rate $\lambda$, respectively. From the Binder
cumulant crossing, we can estimate the epidemic threshold at
$\lambda_{c}=0.228(4)$. The inset in panel (b) is a refinement of the
calculations for $U$ inside the critical region. Panels (d), (e), and (f)
display the log-log plot of the quantities $\phi$, $P$, and $\chi$ calculated
at $\lambda_{c}$ as a function of the linear size of the system $L$,
respectively. Red straight lines are the best linear fit to the corresponding
data.

Figure 4: Static quantities for the non-absorbing SIR stochastic lattice gas
model on square lattices with first- and second-neighbor interactions (case
$p=1$). Panels (a), (b), and (c) display the order parameter $P$, Binder
cumulant $U$, and susceptibility $\chi$ as a function of the recovery rate
$\lambda$, respectively. From the Binder cumulant crossing, we can estimate
the epidemic threshold at $\lambda_{c}=0.275(0)$. The inset in panel (b) is a
refinement of the calculations for $U$ inside the critical region. Panels (d),
(e), and (f) show the log-log plot of the quantities $\phi$, $P$, and $\chi$
calculated at $\lambda_{c}$ as a function of the linear size of the system
$L$, respectively. Red straight lines are the best linear fit to the
corresponding data.

In this section, we show our numerical results of the non-absorbing SIR model
coupled to hybrid lattices. In order to determine both the critical region and
the order of the phase transition in this model on hybrid lattices, we
calculated the order parameter $P$, Binder cumulant $U$, and susceptibility
$\chi$ for each $p$ case in a wide range of the parameter $\lambda$. These
quantities were averaged over at least $10^{5}$ different dynamic realizations
of the SIR model. Furthermore, for the case $p=1/2$ we consider also different
lattice configurations upon taking the averages. We deal with several lattice
sizes, ranging from $N=400$ up to $N=10000$. Let us first discuss the case
$p=0$.

###  4.1  Case $p=0$

In Figs. 2, 2 and 2 are shown the order parameter, Binder cumulant, and
susceptibility as a function of the recovery rate $\lambda$ for the case
$p=0$, respectively. As one can see from Fig. 2, a typical second-order phase
transition takes place for the case $p=0$. As we shall see, the same
conclusion can be drawn for the remaining cases, where a typical sigmoid-
shaped curve also occurs. From Binder cumulant crossings, we can estimate the
corresponding epidemic thresholds for each case. In the inset of Fig. 2 is
shown a refinement of the calculations for $U$ inside the critical region. The
critical thresholds were estimated with five significant figures. For the case
$p=0$, we obtained $\lambda_{c}=0.176(6)$. As expected, it is approximately
equal to that observed in the absorbing version of the SIR model. To estimate
the reciprocal correlation-length exponents $1/\nu$ for each case, we
calculated the modulus of the logarithmic derivative of $P$ at the critical
threshold point $\lambda_{C}$

$\phi\equiv\bigg{|}\frac{d}{{d\lambda}}\ln(P)\bigg{|}_{\lambda=\lambda_{c}},$
(17)

where the derivative of the function $f\equiv\ln(P)$ was evaluated numerically
by using a finite central difference scheme in the form

$\frac{{df}}{{d\lambda}}\simeq\frac{1}{{2h}}\left({f(\lambda+h)-f(\lambda-h)}\right),$
(18)

which has an truncation error of the order of $O(h)^{2}$. The error function
$\delta f^{\prime}$ of $df/d\lambda$ was obtained via error propagation from
the uncertainties in the values of $f$ ($\delta f$), being expressed by

$\delta f^{\prime}=\frac{1}{{2h}}\sqrt{(\delta f(\lambda+h))^{2}+(\delta
f(\lambda-h))^{2}}.$ (19)

In our computations $h$ was taken equal to $2.0\times 10^{-3}$. By taking the
slope of the log-log plot of the quantity $\phi$ versus the linear size of the
system $L$, we can get an estimate for $1/\nu$. Fig. 2 shows the best linear
fit to data for the case $p=0$. We obtained $1/\nu=0.74\pm 0.05$. Error bars
were estimated by using Eq. (19). This value of $1/\nu$ is within one standard
deviation from the exact critical exponent ratio $1/\nu=3/4$ of 2D dynamical
percolation and it deviates only 1% from this exact value. Similarly, a
finite-size scaling analysis of the magnitudes of the order parameter $P$ and
the susceptibility $\chi$ at $\lambda_{c}$ by using Eqs. (11) and (12) yield,
respectively, the critical exponent ratios $\beta/\nu$ and $\gamma/\nu$. Figs.
2 and 2 show the log-log plot of $P$ and $\chi$ (both calculated at
$\lambda_{c}$) against $L$, respectively. The red straight lines in those
figures are the best linear fit to Eqs. (11) and (12), respectively. We
obtained $\beta/\nu=0.105\pm 0.002$ and $\gamma/\nu=1.765\pm 0.006$. These
values are also in very good agreement with the exact critical exponent ratios
of 2D dynamical percolation, namely, $\beta/\nu=3/4$ and $\gamma/\nu=43/24$.
These estimates of the critical exponents ratios and critical threshold
$\lambda_{c}$ for $p=0$ case are found also in Table 1.

###  4.2  Cases $p\neq 0$

An analogous analysis can be done for the cases $p=1/2$ and $p=1.0$. Figs. 3
and 3 display the order parameter $P$ as a function of the recovery rate
$\lambda$ for the cases $p=1/2$ and $p=1$, respectively. As already remarked,
both systems undergo also a second-order transition with their respective $P$
curves exhibiting a typical sigmoidal shape. While Figs. 3 and 3 show the
Binder cumulant for $p=1/2$ and $p=1$, respectively. From these figures, one
see that the crossing points are located at $\lambda_{c}=0.228(4)$ for $p=1/2$
and $\lambda_{c}=0.275(0)$ for $p=1$. Likewise, we took the slope of the log-
log of $\phi$ defined by Eq. (17) versus $L$ to estimate the exponent ratio
$1/\nu$ for $p=1/2$ and $p=1$ cases. Figs. 3 and 4 display the linear curve
fitting to the data for $p=1/2$ and $p=1$ cases, respectively. Error bars were
calculated by using Eq. (19). We got $1/\nu=0.663\pm 0.029$ for $p=1/2$ and
$1/\nu=0.683\pm 0.046$ for $p=1$. Such results strongly suggest that both
cases are not in the same universality class of 2D dynamic percolation.
Moreover, these results are fairly close to each other, as one falls within
only one standard deviation from the other. Similarly, we took the slopes of
the log-log plots of $P$ and $\chi$ versus $L$ for both cases. Again, the red
straight lines are linear regressions to the corresponding data. From Figs. 3
and 4, we obtained $\beta/\nu=0.097\pm 0.003$ for $p=1/2$ and
$\beta/\nu=0.084\pm 0.004$ for $p=1$, respectively. While from Figs. 3 and 4,
we got $\gamma/\nu=1.770\pm 0.006$ for $p=1/2$ and $\gamma/\nu=1.809\pm 0.022$
for $p=1$, respectively. Except for the estimates of $\beta/\nu$, one can see
that the estimates of $\gamma/\nu$ for both case are rather close to each
other. The estimates of the critical exponent ratios and the critical
threshold $\lambda_{c}$ for all treated $p$ cases are summarized in Table 1.
These results clearly indicate that the non-absorbing SIR model both on hybrid
lattices ($p=1/2$) and on square lattices with first- and second-neighbor
interactions (case $p=1$) belong each one to a new universality class
different from that of 2D dynamic percolation. In addition, it was found that
the critical threshold in the non-absorbing SIR model increases with
increasing lattice connectivity.

##  5  Conclusions

We performed Monte Carlo simulations of the non-absorbing SIR stochastic
lattice gas model on hybrid lattices to study the critical behavior presented
by these systems. Both the critical threshold and leading critical exponent
ratios were estimated for different cases. Our analysis has revealed that for
the case with connection probability equal to zero ($p=0$), the system falls
into the same universality class as that of two-dimensional dynamical
percolation. Besides that, the critical threshold for this case is found to be
the same as that of the absorbing version of the model. While for the cases
with connection probability nonzero, held on hybrid lattices ($p=1/2$) and on
square lattices with first- and second-neighbor interactions ($p=1$), belong
each one to a new universality class different from that of two-dimensional
dynamic percolation. Furthermore, it was found that the critical threshold in
the non-absorbing SIR model increases with increasing lattice connectivity.

This study has wide applications to many issues, including not only the spread
of infectious diseases, but also in general diffusion processes, damage
propagation in random networks, and performance optimization in multi-core
architectures. Finally, we expect that the results presented in this paper can
be helpful to understand how increased connectivity effects can influence the
critical properties of more complex systems.

## References

  * Kermack and Mckendrick [1927] W. . Kermack, A. G. Mckendrick, Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character 115 (1927) 700–721.
  * Hartley et al. [2005] D. M. Hartley, J. G. M. Jr, D. L. Smith, PLOS Medicine 3 (2005) e7.
  * Pascual et al. [2006] M. Pascual, K. Koelle, A. P. Dobson, PLoS Medicine 3 (2006) 0931–0933.
  * Kassem and Ndam [2010] T. Kassem, J. Ndam, Science World Journal 3 (2010) 29–32.
  * Tilahun et al. [2020] G. T. Tilahun, S. Demie, A. Eyob, Infectious Disease Modelling 5 (2020) 478–494.
  * Amaku et al. [2003] M. Amaku, F. A. B. Coutinho, R. S. Azevedo, M. N. Burattini, L. F. Lopez, E. Massad, Physical Review E 67 (2003) 051907.
  * Prawoto et al. [2020] B. P. Prawoto, Abadi, R. Artiono, AIP Conference Proceedings 2264 (2020) 020005.
  * Zou et al. [2010] L. Zou, W. Zhang, S. Ruan, Journal of Theoretical Biology 262 (2010) 330–338.
  * Khan et al. [2013] M. A. Khan, S. Islam, M. Arif, Z. U. Haq, BioMed Research International 2013 (2013) 1–10.
  * Kim et al. [2020] Y. Kim, A. V. Barber, S. Lee, PLOS ONE 15 (2020) e0232580.
  * Huang et al. [2011] Z. Huang, Q. Yang, J. Cao, Chaos, Solitons & Fractals 44 (2011) 954–963.
  * Wang et al. [2017] Y. Wang, D. Jiang, T. Hayat, B. Ahmad, Applied Mathematics and Computation 315 (2017) 477–493.
  * Atkeson et al. [2020] A. Atkeson, K. Kopecky, T. Zha, NBER Working Paper No. 27335 (2020).
  * Yang et al. [2021] H. M. Yang, L. P. L. Junior, F. F. M. Castro, A. C. Yang, PLOS ONE 16 (2021) e0252271.
  * Gounane et al. [2021] S. Gounane, Y. Barkouch, A. Atlas, M. Bendahmane, F. Karami, D. Meskine, Epidemiologic Methods 10 (2021) 20200044.
  * Wang et al. [2016] Z. Wang, C. T. Bauch, S. Bhattacharyya, A. d’Onofrio, P. Manfredi, M. Perc, N. Perra, M. Salathé, D. Zhao, Physics Reports 664 (2016) 1–113.
  * Shao and Han [2022] Q. Shao, D. Han, Applied Mathematics and Computation 412 (2022) 126559\.
  * Helbing et al. [2015] D. Helbing, D. Brockmann, T. Chadefaux, K. Donnay, U. Blanke, O. Woolley-Meza, M. Moussaid, A. Johansson, J. Krause, S. Schutte, M. Perc, Journal of Statistical Physics 158 (2015) 735–781.
  * Tomé and Ziff [2010] T. Tomé, R. M. Ziff, Physical Review E 82 (2010) 051921.
  * Souza et al. [2011] D. R. D. Souza, T. Tomé, R. M. Ziff, Journal of Statistical Mechanics: Theory and Experiment 2011 (2011) P03006.
  * Pastor-Satorras et al. [2015] R. Pastor-Satorras, C. Castellano, P. V. Mieghem, A. Vespignani, Reviews of Modern Physics 87 (2015) 925.
  * Satulovsky and Tomé [1994] J. E. Satulovsky, T. Tomé, Physical Review E 49 (1994) 5073.
  * Antal et al. [2001] T. Antal, M. Droz, A. Lipowski, G. Ódor, Physical Review E 64 (2001) 036118.
  * Lugo and McKane [2008] C. A. Lugo, A. J. McKane, Physical Review E 78 (2008) 051911.
  * de Souza and Tomé [2010] D. R. de Souza, T. Tomé, Physica A: Statistical Mechanics and its Applications 389 (2010) 1142–1150.
  * Sukop and Thorne [2007] M. C. Sukop, D. T. Thorne, in: Lattice Boltzmann Modeling, Springer, Berlin, Heidelberg, 2007, pp. 13–26.
  * Frisch et al. [1986] U. Frisch, B. Hasslacher, Y. Pomeau, Physical Review Letters 56 (1986) 1505.
  * D’humieres and Lallemand [1987] D. D’humieres, P. Lallemand, Complex Systems 1 (1987) 599–632.
  * Kehr and Binder [1984] K. W. Kehr, K. Binder, Simulation of Diffusion in Lattice Gases and Related Kinetic Phenomena, Springer, Berlin, Heidelberg, 1984\.
  * Ferraz and Herrmann [2008] C. H. A. Ferraz, H. J. Herrmann, Physica A: Statistical Mechanics and its Applications 387 (2008) 5689–5695.
  * Ferraz and Lima [2018] C. H. A. Ferraz, J. L. S. Lima, Brazilian Journal of Physics 48 (2018) 531–538.
  * Deutsch [1996] A. Deutsch, Nonlinear Physics of Complex Systems (1996) 353–362.
  * Newman and Ziff [2001] M. E. J. Newman, R. M. Ziff, Physical Review E 64 (2001) 016706.
  * Ferraz and Lima [2017] C. H. A. Ferraz, J. L. S. Lima, Physica A: Statistical Mechanics and its Applications 484 (2017) 488–498.
  * Sen [2011] P. Sen, International Journal of Modern Physics C 8 (2011) 229–237.
  * Shante and Kirkpatrick [1971] V. K. Shante, S. Kirkpatrick, Advances in Physics 20 (1971) 325–357.
  * Hoshen et al. [1979] J. Hoshen, D. Stauffer, G. H. Bishop, R. J. Harrison, G. D. Quinn, Journal of Physics A: Mathematical and General 12 (1979) 1285\.
  * Stauffer and Aharony [2014] D. Stauffer, A. Aharony, Introduction to Percolation Theory, Taylor & Francis Group, 2014.
  * Stauffer [1979] D. Stauffer, Physics Reports 54 (1979) 1–74.
  * Selke and Shchur [2005] W. Selke, L. N. Shchur, Journal of Physics A: Mathematical and General 38 (2005) L739.
  * Ferraz [2015] C. H. A. Ferraz, Physica A: Statistical Mechanics and its Applications 440 (2015) 90–99.