††thanks: On leave of absence from Rzhanov Institute of Semiconductor Physics
and the Novosibirsk State University, Russia.

# Parametrization of the Charge-Carrier Mobility in Organic Disordered
Semiconductors. APAE against EGDM.

S. D. Baranovskii  Department of Physics and Material Sciences Center,
Philipps-University, D-35032 Marburg, Germany  Department für Chemie,
Universität zu Köln, Luxemburger Strasse 116, 50939 Köln, Germany  A. V.
Nenashev  Department of Physics and Material Sciences Center, Philipps-
University, D-35032 Marburg, Germany  D. Hertel  Department für Chemie,
Universität zu Köln, Luxemburger Strasse 116, 50939 Köln, Germany  K. Meerholz
Department für Chemie, Universität zu Köln, Luxemburger Strasse 116, 50939
Köln, Germany  F. Gebhard  Department of Physics and Material Sciences Center,
Philipps-University, D-35032 Marburg, Germany

###### Abstract

An appropriately parameterized compact analytical equation (APAE) is suggested
to account for charge carrier mobility in organic disordered semiconductors
(ODSs). This equation correctly reproduces the effects of temperature $T$,
carrier concentration $n$, and electric field $F$ on the carrier mobility
$\mu(T,F,n)$, as evidenced by comparison with analytical theories and Monte
Carlo simulations. The set of material parameters responsible for charge
transport is proven to be at varience to those used in the so-called extended
Gaussian disorder model (EGDM) approach, which is widely exploited in
commercially distributed device–simulation algorithms. While EGDM is only
valid for cubic lattices with a specific choice of parameters, APAE describes
charge transport in systems with spatial disorder in a wide range of
parameters. APAE is user-friendly and, thus, suitable for incorporation into
device-simulation algorithms.

##  I Introduction

A correct parametrization of carrier mobility in organic disordered
semiconductors (ODSs) with hopping charge transport is of vital importance for
the development of algorithms desired to simulate devices based on ODSs, such
as organic light–emitting diodes, organic solar cells, and organic field-
effect transistors. Theoretical equations for the charge carrier mobility
$\mu(T,n,F)$, dependent on the concentration of carriers $n$, on temperature
$T$, and on the applied electric field $F$, are at the heart of the device
simulation algorithms. The choice of the appropriate theoretical description
for $\mu(T,n,F)$ suitable for applications in the device simulation software
has been recently addressed in several studies Upreti  _et al._ (2019); Lee
_et al._ (2021); Sun  _et al._ (2021); Oelerich  _et al._ (2012); Nenashev
_et al._ (2017); Oelerich  _et al._ (2017). In spite of the progress achieved
in these studies, the parametrization of the dependences $\mu(T,n,F)$ still
needs improvement.

So far, most of the device simulation algorithms, including commercially
distributed software packages, are based on the so-called EGDM equation
Pasveer  _et al._ (2005). This equation has been suggested Pasveer  _et al._
(2005) as a “unified description of charge-carrier mobilities in disordered
semiconducting polymers”. Several review papers Coehoorn and Bobbert (2012);
Kuik  _et al._ (2014) promote the EGDM equation as gold standard for the
description of carrier mobilities in ODSs, and the EGDM is the basis of
commercially distributed device-simulation packages Coehoorn and Bobbert
(2012); Lee  _et al._ (2021).

However, the EGDM equation Pasveer  _et al._ (2005) opposes the basic
theoretical concepts developed so far for hopping transport in disordered
materials Mott (1969); Shklovskii (1973); Baranovskii (2018). First, the EGDM
equation is based on an irrelevant parametrization Mott (1969); Shklovskii
(1973); Baranovskii (2018). Some parameters responsible for $\mu(T,n,F)$ are
missed in the EGDM, while some parameters present in the EGDM are not
responsible for the effects. Second, the EGDM equation was formulated Pasveer
_et al._ (2005) to fit simulation data on regular cubic lattices without
spatial disorder. However, hopping mobility on regular lattices is known to
deviate significantly from that in materials with spatial disorder Oelerich
_et al._ (2017); Baranovskii (2018).

Therefore, it is necessary to figure out whether the EGDM based on regular
latices with deficient parametrization could be of use for systems with
spatial disorder.

Regular cubic lattice used in EGDM is a specific case of the Gaussian-
disorder-model (GDM), in which charge transport is due to incoherent hopping
of carriers via random spatially distributed localized states with Gaussian
energy spectrum. This model seems valid by its ability to account for
experimentally observed dependences $\mu(T,n,F)$ Bässler (1993); Schmechel
(2002); Tessler  _et al._ (2009). A parametrization of the carrier mobility
$\mu(T,n,F)$ in the framework of GDM is given in Sec. II.

While a lot of effort has been focused on computer simulations, the transport
problem in the GDM can be easily solved analytically in the form of a simple
closed-form system of equations Baranovskii  _et al._ (2000, 2002); Rubel  _et
al._ (2004); Baranovski (2006); Oelerich  _et al._ (2012); Baranovskii (2014);
Nenashev  _et al._ (2015); Baranovskii (2018). This solution is formulated in
Sec. III. Although this simple closed-form system of equations can easily be
solved numerically, it has not yet become a state of the art for the device
simulation community. Single user-friendly equations, like the EGDM, can be
easier implemented into device-simulation software than systems of
interconnected analytical equations, even so the user-friendliness is achieved
at the cost of accuracy.

The challenging task is to replace the system of analytical equations for
$\mu(T,n,F)$ by a single user-friendly appropriately parameterized analytical
equation (APAE) that can be easily embedded into device-simulation software.
Our APAE is formulated in Sec. IV, providing the main result of our paper.

In Sec. V, the validity of the APAE is proven by comparison with the results
of computer simulations. The agreement of this heuristic equation with
analytical theories and with computer simulations suggests the APAE for using
in the device-simulation algorithms.

In Sec. VI, the widely used alternative description of $\mu(T,n,F)$ based on
the EGDM is discussed. In spite of the insufficient parametrization and of the
reduction of simulated systems to regular lattices, the EGDM could be
applicable to some systems with spatial disorder. Sec. VII is dedicated to
comparison with recent studies in the literature, highlighting the advantages
of the APAE.

##  II Crucial parameters for the carrier mobility

###  II.1 Gaussian Disorder Model (GDM)

It has been established that charge transport in single-component and
multicomponent ODSs is due to incoherent hopping of carriers via randomly
distributed localized states with Gaussian energy spectrum Silinsh (1970);
Bässler (1993); Pasveer  _et al._ (2005); Baranovski (2006); van Mensfoort
_et al._ (2008); Germs  _et al._ (2011); Oelerich  _et al._ (2012); Coehoorn
and Bobbert (2012); Kuik  _et al._ (2014); Baranovskii (2014); Nenashev  _et
al._ (2015); Baranovskii (2018); Upreti  _et al._ (2019); Sun  _et al._
(2021); Lee  _et al._ (2021),

$g(\varepsilon)=\frac{N}{\sigma\sqrt{2\pi}}\exp\left(-\frac{\varepsilon^{2}}{2\sigma^{2}}\right).$
(1)

Here, $\sigma$ is the energy scale of the density of states (DOS) and $N$ is
the concentration of randomly distributed localized states, called henceforth
hopping “sites”. The estimates for $\sigma$ vary between $\sigma\simeq 0.05$
eV and $\sigma\simeq 0.15$ eV Bässler (1993); Pasveer  _et al._ (2005); Upreti
_et al._ (2019), the estimates for $N$ vary between $N\simeq 1.7\times
10^{20}$ cm-3 and $N\simeq 4.6\times 10^{21}$ cm-3, depending on the material
Bässler (1993); Pasveer  _et al._ (2005); Lee  _et al._ (2021). For the sake
of simplicity, we follow most previous studies considering the GDM without
correlations between the spatial positions of hopping sites and their energies
Bässler (1993); Pasveer  _et al._ (2005); Baranovski (2006); van Mensfoort
_et al._ (2008); Germs  _et al._ (2011); Oelerich  _et al._ (2012); Coehoorn
and Bobbert (2012); Kuik  _et al._ (2014); Baranovskii (2014); Nenashev  _et
al._ (2015); Baranovskii (2018); Upreti  _et al._ (2019); Sun  _et al._
(2021); Lee  _et al._ (2021).

Phonon-assisted hopping suggested by Miller and Abrahams Miller and Abrahams
(1960) is usually considered as the dominant charge transport mechanism in
ODSs Bässler (1993); Pasveer  _et al._ (2005); Baranovski (2006); van
Mensfoort  _et al._ (2008); Germs  _et al._ (2011); Oelerich  _et al._ (2012);
Coehoorn and Bobbert (2012); Kuik  _et al._ (2014); Baranovskii (2014);
Nenashev  _et al._ (2015); Baranovskii (2018); Upreti  _et al._ (2019); Sun
_et al._ (2021); Lee  _et al._ (2021). The expression for the rate of carrier
transfer from the occupied site with energy $\varepsilon_{i}$ to the empty
site with energy $\varepsilon_{j}$ over the distance $r_{ij}$ has the form

$\nu_{ij}=\nu_{0}\exp\left(-\frac{2|\mathbf{r}_{ij}|}{a}\right)\chi\left(\frac{\varepsilon_{j}-\varepsilon_{i}-e\mathbf{F}\cdot\mathbf{r}_{ij}}{kT}\right)$
(2)

with

$\chi(X)=\begin{cases}\exp(-X),&\text{if $X>0$,}\\\ 1,&\text{if $X\leq
0$}.\end{cases}$

Here, $\mathbf{F}$ is the applied electric field, $e$ is the elementary
charge, $a$ is the localization length of charge carriers in the localized
states, $k$ is the Boltzmann constant, and $T$ is temperature. Estimates of
$a$ in the range $0.1\leq a\leq 0.75$ nm have been suggested in the literature
Gill (1972); Rubel  _et al._ (2004); Pasveer  _et al._ (2005); Lee  _et al._
(2021). The prefactor in Eq. (2) is usually described by a single parameter,
the so-called attempt-to-escape frequency $\nu_{0}$. Precise quantum-
mechanical calculation of $\nu_{0}$ can be found elsewhere Miller and Abrahams
(1960); Shklovskii and Efros (1984). The energy $\varepsilon_{i}$ of the
starting site and the energy $\varepsilon_{j}$ of the target site in Eq. (2)
are counted without contributions of the applied electric field. The effect of
the electric field $\mathbf{F}$ on the hopping rates is expressed explicitly
by the term $e\mathbf{F}\cdot\mathbf{r}_{ij}$ in the exponent of the r.h.s. in
Eq. (2).

The validity of the GDM determined by Eqs. (1) and (2) is justified by its
ability to account for a broad variety of experimental observations. Among
those is the transition from the dependence $\ln[\mu(T)]\propto 1/T^{2}$ to
the dependence $\ln[\mu(T)]\propto 1/T$ with rising carrier
concentrationBaranovskii  _et al._ (2002); Baranovski (2006); Baranovskii
(2014). Another pronounced phenomenon predicted by Eqs. (1) and (2) is the
transition from the mobility $\mu$ independent on carrier concentration $n$ at
small $n$ to the mobility $\mu$ strongly dependent on $n$ at large $n$ values
Baranovskii  _et al._ (2002); Baranovski (2006); Oelerich  _et al._ (2012);
Baranovskii (2014); Nenashev  _et al._ (2015); Baranovskii (2018).
Experimental data to the latter effect Tanase  _et al._ (2003) are indicative
Oelerich  _et al._ (2012); Baranovskii (2014); Nenashev  _et al._ (2015);
Baranovskii (2018) for the Gaussian shape of the DOS given by Eq. (1).
Therefore, our consideration is based on Eqs. (1) and (2), in agreement with
several recent studies Baranovskii (2018); Upreti  _et al._ (2019); Sun  _et
al._ (2021); Lee  _et al._ (2021).

###  II.2 Parametrization of the mobility $\mu(T,n,F)$

For the discusion of the proper parametrization of $\mu(T,n,F)$, we focus on
strong exponential dependencies of the carrier mobility $\mu(T,n,F)$, namely,
on temperature $T$, electric field $F$, and carrier concentration $n$.

####  II.2.1 Parametrization at small electric fields

As evident from the exponents in Eqs. (1), (2), the transport problem at small
electric fields, $F\rightarrow 0$, is determined by only two dimensionless
parameters: $\widetilde{\alpha}=kT/\sigma$ and $\beta=N^{-1/3}/a$. This is
true at low carrier concentrations Baranovskii  _et al._ (2000, 2002); Rubel
_et al._ (2004), $n\ll N$, when the mobility $\mu$ does not depend on $n$. At
elevated $n$, a third dimensionless parameter, $\gamma=n/N$, enters the
carrier mobility $\mu(T,n)$.

Remarkably, the dependence of the carrier mobility $\mu$ on temperature $T$ is
affected by the parameter ($N^{-1/3}/a$).Mott (1969) This is true in the so-
called variable-range hopping (VRH) regime, when the characteristic hopping
length depends on $T$. Charge transport in ODSs is dominated by the VRH
process, as has been proven by Monte Carlo simulations and analytical
calculations Baranovskii  _et al._ (2000, 2002); Rubel  _et al._ (2004);
Nenashev  _et al._ (2017); Oelerich  _et al._ (2017). Therefore, the
temperature dependence $\mu(T,n)$ is sensitive to $\beta$. Baranovskii  _et
al._ (2000, 2002); Rubel  _et al._ (2004); Nenashev  _et al._ (2017); Oelerich
_et al._ (2017) This effect is often overlooked Bässler (1993); Pasveer  _et
al._ (2005).

####  II.2.2 Parametrization of the field dependence

The appropriate parametrization of the dependence $\mu(F)$ in hopping
transport has been revealed in 1973 by Shklovskii Shklovskii (1973), who
concluded that the effect of the electric field $F$ on the carrier mobility
$\mu$ is determined by the product $eaF$, where $a$ is the localization
length. Shklovskii considered for simplicity the case $T=0$, recognizing that
a charge carrier gains the amount of energy $\Delta=eFx$ tunnelling in the
field direction over some distance $x$. The tunneling rate
$\nu(x)\propto\exp(-2x/a)$ can be rewritten as
$\nu(x)\propto\exp(-\Delta/kT_{\mathrm{eff}})$ with $T_{\mathrm{eff}}\simeq
eFa/2$. Apparently, the field-dependent effective temperature
$T_{\mathrm{eff}}\simeq eFa/2$ accounts for the effect of electric field $F$
on hopping transport at $T=0$.

For the case $T\neq 0$, Marianer and Shklovskii Marianer and Shklovskii (1992)
suggested that the combined effects of electric field $F$ and temperature $T$
can be expressed in the form of the effective temperature

$T_{\text{eff}}=T\left[1+\left(c_{1}\frac{eFa}{kT}\right)^{2}\right]^{1/2}$
(3)

with $c_{1}\approx 0.67$. Several studies performed by numerical simulations
Baranovskii  _et al._ (1993); Cleve  _et al._ (1995); Jansson  _et al._
(2008); Nenashev  _et al._ (2017) confirmed the validity of this approach with
values $c_{1}$ distributed in the range $0.5\leq c_{1}\leq 0.9$.

Apparently, the localization length $a$ and not the intersite distance
$N^{-1/3}$ governs the effect of electric field on the hopping conductivity.
This fact is non-trivial because the electric field enters the theory only via
the combination $e\mathbf{F}\cdot\mathbf{r}_{ij}$, in which the length of a
hop $|\mathbf{r}_{ij}|$ is of the order of the intersite distance $N^{-1/3}$.
Therefore, one might expect the combination of parameters $eN^{-1/3}F$ to be
essential for the field-dependent mobility. However, it has been rigorously
proven by straightforward computer simulations Nenashev  _et al._ (2017) that
the localization length $a$, i.e., the feature of a single localized state,
and not the intersite distance $N^{-1/3}$, is responsible for $\mu(F)$. This
counterintuitive result has not yet been adopted by the broad scientific
community in spite of its rigorous proof Nenashev  _et al._ (2017);
Baranovskii (2018). Equation (3) along with Eqs. (1) and (2) implies that the
combined effects of electric field $\mathbf{F}$ and temperature $T$ are
described by a single parameter $\alpha=kT_{\text{eff}}/\sigma$. It means that
the effect of the electric field on the hopping conductivity is governed in
accord with Eq. (3) by the parameter

$\delta=eFa/(kT)\,.$ (4)

Herewith, only three parameters are responsible for $\mu(T,n,F)$, namely,

$\alpha=kT_{\text{eff}}/\sigma;\ \beta=N^{-1/3}/a;\ \gamma=n/N\ .$ (5)

##  III Analytical description of $\mu(T,n,F)$

The analytical theory for the description of hopping transport in amorphous
materials with strongly energy-dependent DOS $g(\varepsilon)$ is known for
decades Shklovskii  _et al._ (1990); Baranovskii  _et al._ (2005); Baranovski
(2006). In particular, it has been proven Oelerich  _et al._ (2012);
Baranovskii (2014); Nenashev  _et al._ (2015); Baranovskii (2018) that the
charge carrier mobility can be desribed in the framework of the GDM as

$\displaystyle\begin{split}\mu=\mu_{0}\gamma^{-1}\exp\left(-\frac{2B_{c}^{1/3}}{a}r(\varepsilon_{t})-\frac{\varepsilon_{t}-\varepsilon_{F}}{kT}\right)\,,\end{split}$
(6)

where

$\displaystyle
r(\varepsilon_{t})=\left[\frac{4\pi}{3}\int\limits_{-\infty}^{\varepsilon_{t}}g(\varepsilon^{\prime})[1-f(\varepsilon^{\prime},\varepsilon_{F})]d\varepsilon^{\prime}\right]^{-1/3}\,,$
(7)

and the transport energy, $\varepsilon_{t}$, is calculated from Oelerich  _et
al._ (2012); Baranovskii (2014, 2018)

$\displaystyle\begin{split}\frac{2}{3}\left(\frac{4\pi}{3B_{c}}\right)^{-\frac{1}{3}}\frac{kT}{a}&\left[\int_{-\infty}^{\varepsilon_{t}}[1-f(\varepsilon,\varepsilon_{F})]g(\varepsilon)d\varepsilon\right]^{-\frac{4}{3}}\\\
&\times[1-f(\varepsilon_{t},\varepsilon_{F})]g(\varepsilon_{t})=1\,,\end{split}$
(8)

where $f(\varepsilon,\varepsilon_{F})$ is Fermi function

$\displaystyle
f(\varepsilon,\varepsilon_{F})=\left[1+\exp{\frac{(\varepsilon-\varepsilon_{F})}{kT}}\right]^{-1}\,$
(9)

with the Fermi energy $\varepsilon_{F}$ determined by the relation

$\displaystyle\int_{-\infty}^{\infty}g(\varepsilon)f(\varepsilon,\varepsilon_{F})d\varepsilon=n\,,$
(10)

which accounts for the finite charge carrier concentration $n$. The
coefficient $B_{c}\simeq 2.7$ is due to the percolation nature of the hopping
transport Rubel  _et al._ (2004).

The preexponential factor in Eq. (6) has the value

$\displaystyle\mu_{0}=B\frac{e\nu_{0}}{kTN^{2/3}}\,,$ (11)

where $B$ is a numerical factor, which can be determined by comparison with
computer simulations.

The structure of prefactor $\mu_{0}$ in Eq. (11) relies on the conventional
form of the Einstein relation between the carrier mobility and the carrier
diffusion coefficient that is valid for non-degenerate systems, i.e., at low
carrier concentrations $n\ll N$. In general, the ratio $e/kT$ should be
replaced by the generalized Einstein relation Roichman and Tessler (2002);
Tessler  _et al._ (2009); Baranovskii (2014). We nevertheless leave Eq. (11)
in the given form because the final result for $\mu(T,n,F)$ in the next Sec.
IV is obtained from a comparison of the parameters in the analytical theory
with numerical data from computer simulations. In this respect, our approach
is essentially similar to that exploited in the recent studies by Upreti et
al. Upreti  _et al._ (2019) and by Lee et al. Lee  _et al._ (2021), who
calibrated parameters of a similar analytical theory by comparison with
numerical simulations.

##  IV Appropriately parameterized analytical equation (APAE) for $\mu(T,n,F)$

The set of Eqs. (6)-(11) determines $\mu(T,n,F)$ in terms of the relevant
parameters given by Eq. (5). Remarkably, the function $\mu(T,n,F)$ obtained
from this set of equations can be approximated in the realistic parameter
ranges by a single appropriately parameterized analytical equation (APAE)

$\displaystyle\begin{split}\mu(T,n,F)=\frac{e\nu_{0}}{\sigma
N^{2/3}}\Phi(\gamma)\exp\left[-A(\beta)-\alpha^{-2}C(\beta)\right]\,,\end{split}$
(12)

where

$\displaystyle A(\beta)=c_{2}\beta^{2}+c_{3}\beta+c_{4}\,,$ (13)
$\displaystyle C(\beta)=c_{5}\beta^{-2}+c_{6}\beta^{-1}+c_{7}\,.$ (14)

The constant $\exp(c_{4})$ replaces the factor $B$ in Eq. (11).

The function $\Phi(\gamma)$ in Eq. (12), responsible for the dependence
$\mu(n)$, is determined by the equation

$\displaystyle\Phi(\gamma)=\exp\left[(c_{8}\ln\alpha^{-1}-c_{9})\xi^{2}\Theta(\xi)\right]\,,$
(15)

where

$\displaystyle\xi=\frac{\alpha^{-1}}{2}+2-\left|\ln\left(\gamma-\frac{c_{10}\gamma^{2}}{\max(\beta,c_{11})}\right)\right|^{1/2}\,,$
(16) $\Theta(x)=\begin{cases}1,&\text{if $x\geq 0$,}\\\ 0,&\text{if
$x<0$}.\end{cases}$

Numerical parameters $c_{1}$ to $c_{11}$ were further optimized by comparison
with computer simulations attaining the values

$\begin{split}c_{1}&=0.6\\\ c_{2}&=-0.066\\\ c_{3}&=2.65\\\ c_{4}&=-1.35\\\
c_{5}&=0.89\\\ c_{6}&=-0.86\\\ c_{7}&=0.54\\\ c_{8}&=0.6\\\ c_{9}&=0.15\\\
c_{10}&=13\\\ c_{11}&=3\end{split}$ (17)

Recall that $\alpha=kT_{\text{eff}}/\sigma,\beta=N^{-1/3}/a,\gamma=n/N$, see
Eq. (5), and $T_{\text{eff}}$ is determined by Eq. (3) with parameter $c_{1}$
given by Eq. (17). Equation (12) is the central result of our work.

##  V Testing the APAE by computer simulations

###  V.1 APAE and $\mu(F)$

In Fig. 1, the dependences of the normalized mobility $\mu(F)/\mu(0)$ obtained
by Monte Carlo (MC) simulations at different values of the parameter
$\beta=N^{-1/3}/a$ are depicted by symbols. Simulations were carried out in
the framework of the model with spatial and energy disorder formulated in Sec.
II.1. The simulation algorithm is the same as the one used in previous studies
Nenashev  _et al._ (2017). In particular, we used the HopHop algorithm
(https://github.com/janoliver/hophop). Simulation data at $\beta^{-1}\leq 0.3$
agree with those from the previous studies Nenashev  _et al._ (2017), while
the data at $\beta=2$ are novel. At all values of $\beta$ in Fig. 1,
$\mu(F)/\mu(0)$ steeply increases with rising field $F$ at $F\leq
2\sigma/eN^{-1/3}$. At larger fields, the field dependences in Fig. 1
saturate. Solid lines in Fig. 1 depict the results of the APAE. Apparently,
APAE appropriately describes the simulation data at $F\leq 2\sigma/eN^{-1/3}$,
while at larger fields, APAE predicts a stronger dependence $\mu(F)/\mu(0)$
than the one yielded by simulations.

The weakening of the field dependency at high fields in Fig. 1 is not
surprising. There are several mechanisms contributing to this effect
Baranovskii (2014, 2018). In fact, $\mu(F)$ in hopping transport should even
decrease at large electric fields Nguyen Van Lien and Shklovskii (1981); Levin
and Shklovskii (1988); Nenashev  _et al._ (2008). Such effects at very high
electric fields cannot be described by the effective temperature in Eq. (3).

It is seen in Fig. 1 that APAE based on Eq. (3) describes the dependence
$\mu(F)$ only at electric fields smaller than $F^{*}\simeq 2\sigma/eN^{-1/3}$.
Taking relevant parameter values for ODSs Pasveer  _et al._ (2005); Lee  _et
al._ (2021), $N\approx 10^{21}$ cm-3 and $\sigma\approx 0.1$ eV, one obtains
$F^{*}\approx 2\times 10^{6}$ Vcm-1. The range of electric fields $F<2\times
10^{6}$ Vcm-1 is relevant to experimental studies Hirao  _et al._ (1995);
Mozer  _et al._ (2005) and to device applications Upreti  _et al._ (2019).
Hence, the APAE can be used to interpret experimental data and it can be
applied in device simulation algorithms.

Figure 1: Normalized mobility $\mu(F)/\mu(0)$ at $\sigma/kT=4$ and
$n\rightarrow 0$ as a function of $eFN^{-1/3}/\sigma$ at different values of
$\beta=N^{-1/3}/a$. Symbols: results of MC simulations; solid lines: results
of the APAE equation (12).

###  V.2 APAE and $\mu(T)$

The scale of energy disorder $\sigma$ is usually estimated in ODSs between
$\sigma\simeq 50$ meV and $\sigma\simeq 150$ meV Bässler (1993); Upreti  _et
al._ (2019); Lee  _et al._ (2021). This gives for the ratio $\sigma/(kT)$ at
room temperature $T\thickapprox 300$ K the values between $\sigma/(kT)\simeq
2$ and $\sigma/(kT)\simeq 6$.

In Fig. 2, the mobility $\mu(F\rightarrow 0,n\rightarrow 0)$, as obtained by
computer simulations, is plotted as a function of $\sigma/(kT$) at different
values of $\beta=N^{-1/3}/a$. The results of the APAE are in perfect agreement
with the simulation data. In contrast to the MC simulations used for the data
in Fig. 1, numerical calculations at $F\rightarrow 0$ were performed by
solution of Kirchhoff equations in a resistor network proposed by Miller and
Abrahams Shklovskii and Efros (1984).

Figure 2: Mobility $\mu(T)$ at $n\rightarrow 0$, $F\rightarrow 0$ as a
function of $\alpha^{-2}=\sigma^{2}/(kT)^{2}$ for different values of
$\beta=N^{-1/3}/a$. Symbols: results of simulations; solid lines: results of
the APAE.

The data in Fig. 2 correspond to the well-known dependence $\mu(T)$, which can
be approximated by the expression

$\mu(T)\propto\exp\left[-C\left(\frac{\sigma}{kT}\right)^{2}\right].$ (18)

By fitting the Monte Carlo simulation data for $\mu(T)$ obtained on a cubic
lattice at $\beta=N^{-1/3}/a=5$, Bässler Bässler (1993) suggested the value
$C=(2/3)^{2}\simeq 0.44$. Using analytical theory based on the approach of
transport energy given by Eq. (8), it was shown Baranovskii  _et al._ (2000)
that in the GDM with spatial and energy disorder described in Sec. II.1,
coefficient $C$ depends on parameter $\beta=N^{-1/3}/a$ having the values
$C\simeq 0.46$ at $\beta=10$, $C\simeq 0.41$ at $\beta=5$ and $C\simeq 0.38$
at $\beta\simeq 3.7$. These values have been confirmed by analytical
calculations based on the percolation theory Baranovskii  _et al._ (2002). The
data depicted in Fig. 2 are in good agreement with these previous results
yielding the values $C$ specified in the figure.

###  V.3 APAE and $\mu(n)$

In Fig. 3, the normalized mobility $\mu(n)/\mu(n\rightarrow 0)$ at
$F\rightarrow 0$ is plotted as a function of $\gamma=n/N$ at a typical ODSs
parameter Bässler (1993) $\beta=5$ and different values of $\alpha=kT/\sigma$.
The choice of $\beta\simeq 5$ is further supported by the estimates Rubel  _et
al._ (2004) $0.1$ nm $\leq a\leq 0.3$ nm and $0.5$ nm $\leq N^{-1/3}\leq 1.8$
nm in ODSs Rubel  _et al._ (2004); Pasveer  _et al._ (2005).

Figure 3: Normalized mobility $\mu(n)/\mu(n\rightarrow 0)$ at $F\rightarrow 0$
as a function of $\gamma=n/N$ at $\beta=5$ and different values of
$\alpha=kT/\sigma$. Symbols: results of simulations; solid lines: results of
the APAE.

Apparently, the computer simulations confirm the capability of the APAE to
describe the dependences $\mu(T,n,F)$ in the framework of the GDM, i.e., for
hopping transport via randonly distributed localized states with Gaussian
energy spectrum.

Most algorithms developed so far for the simulation of devices based on ODSs,
such as organic light–emitting diodes, organic solar cells, and organic field-
effect transistors, are based on the EGDM equation Pasveer  _et al._ (2005);
Coehoorn and Bobbert (2012). In the following sections, we analyze the EGDM
equation and compare it with APAE.

##  VI Comparison with EGDM

A reduced version of the GDM without spatial disorder is, in fact, what one
calls EGDM Coehoorn and Bobbert (2012). Charge transport on a regular cubic
lattice was simulated for the case $b=10a$, where $b$ is the lattice spacing
and $a$ is the localization length in Eq. (2). To make a bridge to the case of
randomly placed localized states, we note that $b=N^{-1/3}$. Simulation
results for the carrier mobility $\mu(T,n,F)$ were fitted in the form Pasveer
_et al._ (2005)

$\displaystyle\mu(T,n,F)\approx\widetilde{\mu}(T,n)\phi(T,F)\,,$ (19)

where $\phi(T,F)$ is

$\displaystyle\begin{split}\phi(T,F)=\exp\Biggl{\\{}&0.44\left[\widehat{\sigma}^{3/2}-2.2\right]\\\
&\times\left[\sqrt{1+0.8\left(\frac{Feb}{\sigma}\right)^{2}}-1\right]\Biggr{\\}}\,\end{split}$
(20)

with $\widehat{\sigma}=\sigma/(kT)$. The function $\widetilde{\mu}(T,n)$ is
given by Pasveer  _et al._ (2005)

$\displaystyle\widetilde{\mu}(T,n)=\mu_{0}(T)\exp[(\widehat{\sigma}^{2}-\widehat{\sigma})(2nb^{3})^{\widetilde{\delta}}/2],$
(21a)
$\displaystyle\mu_{0}(T)=\mu_{0}\mathfrak{c}_{1}\exp[-\mathfrak{c}_{2}\widehat{\sigma}^{2}],$
(21b) $\displaystyle\widetilde{\delta}\equiv
2\frac{\ln(\widehat{\sigma}^{2}-\widehat{\sigma})-\ln(\ln
4)}{\widehat{\sigma}^{2}},\,\mu_{0}\equiv\frac{b^{2}\nu_{0}e}{\sigma},$ (21c)

with $\mathfrak{c}_{1}=1.8\times 10^{-9}$ and $\mathfrak{c}_{2}=0.42$.
Equations (19) – (21c) were named Coehoorn and Bobbert (2012) the extended
Gaussian disorder model (EGDM). These equations fit numerical data for
$\mu(T,n,F)$ obtained on a cubic lattice at $b/a=10$.

Remarkably, the localization length of charge carriers, $a$, does not enter
the EGDM equations at all, although $a$ is known Mott (1969); Shklovskii
(1973) to determine the temperature- and the field-dependencies of the hopping
mobility via parameters $\alpha=kT_{\text{eff}}/\sigma$, $\beta=N^{-1/3}/a$
present in Eqs. (3), (5).

The absence of $a$ in the EGDM, was supported by two arguments. First, it was
stated that presumably in all published 3D-modeling work on the (E)GDM a fixed
value $N^{-1/3}/a=10$ had been used, suggesting $N^{-1/3}/a=10$ as a standard
one for ODSs Coehoorn and Bobbert (2012). In this respect, one has to remark,
however, that published 3D-modeling work on the GDM had mostly used a fixed
value $N^{-1/3}/a=5$ instead Bässler (1993).

Second, it was stated, referring to unpublished work, that varying the ratio
$N^{-1/3}/a$ has no significant effect on the temperature, field, or carrier
density dependence of the mobility Coehoorn and Bobbert (2012). Fig. 1
evidences, however, that varying the ratio $N^{-1/3}/a$ does have a
significant effect on the field dependence of the mobility in the GDM with
spatial disorder. A similar effect on cubic lattices has been proven elsewhere
Nenashev  _et al._ (2017); Baranovskii (2018).

This analysis raises a question of whether EGDM can be used for description of
$\mu(T,n,F)$ in ODSs. The answer to this question is of high importance not
only for academic researchers, but particularly for the community dealing with
device simulations, where EGDM is, so far, considered as the state of the art.
It is a lucky coincidence that the EGDM equation could be, in some cases,
applicable to ODSs.

In Fig. 4, we depict by symbols the results of Monte Carlo simulations
Nenashev  _et al._ (2017); Oelerich  _et al._ (2017) for the system of random
sites at $0.1\leq\beta^{-1}\leq 0.3$ along with the result of EGDM (dashed
line), which fits Monte Carlo simulations Pasveer  _et al._ (2005) on a cubic
lattice at $\beta^{-1}=0.1$. The case $\sigma/kT=4$ is considered as typical
for ODSs at room temperature, since $\sigma$ is usually estimated as
$\sigma\simeq 0.1$ eV Pasveer  _et al._ (2005); Lee  _et al._ (2021).

Apparently, $\mu(F)$ on the lattice unequals $\mu(F)$ in a system of random
sites. Equation (20) fits the simulation data on a cubic lattice at
$\beta^{-1}=0.1$, though it does not fit the simulation data in a spatially
disordered system at $\beta^{-1}=0.1$ (black circles). The physical mechanism
responsible for the drastic difference in $\mu(F)$ between lattices and random
sites has been discussed in detail elsewhere Oelerich  _et al._ (2017).

From Fig. 4 we, however, learn that the EGDM, though desired to fit $\mu(F)$
on a cubic lattice at $\beta^{-1}=0.1$, occasionally fits well $\mu(F)$ on
random sites at $\beta^{-1}=0.18$. While the value $\beta=N^{-1/3}/a=10$ was
standardised in the EGDM Coehoorn and Bobbert (2012), the previously used
value $\beta=N^{-1/3}/a=5$ looks more relevant to ODSs Bässler (1993). This
conclusion is supported by the estimates for $a$ Gill (1972); Rubel  _et al._
(2004); Pasveer  _et al._ (2005); Lee  _et al._ (2021), $0.1$ nm $\leq a\leq
0.75$ nm and estimates for $N$ between $N\simeq 1.7\times 10^{20}$ cm-3 and
$N\simeq 4.6\times 10^{21}$ cm-3 depending on the material Bässler (1993);
Pasveer  _et al._ (2005); Lee  _et al._ (2021). The data in Fig. 4 suggest,
therefore, that EGDM could be applied to ODSs with a realistic value
$\beta^{-1}=a/N^{-1/3}\simeq 0.18$. For materials with $\beta^{-1}>0.18$, one
should, instead, use the APAE for $\mu(F)$, as illustrated in Fig. 4.

Figure 4: Normalized mobility $\mu(F)/\mu(0)$ at $n\rightarrow 0$ for
$\sigma/kT=4$ as a function of $\delta\beta$ at different values of $\beta$.
Symbols: Monte Carlo simulation data obtained for random cites at
$\beta^{-1}=0.1$ (black), $\beta^{-1}=0.18$ (blue), $\beta^{-1}=0.22$ (green),
and $\beta^{-1}=0.3$ (red). The orange dashed line is the result of EGDM,
i.e., the data for a cubic lattice at $\beta^{-1}=0.1$.

Furthermore, EGDM occasionally appears capable to account, in some cases, for
the dependece $\mu(T)$ on a system with spatial disorder. In fact, the
localization length $a$ affects not only the field dependence $\mu(F)$, but
also the temperature dependence $\mu(T)$ Baranovskii  _et al._ (2000, 2002);
Lee  _et al._ (2021); Sun  _et al._ (2021). This effect is neglected in the
EGDM. The effect of $a$ on the dependence $\mu(T)$ can be taken into account,
for instance,Baranovskii  _et al._ (2000, 2002); Lee  _et al._ (2021); Sun
_et al._ (2021) by replacing the constant $\mathfrak{c}_{2}$ in Eq. (21b) by
the appropriate functions of the parameter $\beta=N^{-1/3}/a$. The
illustration is provided by Fig. 2, where the slopes of the straight lines
depend on $\beta=N^{-1/3}/a$.

At low carrier concentrations and low electric fields, EGDM describes the
dependence $\mu(T)$ by Eq. (18) with $C=0.42$. As evident in Fig. 2, the
dependence $\mu(T)$ for carrier mobility on random sites with $\beta\simeq 5$,
typical for ODSs Bässler (1993), is described by Eq. (18) with $C\simeq 0.40$,
which is very close to $C\simeq 0.42$ given by the EGDM. At $\beta$ values
different to $\beta\simeq 5$, the APAE should be preferred.

##  VII Discussion

###  VII.1 Comparison with recent results in the literature

Several recent studies Upreti  _et al._ (2019); Lee  _et al._ (2021); Sun  _et
al._ (2021) were dedicated to the appropriate theoretical description of the
dependences $\mu(T,n,F)$ suitable for incorporation into device simulation
algorithms.

Upreti et al. Upreti  _et al._ (2019) considered an analytical theory for
$\mu(T,n,F)$ similar to that described in Sec. III and carried out Monte Carlo
simulations used to calibrate parameters in the analytical theory. Their
theory was further implemented into the drift-diffusion solver. Along with the
theoretical development, Upreti et al. Upreti  _et al._ (2019) fabricated
hole-only and electron-only organic devices and performed measurements of
temperature-dependent space-charge-limited currents. It was found that the
suggested theory can adequately describe electron and hole transport in a wide
variety of organic semiconductor blends, which are used as the active layer in
typical bulk heterojunction organic solar cells Deibel and Dyakonov (2010). It
was also recognised that EGDM fails to produce an acceptable fit for
experimental data obtained on the electron-only devices Upreti  _et al._
(2019). Upreti et al. supposed that the reason for the failure of the EGDM for
electron-only devices is a small value of parameter $\beta=N^{-1/3}/a=2$,
which favors the variable-range hopping (VRH), while EGDM with the large fixed
parameter $\beta=10$ presumably favors the nearest-neighbor hopping mechanism.
However, Monte Carlo simulations have proven Oelerich  _et al._ (2017) that
EGDM is, in fact, based on the VRH transport regime, though on a cubic lattice
without spatial disorder. Therefore, more study would be needed to clarify why
the EGDM fails to account for experimental data in electron-only devices,
though being capable to describe the data in hole-only devices reported by
Upreti et al. Upreti  _et al._ (2019).

Lee et al. Lee  _et al._ (2021) also addressed an analytical theory for
$\mu(T,n,F)$ described in Sec. III. They implemented the results in a
technology computer-aided design (TCAD) simulation tool, ATLAS, from Silvaco
(Silvaco Inc., Atlas Ver. 5.30.0.R (2020)) and deduced parameters of the GDM
by comparison with experimental data for ODSs with high carrier mobilities.
Only the case of low electric fields was considered by Lee et al. Lee  _et
al._ (2021), who highlighted that the GDM with spatial and energy disorder is
superior to the EGDM that lacks spatial disorder and uses only one fixed value
$\beta=N^{-1/3}/a=10$ on a simple cubic lattice.

Neither Upreti et al. Upreti  _et al._ (2019), nor Lee et al. Lee  _et al._
(2021) attempted to replace the set of analytical equations by a single
closed-form expression for $\mu(T,n,F)$ as, for example, APAE derived in Sec.
IV.

Sun et al. Sun  _et al._ (2021) improved the EGDM, though not going beyond the
model based on a cubic lattice. One improvement is the replacement of
constants $\mathfrak{c}_{1}$ and $\mathfrak{c}_{2}$ in Eq. (21b) by functions
of the parameter $\beta=N^{-1/3}/a$. The other improvement is the replacement
of the lattice constant $b$ in Eq. (20) by $b=10a$ for using the localization
length $a$ as adjustable parameter Sun  _et al._ (2021). Sun et al. Sun  _et
al._ (2021) solved analytically the degenerate drift-diffusion equation and
extracted model parameters for several organic diodes by comparison with
experimental data.

Figure 5: Proof that $\alpha=kT/\sigma$ and $\beta=N^{-1/3}/a$ are the cricial
parameters of the GDM and not separately $kT$, $\sigma$, $N^{-1/3}$, $a$. (a)
parameter $\alpha=kT/\sigma$ has the same value $\alpha\simeq 0.32$ for the
both sets of data; (b) parameter $\beta=N^{-1/3}/a$ has the same value
$\beta=5$ for the both sets of data.

These studies have not revealed the combinations of parameters given in Eq.
(5) as parametrization of the theoretical model. Instead, parameters
$N^{-1/3}$ and $a$ and parameters $\sigma$ and $kT$ were treated separately
from each other. In order to emphasize once again the validity of the
parametrization given in Sec. II.2, we use in Fig. 5 the data of Upreti et al.
Upreti  _et al._ (2019) obtained by Monte Carlo simulations for the dependence
$\mu(n)$ at different values of parameters $N^{-1/3}$, $a$, $\sigma$, and
$kT$. While Upreti et al. Upreti  _et al._ (2019) plotted their numerical data
for $\mu(n)$ as functions of $N^{-1/3}$ at fixed $a$ and as functions of $a$
at fixed $N^{-1/3}$, as well as functions of $\sigma$ at fixed $T$ and
functions of $T$ at fixed $\sigma$, we have chosen their combinations of
$N^{-1/3}$ and $a$, as well as the combinations of $T$ and $\sigma$, which
correspond to the same values of the ratios $\alpha=kT/\sigma$ and
$\beta=N^{-1/3}/a$. The plots illustrate that, in fact, $\alpha=kT/\sigma$ and
$\beta=N^{-1/3}/a$ control $\mu(n)$ and not separately $N^{-1/3}$, $a$,
$\sigma$, and $kT$. This is the reason why the APAE introduced in Sec. IV is
formulated in terms of only three parameters $\alpha$, $\beta$, $\gamma$
determined in Eq. (5).

###  VII.2 Conclusions

The APAE formulated for $\mu(T,n,F)$ in Sec. IV is well-justified by classical
theories of hopping transport. It is calibrated by Monte Carlo simulations in
the framework of the GDM. It is user-friendly and can be easily implemented
into device-simulation software.

EGDM Pasveer  _et al._ (2005) desired to fit simulation data for $\mu(T,n,F)$
in the framework of the GDM on a cubic lattice at fixed value
$\beta=N^{-1/3}/a=10$ occasionally agrees with simulation data on a system of
random sites, though at $\beta\simeq 5$. For other $\beta$ values than
$\beta\simeq 5$, APAE should be preferred.

###### Acknowledgements.

A.N. thanks the Faculty of Physics of the Philipps University Marburg for the
kind hospitality during his research stay. S.D.B. and K.M. acknowledge
financial support by the Deutsche Forschungsgemeinschaft (Research Training
Group “TIDE”, RTG2591).

## References

  * Upreti  _et al._ (2019) T. Upreti, Y. Wang, H. Zhang, D. Scheunemann, F. Gao, and M. Kemerink, Phys. Rev. Applied 12, 064039 (2019).
  * Lee  _et al._ (2021) Y. Lee, S. Jung, A. Plews, A. Nejim, O. Simonetti, L. Giraudet, S. D. Baranovskii, F. Gebhard, K. Meerholz, S. Jung, G. Horowitz, and Y. Bonnassieux, Phys. Rev. Applied 15, 024021 (2021).
  * Sun  _et al._ (2021) J.-X. Sun, H.-C. Yang, Y. Li, and H.-J. Cui, Phys. Rev. Applied 16, 034037 (2021).
  * Oelerich  _et al._ (2012) J. O. Oelerich, D. Huemmer, and S. D. Baranovskii, Phys. Rev. Lett. 108, 226403 (2012).
  * Nenashev  _et al._ (2017) A. V. Nenashev, J. O. Oelerich, A. V. Dvurechenskii, F. Gebhard, and S. D. Baranovskii, Phys. Rev. B 96, 035204 (2017).
  * Oelerich  _et al._ (2017) J. O. Oelerich, A. V. Nenashev, A. V. Dvurechenskii, F. Gebhard, and S. D. Baranovskii, Phys. Rev. B 96, 195208 (2017).
  * Pasveer  _et al._ (2005) W. F. Pasveer, J. Cottaar, C. Tanase, R. Coehoorn, P. A. Bobbert, P. W. M. Blom, D. M. de Leeuw, and M. A. J. Michels, Phys. Rev. Lett. 94, 206601 (2005).
  * Coehoorn and Bobbert (2012) R. Coehoorn and P. A. Bobbert, Phys. Status Solidi A 209, 2354 (2012).
  * Kuik  _et al._ (2014) M. Kuik, G.-J. A. H. Wetzelaer, H. T. Nicolai, N. I. Craciun, D. M. De Leeuw, and P. W. M. Blom, Adv. Mater. 26, 512 (2014).
  * Mott (1969) N. F. Mott, Phil. Mag. 19, 835 (1969).
  * Shklovskii (1973) B. I. Shklovskii, Sov. Phys. Semicond. 6, 1964 (1973).
  * Baranovskii (2018) S. D. Baranovskii, Phys. Status Solidi A 215, 1700676 (2018).
  * Bässler (1993) H. Bässler, Phys. Status Solidi B 175, 15 (1993).
  * Schmechel (2002) R. Schmechel, Phys. Rev. B 66, 235206 (2002).
  * Tessler  _et al._ (2009) N. Tessler, Y. Preezant, N. Rappaport, and Y. Roichman, Adv. Mater. 21, 2741 (2009).
  * Baranovskii  _et al._ (2000) S. D. Baranovskii, H. Cordes, F. Hensel, and G. Leising, Phys. Rev. B 62, 7934 (2000).
  * Baranovskii  _et al._ (2002) S. D. Baranovskii, I. P. Zvyagin, H. Cordes, S. Yamasaki, and P. Thomas, Phys. Status Solidi B 230, 281 (2002).
  * Rubel  _et al._ (2004) O. Rubel, S. D. Baranovskii, P. Thomas, and S. Yamasaki, Phys. Rev. B 69, 014206 (2004).
  * Baranovski (2006) S. Baranovski, ed.,  _Charge Transport in Disordered Solids with Applications in Electronics_ (John Wiley & Sons, Ltd, Chichester, 2006).
  * Baranovskii (2014) S. D. Baranovskii, Phys. Status Solidi B 251, 487 (2014).
  * Nenashev  _et al._ (2015) A. V. Nenashev, J. O. Oelerich, and S. D. Baranovskii, J. Phys.: Condens. Matter 27, 093201 (2015).
  * Silinsh (1970) E. A. Silinsh, Phys. Status Solidi A 8, 817 (1970).
  * van Mensfoort  _et al._ (2008) S. L. M. van Mensfoort, S. I. E. Vulto, R. A. J. Janssen, and R. Coehoorn, Phys. Rev. B 78, 085208 (2008).
  * Germs  _et al._ (2011) W. C. Germs, J. J. M. van der Holst, S. L. M. van Mensfoort, P. A. Bobbert, and R. Coehoorn, Phys. Rev. B 84, 165210 (2011).
  * Miller and Abrahams (1960) A. Miller and E. Abrahams, Phys. Rev. 120, 745 (1960).
  * Gill (1972) W. Gill, J. Appl. Phys. 43, 5033 (1972).
  * Shklovskii and Efros (1984) B. I. Shklovskii and A. L. Efros,  _Electronic Properties of Doped Semiconductors_ (Springer, Berlin, 1984).
  * Tanase  _et al._ (2003) C. Tanase, E. Meijer, P. Blom, and D. de Leeuw, Phys. Rev. Lett. 91, 216601 (2003).
  * Marianer and Shklovskii (1992) S. Marianer and B. I. Shklovskii, Phys. Rev. B 46, 13100 (1992).
  * Baranovskii  _et al._ (1993) S. D. Baranovskii, B. Cleve, R. Hess, and P. Thomas, J. Non-Cryst. Solids 437–440, 164 (1993).
  * Cleve  _et al._ (1995) B. Cleve, B. Hartenstein, S. D. Baranovskii, M. Scheidler, P. Thomas, and H. Bässler, Phys. Rev. B 51, 16705 (1995).
  * Jansson  _et al._ (2008) F. Jansson, S. D. Baranovskii, F. Gebhard, and R. Österbacka, Phys. Rev. B 77, 195211 (2008).
  * Shklovskii  _et al._ (1990) B. I. Shklovskii, E. I. Levin, H. Fritzsche, and S. D. Baranovskii, Advances in disordered semiconductors, vol. 3, ed. fritzsche, h. (World Scientific, Singapore, 1990) Chap. Hopping photoconductivity in amorphous semiconductors: dependence on temperature, electric field and frequency, pp. 161–191.
  * Baranovskii  _et al._ (2005) S. D. Baranovskii, O. Rubel, and P. Thomas, Thin Solid Films 487, 2 (2005).
  * Roichman and Tessler (2002) Y. Roichman and N. Tessler, Appl. Phys. Lett. 80, 1948 (2002).
  * Nguyen Van Lien and Shklovskii (1981) Nguyen Van Lien and B. I. Shklovskii, Solid State Commun. 38, 99 (1981).
  * Levin and Shklovskii (1988) E. I. Levin and B. I. Shklovskii, Solid State Commun. 67, 233 (1988).
  * Nenashev  _et al._ (2008) A. V. Nenashev, F. Jansson, S. D. Baranovskii, R. Österbacka, A. V. Dvurechenskii, and F. Gebhard, Phys. Rev. B 78, 165207 (2008).
  * Hirao  _et al._ (1995) A. Hirao, H. Nishizawa, and M. Sugiuchi, Phys. Rev. Lett. 75, 1787 (1995).
  * Mozer  _et al._ (2005) A. J. Mozer, N. S. Sariciftci, A. Pivrikas, R. Österbacka, G. Juška, L. Brassat, and H. Bässler, Phys. Rev. B 71, 035214 (2005).
  * Deibel and Dyakonov (2010) C. Deibel and V. Dyakonov, Reports on Progress in Physics 73, 096401 (2010).


# Global dynamics of epidemic network models via construction of Lyapunov
functions

Rachidi B<EMAIL_ADDRESS>Department of Mathematical
Sciences, University of Nevada Las Vegas, Las Vegas, USA and Yixiang
<EMAIL_ADDRESS>Department of Mathematical Sciences, Middle
Tennessee State University, Murfreesboro, Tennessee 37132, USA

###### Abstract

In this paper, we study the global dynamics of epidemic network models with
standard incidence or mass-action transmission mechanism, when the dispersal
of either the susceptible or the infected people is controlled. The
connectivity matrix of the model is not assumed to be symmetric. Our main
technique to study the global dynamics is to construct novel Lyapunov type
functions.

##  1 Introduction

In this paper, we study the global dynamics of the following Susceptible-
Infected-Susceptible (SIS) epidemic network (patch) model

$\left\\{\begin{array}[]{lll}\displaystyle
S_{i}^{\prime}&=d_{S}\displaystyle\sum_{j\in\Omega,j\neq
i}(L_{ij}S_{j}-L_{ji}S_{i})-f_{i}(S_{i},I_{i})+\gamma_{i}I_{i},&i\in\Omega,\
t>0,\cr\displaystyle I_{i}^{\prime}&=d_{I}\displaystyle\sum_{j\in\Omega,j\neq
i}(L_{ij}I_{j}-L_{ji}I_{i})+f_{i}(S_{i},I_{i})-\gamma_{i}I_{i},&i\in\Omega,\
t>0,\end{array}\right.$ (1.1)

where $\bm{S}=(S_{1},\dots,S_{n})$ and $\bm{I}=(I_{1},\dots,I_{n})$ denote the
number of susceptible and infected people living in $n-$patches (e.g., cities,
countries, etc.) respectively; $n\geq 2$ is the number of patches;
$\Omega=\\{1,2,\cdots,n\\}$ is the collection of patch labels; $\gamma_{i}$ is
the disease recovery rate; and $f_{i}(S_{i},I_{i})$ describes the interaction
of susceptible and infected people. The coefficient $L_{ij}$, $j\neq i$, is
the degree of movement of individuals from patch $j$ to $i$, and $d_{S}$ and
$d_{I}$ are the dispersal rates of susceptible and infected people,
respectively. For simplicity, let $L_{ii}=-\sum_{j\in\Omega,j\neq i}L_{ji}$,
which is the total degree of movement out from patch $i$. The matrix
$\mathcal{L}=(L_{ij})$ is called the connectivity matrix, which will not be
assumed to be symmetric. We can rewrite model (1.1) as

$\left\\{\begin{array}[]{lll}\displaystyle
S_{i}^{\prime}&=d_{S}\displaystyle\sum_{j\in\Omega}L_{ij}S_{j}-f_{i}(S_{i},I_{i})+\gamma_{i}I_{i},&i\in\Omega,\
t>0,\cr\displaystyle
I_{i}^{\prime}&=d_{I}\displaystyle\sum_{j\in\Omega}L_{ij}I_{j}+f_{i}(S_{i},I_{i})-\gamma_{i}I_{i},&i\in\Omega,\
t>0.\end{array}\right.$ (1.2)

Adding up all the equations in (1.1), we obtain
$\frac{d}{dt}\sum_{i\in\Omega}(S_{i}+I_{i})=0$. This means that the total
population $N:=\sum_{i\in\Omega}(S_{i}+I_{i})$ remains a constant for all
$t\geq 0$.

Model (1.1) when $d_{S}>0$ and $d_{I}>0$ has been studied recently in order to
understand the joint impact of environmental heterogeneity and population
movement on disease transmission. In [1], Allen et al. studied (1.1) with
standard incidence transmission mechanism
$f_{i}=\beta_{i}S_{i}I_{i}/(S_{i}+I_{i})$ ($\beta_{i}$ is called the disease
transmission rate) and symmetric connectivity matrix $\mathcal{L}$. In
particular, the authors defined a basic reproduction number $\mathcal{R}_{0}$
and showed that $\mathcal{R}_{0}$ is a threshold parameter: the solution
converges to the disease free equilibrium if $\mathcal{R}_{0}<1$ and the model
has a unique endemic equilibrium (EE) (i.e. positive equilibrium) if
$\mathcal{R}_{0}>1$. Moreover, they showed that the disease component of the
EE converges to zero if $d_{S}\to 0$ when $\beta_{i}-\gamma_{i}$ changes sign
for $i\in\Omega$. Biologically, this result indicates that one may control the
disease by limiting the movement of susceptible people. In contrast, Li and
Peng [11] show that one may not eliminate the disease by limiting the movement
rate of infected people $d_{I}$. In [4], Chen et al. show that the above
results remain true if the connectivity matrix $\mathcal{L}=(L_{ij})$ is not
symmetric. Model (1.1) with the mass-action incidence transmission mechanism
$f_{i}=\beta_{i}S_{i}I_{i}$ has also received some attention recently (see
[12] for symmetric $\mathcal{L}$ and [16] for asymmetric $\mathcal{L}$). These
results show that whether the disease can be eliminated by limiting $d_{S}$ is
depending on the total population size $N$ and the infected people may
concentrate on the patches of the highest risk if $d_{I}\to 0$. Moreover, it
has been shown in [16] that the model with mass-action incidence transmission
mechanism may have multiple EE solutions when $\mathcal{R}_{0}<1$. Therefore,
the predictions on the consequences of disease control by limiting population
movement are depending on the chosen transmission mechanism, the population
size, and the control strategy. We remark that Gao et al. [6, 7] studied the
impact of population movement on the number of infected cases using model
(1.1) with the standard incidence transmission mechanism, and the related
reaction-diffusion models have attracted many research interests (see, e.g.,
[2, 13, 19, 3]).

The global dynamics of epidemic network models have attracted tremendous
research interests recently. In particular, a graph theoretical method based
on Kirkland’s matrix tree theorem to construct Lyapunov functions for network
models has been proposed in [8, 14]. The idea of this method is to take a
weighted combination of the Lyapunov functions for the decoupled local models
(i.e., the system without dispersal), and the weights are determined by a
Laplacian matrix that is related to the network structure. This method has
been applied to obtain the global stability of positive equilibrium for many
biological network models (see, e.g., [18, 20, 17, 9]).

The main objective of the current paper is to investigate the global dynamics
of model (1.1) with $d_{S}=0$ or $d_{I}=0$. As pointed out in the recent work
[15] on the reaction-diffusion counterpart of model (1.1), the asymptotic
profiles of the EE can reflect the consequences of disease control by limiting
population movement if the solutions of (1.1) always establish at some EE and
the disease dynamics happens in a faster time scale than the disease control.
However, if the disease control happens in a faster time scale, then the
global dynamics of (1.1) with $d_{S}=0$ or $d_{I}=0$ may better tell the
impact of the disease control strategies. As in [15], we will construct
Lyapunov type functions to study the global dynamics of the solutions. The
Lyapunov functions constructed by us do not use the graph-theoretical method,
and we will conquer the difficulties arising from asymmetric $\mathcal{L}$ by
proving a useful result, i.e., Theorem 2.4.

The rest of the paper is organized as follows. In section 2, we present some
preliminary results. In section 3, we study the global dynamics of (1.1) with
the mass-action incidence mechanism. In section 4, we study the global
dynamics of the model with the standard incidence mechanism.

##  2 Preliminaries

###  2.1 Notations

We will use a bold letter to represent a column vector in $\mathbb{R}^{n}$,
and the corresponding no-bold form with a subscript $j$ will be the $j$th-
component of the vector. Let $\bm{0}:=(0,\dots,0)^{T}$ and
$\bm{1}:=(1,\dots,1)^{T}$. For $\bm{X}\in\mathbb{R}^{n}$, define

${\bm{X}}_{m}:=\min_{j\in\Omega}X_{j},\quad{\bm{X}}_{M}:=\max_{j\in\Omega}X_{j},\quad\|\bm{X}\|_{1}:=\sum_{j=1}^{n}|X_{i}|,\quad\text{and}\quad\|\bm{X}\|_{\infty}:=\max_{j\in\Omega}|X_{j}|.$

We adopt the following convention: $\bm{X}\geq\bm{0}$ if $X_{i}\geq 0$ for all
$i\in\Omega$, $\bm{X}\gg\bm{0}$ if $X_{i}>0$ for all $i\in\Omega$, and
$\bm{X}>\bm{0}$ if $\bm{X}\geq\bm{0}$ and $\bm{X}\neq\bm{0}$. For
$\bm{X},\bm{Y}\in\mathbb{R}^{n}$, the Hadamard product of them is
$\bm{X}\circ\bm{Y}:=(X_{1}Y_{1},\cdots,X_{n}Y_{n})^{T}$ and we define
${\bm{X}}/{\bm{Y}}=(X_{1}/Y_{1},\cdots,$ $X_{n}/Y_{n})^{T}$ if $Y_{i}\neq 0$
for all $i\in\Omega$.

For $n\times n$ square matrix $A$, let $\sigma(A)$ be the spectrum of $A$,
$s(A)$ be the spectral bound of $A$, i.e.,

$s(A):=\max\\{\mathfrak{R}e(\lambda)\ :\ \lambda\in\sigma(A)\\},$

where $\mathfrak{R}e(\lambda)$ is the real part of $\lambda\in\mathbb{C}$, and
$\rho(A)$ be the spectral radius of $A$, i.e.,

$\rho(A):=\max\\{|\lambda|\ :\ \lambda\in\sigma(A)\\}.$

Throughout this work, we shall suppose that the following assumptions hold.

(A1) The matrix $\mathcal{L}=(L_{ij})_{i,j=1}^{n}$ is quasi-positive (i.e.,
$L_{ij}\geq 0$ for any $j\neq i$) and irreducible with
$L_{ii}=-\sum_{j\in\Omega,j\neq i}L_{ji}$ for each $i\in\Omega$.

(A2) The initial data $(\bm{S}^{0},\bm{I}^{0})$ satisfies
$\bm{S}^{0}\geq\bm{0}$, $\bm{I}^{0}>\bm{0}$ and
$\sum_{j\in\Omega}(S_{j}^{0}+I^{0}_{j})=N$, where $N$ is a fixed positive
constant.

(A3) $d_{S}\geq 0,d_{I}\geq 0$, $\bm{\beta}\gg\bm{0}$ and
$\bm{\gamma}\gg\bm{0}$.

By (A1) and the Perron-Frobenius theorem, we have

######  Lemma 2.1.

Suppose that (A1) holds. Then $s(\mathcal{L})=0$ is a simple eigenvalue of
$\mathcal{L}$, and there is an eigenvector $\bm{\alpha}$ associated with
$s(\mathcal{L})$ satisfying

$\mathcal{L}\bm{\alpha}=\bm{0},\quad\sum_{j\in\Omega}{\alpha}_{j}=1,\quad\text{and}\quad{\alpha}_{j}>0,\
\forall\ j\in\Omega,$ (2.1)

and $\bm{\alpha}$ is the unique nonnegative eigenvalue of $\mathcal{L}$
satisfying $\sum_{j\in\Omega}{\alpha}_{j}=1$.

Throughout the paper, $\bm{\alpha}$ is fixed and given as in Lemma 2.1. By
(A1)-(A3), for each initial data $(\bm{S}^{0},\bm{I}^{0})$, model (1.1) has a
unique nonnegative global solution, which leaves the following compact set
invariant

$\mathcal{E}:=\Big{\\{}(\bm{S},\bm{I})\in[\mathbb{R}^{n}_{+}]^{2}\ :\
\sum_{j\in\Omega}(S_{j}+I_{j})=N\Big{\\}}.$

Moreover, if $d_{I}=0$, then $I_{j_{0}}(0)=0$ for some $j_{0}\in\Omega$
implies that $I_{j_{0}}(t)=0$ for all $t>0$. However if $d_{I}>0$, by
$\bm{I}^{0}>\bm{0}$, $I_{j}(t)>0$ for all $t>0$ and $j\in\Omega$.

For convenience, define

${r}_{j}:=\frac{\gamma_{j}}{\beta_{j}},\quad j\in\Omega,$

which describes the local risk of patch $j$ when the network is decoupled.

###  2.2 Some useful results

We present an important result, Theorem 2.4, which will help us to construct
delicate Lyapunov functions later on. An $n\times n$ real matrix $A=(a_{ij})$
is called line-sum symmetric if $\sum_{j}a_{ij}=\sum_{j}a_{ji}$, i.e. the sum
of the $i$-th row equals the sum of the $i$-th column for all $i=1,\dots,n$.

######  Lemma 2.2.

Suppose that (A1) holds and $\mathcal{L}$ is line-sum symmetric. Then
$\sum_{i,j\in\Omega}L_{ij}X_{i}X_{j}\leq 0$ for any $\bm{X}\in\mathbb{R}^{n}$,
where the equality holds if and only if $\bm{X}$ is a multiple of $\bm{1}$.

###### Proof.

Since $\mathcal{L}$ is line-sum symmetric, $-L_{ii}=\sum_{j\neq
i,j\in\Omega}L_{ij}=\sum_{j\neq i,j\in\Omega}L_{ji}$ for any $i=1,\dots,n$.
For any $\bm{X}\in\mathbb{R}^{n}$, we have

$\displaystyle\sum_{i,j\in\Omega}L_{ij}X_{i}X_{j}$
$\displaystyle=\sum_{i,j\in\Omega,j\neq
i}L_{ij}X_{i}X_{j}+\sum_{i\in\Omega}L_{ii}X_{i}^{2}$
$\displaystyle=-\frac{1}{2}\sum_{i,j\in\Omega,j\neq
i}L_{ij}(X_{i}-X_{j})^{2}+\frac{1}{2}\sum_{i,j\in\Omega,j\neq
i}L_{ij}X_{i}^{2}+\frac{1}{2}\sum_{i,j\in\Omega,j\neq
i}L_{ij}X_{j}^{2}+\sum_{i\in\Omega}L_{ii}X_{i}^{2}$
$\displaystyle=-\frac{1}{2}\sum_{i,j\in\Omega,j\neq
i}L_{ij}(X_{i}-X_{j})^{2}+\frac{1}{2}\sum_{i,j\in\Omega,j\neq
i}L_{ij}X_{i}^{2}+\frac{1}{2}\sum_{i,j\in\Omega,j\neq
i}L_{ji}X_{i}^{2}+\sum_{i\in\Omega}L_{ii}X_{i}^{2}$
$\displaystyle=-\frac{1}{2}\sum_{i,j\in\Omega,j\neq
i}L_{ij}(X_{i}-X_{j})^{2}-\frac{1}{2}\sum_{i\in\Omega}L_{ii}X_{i}^{2}-\frac{1}{2}\sum_{i\in\Omega}L_{ii}X_{i}^{2}+\sum_{i\in\Omega}L_{ii}X_{i}^{2}$
$\displaystyle=-\frac{1}{2}\sum_{i,j\in\Omega,j\neq
i}L_{ij}(X_{i}-X_{j})^{2}\leq 0.$

Since $L$ is irreducible, the equality holds if and only if
$X_{1}=\dots=X_{n}$, i.e., $\bm{X}$ is a multiple of $\bm{1}$. ∎

######  Remark 2.3.

By Lemma 2.2, if $\mathcal{L}$ satisfies (A1) and is line-sum symmetric, then
$(\mathcal{L}+\mathcal{L}^{T})/2$ is negative semidefinite.

######  Theorem 2.4.

Suppose that (A1) holds. Let $\bm{\alpha}$ be defined as above and
$\bm{\theta}:=1/\bm{\alpha}$. Then $\sum_{i,j}\theta_{i}L_{ij}X_{i}X_{j}\leq
0$ for any $\bm{X}\in\mathbb{R}^{n}$, where the equality holds if and only if
$\bm{X}$ is a multiple of $\bm{\alpha}$.

###### Proof.

Let $D=\text{diag}(\bm{\alpha})$. Then, we have
$\mathcal{L}D=(L_{ij}\alpha_{j})_{i,j=1}^{n}$. The sum of the $j$-th column of
$\mathcal{L}D$ is $\sum_{i}L_{ij}\alpha_{j}=\alpha_{j}\sum_{i}L_{ij}=0$ for
any $1\leq j\leq n$, since the sum of each column of $\mathcal{L}$ is zero.
The sum of the $i$-th row of $\mathcal{L}D$ is $\sum_{j}L_{ij}\alpha_{j}=0$
for any $1\leq i\leq n$, since $\alpha$ is an eigenvector of $\mathcal{L}$
corresponding with eigenvalue zero. Therefore, $\mathcal{L}D$ is line-sum
symmetric, where the sum of each column and row is zero. By Lemma 2.2, we have
$\sum_{i,j}L_{ij}\alpha_{j}Y_{i}Y_{j}\leq 0$ for any
$\bm{Y}\in\mathbb{R}^{n}$, where the equality holds if and only if $\bm{Y}$ is
a multiple of $\bm{1}$.

We notice that
$\sum_{i,j}L_{ij}\alpha_{j}Y_{i}Y_{j}=\sum_{i,j}\theta_{i}L_{ij}(\alpha_{i}Y_{i})(\alpha_{j}Y_{j})$.
Since $\bm{\alpha}\gg\bm{0}$ and $\bm{Y}\in\mathbb{R}^{n}$ were arbitrary, we
have $\sum_{i,j}\theta_{i}L_{ij}X_{i}X_{j}\leq 0$ for any
$\bm{X}\in\mathbb{R}^{n}$, where the equality holds if and only if $\bm{X}$ is
a multiple of $\bm{\alpha}$. ∎

######  Remark 2.5.

Suppose that $L$ satisfies (A1). Let $\bm{\alpha}$ be a positive eigenvector
of $\mathcal{L}$ and $\bm{\theta}=1/\bm{\alpha}$. Let
$D_{1}=\text{diag}(\bm{\alpha})$ and $D_{2}=\text{diag}(\bm{\theta})$. By
Lemma 2.2 and Theorem 2.4, the two matrices
$(\mathcal{L}D_{1}+(\mathcal{L}D_{1})^{T})/2$ and
$(D_{2}\mathcal{L}+(D_{2}\mathcal{L})^{T})/2$ are negative semidefinite.

Next, we collect some useful results.

######  Lemma 2.6 ([16]).

Suppose that (A1) holds. Let $d$ be a positive number and
$\bm{F}:\mathbb{R}_{+}\to\mathbb{R}^{n}$ be a continuous map satisfying
$\|\bm{F}(t)\|_{1}\to 0$ as $t\to\infty$. Let $\bm{X}(t)$ be a bounded
solution of the system

$\begin{cases}\bm{X}^{\prime}(t)=d\mathcal{L}\bm{X}(t)+\bm{F}(t),\
t>0,\cr\bm{X}(0)=\bm{X}^{0}\in\mathbb{R}^{n}.\end{cases}$

Then
$\lim_{t\to\infty}\|\bm{X}(t)-(\sum_{j\in\Omega}X_{j}(t))\bm{\alpha}\|_{1}=0.$
In particular, if $\bm{F}(t)={\bf 0}$ for all $t\geq 0$, then
$\bm{X}(t)\to(\sum_{j\in\Omega}X^{0}_{j})\bm{\alpha}$ as $t\to\infty$.

The following is a Harnack’s type inequality.

######  Lemma 2.7.

[5, Lemma 3.1] Suppose that (A1) holds. Let $d>0$ and $\bm{M}\in
C(\mathbb{R}_{+}:[\mathbb{R}]^{n})$ such that

$m_{\infty}:=\sup_{t\geq 0}\|\bm{M}(t)\|_{\infty}<\infty.$

Then there is a positive number $c_{d,m_{\infty}}$ such that any nonnegative
solution $\bm{U}(t)$ of

$\bm{U}^{\prime}=d\mathcal{L}\bm{U}+\bm{M}(t)\circ\bm{U},\quad t>0,$ (2.2)

satisfies

$\|\bm{U}(t)\|_{\infty}\leq c_{d,m_{\infty}}\bm{U}_{m}(t),\quad\forall\ t\geq
1.$ (2.3)

The following result is well-known.

######  Lemma 2.8.

Suppose that $\phi:\mathbb{R}_{+}\to\mathbb{R}_{+}$ is Hölder continuous and
$\int_{0}^{\infty}\phi(t)dt<\infty$. Then $\lim_{t\to\infty}\phi(t)=0$.

##  3 Model with the mass-action incidence mechanism

In this section, we focus on the model with the mass-action incidence
mechanism, that is

$\begin{cases}\bm{S}^{\prime}=d_{S}\mathcal{L}\bm{S}-\bm{\beta}\circ\bm{S}\circ\bm{I}+\bm{\gamma}\circ\bm{I},\quad&t>0,\cr\bm{I}^{\prime}=d_{I}\mathcal{L}\bm{I}+\bm{\beta}\circ\bm{S}\circ\bm{I}-\bm{\gamma}\circ\bm{I},\quad&t>0.\end{cases}$
(3.1)

First, we will discuss the scenario of a complete limitation of the movement
of infected population, i.e $d_{I}=0$, and obtain Theorem 3.1. Then, we will
present our result on the scenario of a complete limitation of the movement of
susceptible population, i.e $d_{S}=0$, and establish Theorem 3.3.

We introduce the low-risk, moderate-risk, and high-risk $\tilde{H}^{-}$,
$\tilde{H}^{0}$, and $\tilde{H}^{+}$ patches as follows:

$\tilde{H}^{-}=\\{i\in\Omega:1<r_{i}/(N\alpha_{i})\\},\quad\tilde{H}^{0}:=\\{i\in\Omega:1=r_{i}/(N\alpha_{i})\\},\quad\text{and}\quad\tilde{H}^{+}:=\\{i\in\Omega:1>r_{i}/(N\alpha_{i})\\}.$

The sets $\tilde{H}^{-}$, $\tilde{H}^{0}$, and $\tilde{H}^{+}$ partition
$\Omega$ into disjoint subsets, and either of these sets can be empty. Note
that if the system is decoupled,
$R_{i}=N\alpha_{i}/r_{i}=N\alpha_{i}\beta_{i}/\gamma_{i}$ is the basic
reproduction number of patch $i$.

Given an initial data $(\bm{S}^{0},\bm{I}^{0})$, depending on the initial
distribution of the infected population, we introduce the set

$\Omega^{0}_{\bm{I}^{0}}:=\\{i\in\Omega:\
I^{0}_{i}=0\\}\quad\text{and}\quad\Omega^{+}_{\bm{I}^{0}}:=\\{i\in\Omega:\
I_{i}^{0}>0\\}$

and the quantity

$\tilde{r}_{m}=\min_{i\in\Omega_{\bm{I}^{0}}^{+}}\frac{r_{i}}{N\alpha_{i}}.$

The patches realizing the above minimum will be called the patches of the
highest-risk.

######  Theorem 3.1.

Suppose that (A1)-(A3) holds, $d_{S}>0$ and $d_{I}=0$. Then every solution
$(\bm{S}(t),\bm{I}(t))$ of (3.1) satisfies

$\lim_{t\to\infty}\|\bm{S}(t)-(\sum_{i\in\Omega}S_{i}(t)){\bm{\alpha}}\|_{\infty}=0\quad\text{and}\quad\lim_{t\to\infty}\sum_{j\in\tilde{H}^{-}{\cup\tilde{H}^{0}}\cup\Omega^{0}_{\bm{I}^{0}}}I_{j}(t)=0.$
(3.2)

Furthermore, the following conclusions hold.

  * (i)

If $\tilde{H}^{+}\cap\Omega^{+}_{\bm{I}^{0}}=\emptyset$, then
$\|\bm{S}(t)-{N{\bm{\alpha}}}\|_{\infty}\to 0$ and $\|\bm{I}(t)\|_{\infty}\to
0$ as $t\to\infty$;

  * (ii)

If $\tilde{H}^{+}\cap\Omega^{+}_{\bm{I}^{0}}\neq\emptyset$, then
$\|\bm{S}(t)-N\tilde{r}_{m}{\bm{\alpha}}\|_{\infty}\to 0$,
$\sum_{j\in\Omega}I_{j}(t)dx\to N-|\Omega|\tilde{r}_{m}$, and $I_{j}(t)\to 0$
on $\\{j\in\Omega:\ r_{j}/(N\alpha_{j})\neq\tilde{r}_{m}\\}$ as $t\to\infty$.

###### Proof.

Define

$V(\bm{S},\bm{I})=\sum_{i\in\Omega}\theta_{i}\left(\frac{S_{i}^{2}}{2}+r_{i}I_{i}\right),$

where $\theta_{i}=1/\alpha_{i}$ for $i\in\Omega$. It is easy to check that

$\dot{V}(\bm{S},\bm{I})=d_{S}\sum_{i,j\in\Omega}\theta_{i}L_{ij}S_{i}S_{j}-\sum_{i\in\Omega}\theta_{i}\beta_{i}\left(S_{i}-r_{i}\right)^{2}I_{i}.$

By Theorem 2.4, we have $\sum_{i,j\in\Omega}\theta_{i}L_{ij}S_{i}S_{j}\leq 0$
for all $\bm{S}\in\mathbb{R}^{k}$ and the equality holds if and only if
$\bm{S}$ is a multiple of $\bm{\alpha}$.

Denote $\mathcal{M}:=\\{(\bm{S},\bm{I})\in\mathcal{E}:\
\dot{V}(\bm{S},\bm{I})=0\\}$. Then,

$\mathcal{M}=\left\\{(\bm{S},\bm{I})\in\mathcal{E}:\ S_{i}=r_{i}\ \text{or}\
I_{i}=0\ \ \forall i\in\Omega\ \ \text{and}\ \bm{S}=a\bm{\alpha}\text{ for
some }a\in\mathbb{R}\right\\}.$

The maximal connected invariant subsets of $\mathcal{M}$ are

$\mathcal{M}_{i}:=\left\\{(\bm{S},\bm{I})\in\mathcal{E}:\
\bm{S}=\frac{r_{i}}{\alpha_{i}}\bm{\alpha}\ \ \text{and}\ \ I_{j}=0\ \
\text{if}\ \frac{r_{j}}{\alpha_{j}}\neq\frac{r_{i}}{\alpha_{i}}\ \forall
j\in\Omega\right\\},\ \ i\in\Omega,$

and

$\mathcal{M}_{0}:=\left\\{(\bm{S},\bm{I})\in\mathcal{E}:\
\bm{S}=\frac{N}{|\Omega|}\bm{\alpha}\ \ \text{and}\ \ \bm{I}=\bm{0}\right\\}.$

Hence by the invariance principle [10, Theorem 4.3.4],
$(\bm{S}(t),\bm{I}(t))\to\mathcal{M}_{i}$ for some $i=0,1,...,n$ as
$t\to\infty$. It follows from the equation of $\bm{S}$ in (3.1) and Lemma 2.6
that

$\lim_{t\to\infty}\Big{\|}\bm{S}(t)-{\Big{(}\sum_{j\in\Omega}S_{j}(t)\Big{)}}{\bm{\alpha}}\Big{\|}_{\infty}=0.$
(3.3)

If $\tilde{H}^{-}\neq\emptyset$, by its definition, we can choose
$\varepsilon>0$ such that

$\eta_{\varepsilon}:=\max_{i\in\tilde{H}^{-}}\left\\{\beta_{i}\left((N+\varepsilon)\alpha_{i}-r_{i}\right)\right\\}<0.$

By (3.3) and $\sum_{i\in\Omega}(S_{i}+I_{i})=N$ for all $t\geq 0$, there is
$t_{\varepsilon}>0$ such that

$\bm{S}(t)\leq\Big{(}N-\sum_{i\in\Omega}I_{i}+\varepsilon\Big{)}\bm{\alpha},\
\ \forall\;t\geq t_{\varepsilon}.$ (3.4)

It follows that

$I_{i}^{\prime}=\beta_{i}\left(S_{i}-r_{i}\right)I_{i}\leq\beta_{i}\left((N+\epsilon)\alpha_{i}-r_{i}\right)I_{i}\leq\eta_{\varepsilon}I_{i},\
\ \forall\ t\geq t_{\varepsilon},\ i\in\tilde{H}^{-}.$

Noticing $\eta_{\varepsilon}<0$, we have

$I_{i}(t)\leq I_{i}(t_{\varepsilon})e^{\eta_{\varepsilon}t}\to 0\ \ \text{as}\
t\to\infty,\quad\forall\;i\in\tilde{H}^{-}.$

Next, if $\tilde{H}^{0}\neq\emptyset$, by (3.4),

$I_{i}^{\prime}=\beta_{i}(S_{i}-r_{i})I_{i}\leq\beta_{i}((N-I_{i}+\varepsilon)\alpha_{i}-r_{i})I_{i}=\beta_{j}\alpha_{i}(\varepsilon-
I_{i})I_{i},\quad r\geq t_{\varepsilon},\ i\in\tilde{H}^{0}.$

Hence

$\limsup_{t\to\infty}I_{i}(t)\leq\varepsilon\quad\forall\ i\in\tilde{H}^{0}.$

Letting $\varepsilon\to 0$ in the last inequality, we obtain that $I_{i}(t)\to
0$ as $t\to\infty$ for each $i\in\tilde{H}^{0}.$ Since $I_{i}(t)=0$ for all
$t\geq 0$ for all $i\in\Omega^{0}_{\bm{I}^{0}}$, (3.2) holds.

It remains to prove (i)-(ii). We first show that

$\limsup_{t\to\infty}\sum_{i\in\Omega}I_{i}(t)\leq N(1-\tilde{r}_{m})_{+}.$
(3.5)

Since $I_{i}(t)=0$ for all $t\geq 0$ for all $i\in\Omega^{0}_{\bm{I}^{0}}$, it
suffices to prove (3.5) with $\Omega$ replaced by $\Omega^{+}_{\bm{I}^{0}}$.
To see this, by (3.4),

$\displaystyle I_{i}^{\prime}=\beta_{i}\left(S_{i}-r_{i}\right)I_{i}$
$\displaystyle\leq\beta_{i}\Big{(}\Big{(}N-\sum_{i\in\Omega^{+}_{\bm{I}^{0}}}I_{i}+\epsilon\Big{)}\alpha_{i}-r_{i}\Big{)}I_{i}$
$\displaystyle\leq\alpha_{i}\beta_{i}\Big{(}N(1-\tilde{r}_{m})_{+}+\varepsilon-\sum_{i\in\Omega^{+}_{\bm{I}^{0}}}I_{i}\Big{)}I_{i},\
\ \forall\ t\geq t_{\varepsilon},\ i\in\Omega^{+}_{\bm{I}^{0}}.$

Define

$F(t)=\frac{\sum_{i\in\Omega^{+}_{\bm{I}^{0}}}\alpha_{i}\beta_{i}I_{i}(t)}{\sum_{i\in\Omega^{+}_{\bm{I}^{0}}}I_{i}(t)},\
\ t\geq 0.$

It then follows that

$\frac{d}{dt}\sum_{i\in\Omega^{+}_{\bm{I}^{0}}}I_{i}\leq
F(t)\Big{(}N(1-\tilde{r}_{m})_{+}+\varepsilon-\sum_{i\in\Omega^{+}_{\bm{I}^{0}}}I_{i}\Big{)}\sum_{i\in\Omega^{+}_{\bm{I}^{0}}}I_{i},\
\ \forall\ t\geq t_{\varepsilon}$

Noticing $0<\min\\{\alpha_{i}\beta_{i}:\ x\in\Omega^{+}_{\bm{I}^{0}}\\}\leq
F(t)\leq\max\\{\alpha_{i}\beta_{i}:\ x\in\Omega^{+}_{\bm{I}^{0}}\\}$, we have

$\limsup_{t\to\infty}\sum_{i\in\Omega^{+}_{\bm{I}^{0}}}I_{i}(t)\leq
N(1-\tilde{r}_{m})_{+}+\varepsilon.$

Since $\varepsilon>0$ was arbitrary, (3.5) holds.

(i) If $H^{+}\cap\Omega^{+}_{\bm{I}^{0}}=\emptyset$, then
$(1-\tilde{r}_{m})_{+}=0$. By (3.5), we have
$\lim_{t\to\infty}\|\bm{I}(t)\|_{\infty}=0$. Since
$\sum_{i\in\Omega}(S_{i}+I_{i})=N$ for all $t\geq 0$,
$\sum_{i\in\Omega}S_{i}(t)\to N$ as $t\to\infty$. Then it follows from (3.3)
that $\|\bm{S}(t)-{N{\bm{\alpha}}}\|_{\infty}\to 0$ as $t\to\infty$.

(ii) If $H^{+}\cap\Omega^{+}_{\bm{I}^{0}}\neq\emptyset$, then
$1>\tilde{r}_{m}$. By (3.5) and $\sum_{i\in\Omega}(S_{i}(t)+I_{i}(t))=N$ for
all $t\geq 0$,

$\liminf_{t\to\infty}\sum_{i\in\Omega}S_{i}(t)\geq N{\tilde{r}_{m}}.$

By (3.3), we have

$\liminf_{t\to\infty}S_{i}(t)\geq N{\tilde{r}_{m}}\alpha_{i},\ \ \forall\
i\in\Omega.$

Now we claim that there exists $i_{0}\in\Omega$ such that

$\liminf_{t\to\infty}S_{i_{0}}(t)=N{\tilde{r}_{m}}\alpha_{i_{0}}.$ (3.6)

Assume to the contrary that the claim is false. Then there exist
$\varepsilon^{\prime},t_{\varepsilon^{\prime}}>0$ such that

$S_{i}(t)\geq N{\tilde{r}_{m}}\alpha_{i}+\varepsilon^{\prime},\ \ \forall\
t\geq t_{\varepsilon^{\prime}},\ i\in\Omega.$

Let $i_{0}\in\Omega^{+}_{\bm{I}^{0}}$ be such that
${\tilde{r}_{m}}=r_{i_{0}}/(N\alpha_{i_{0}})$. Then,

$I^{\prime}_{i_{0}}=\beta_{i_{0}}(S_{i_{0}}-r_{i_{0}})I_{i_{0}}\geq\varepsilon^{\prime}\beta_{i_{0}}I_{i_{0}},\
\ \forall\;t\geq t_{\varepsilon^{\prime}}.$

This implies $I_{i_{0}}(t)\geq
I_{i_{0}}(t_{\varepsilon^{\prime}})e^{\varepsilon^{\prime}\beta_{i_{0}}(t-t_{\varepsilon^{\prime}})}$,
$t\geq t_{\varepsilon^{\prime}}$, which contradicts the boundedness of the
solution. This proves (3.6).

By (3.6), $(\bm{S}(t),\bm{I}(t))\to\mathcal{M}_{i_{0}}$ as $t\to\infty$. This
proves (ii). ∎

######  Remark 3.2.

By Theorem 3.1, if $\tilde{H}^{+}\cap\Omega^{+}_{\bm{I}^{0}}=\emptyset$, the
disease can be eliminated by limiting the movement of infected people when the
model has mass-action mechanism. However if
$\tilde{H}^{+}\cap\Omega^{+}_{\bm{I}^{0}}\neq\emptyset$, the infected people
will concentrate on the patches of the highest risk, $\\{j\in\Omega:\
r_{j}/(N\alpha_{j})\neq\tilde{r}_{m}\\}$, when limiting the movement of them.
These predictions largely agree with those in [12, 16] for the model with
mass-action mechanism, which are based on the asymptotic profiles of the EE as
$d_{I}\to 0$. The main difference is that our results are depending on the
distribution of initial data $\bm{I}^{0}$.

Next, we examine the scenario of a total restriction of the movement of the
susceptible population. In this direction, our result reads as follows.

######  Theorem 3.3.

Suppose that (A1)-(A3) holds, $d_{S}=0$ and $d_{I}>0$. Then any solution
$(\bm{S}(t),\bm{I}(t))$ of (3.1) satisfies
$\|\bm{S}(t)-\bm{S}^{*}\|_{\infty}\to 0$ as $t\to\infty$ for some
$\bm{S}^{*}\gg\bm{0}$. Furthermore, exactly one of the following statements
holds:

  * (i)

$\bm{S}^{*}=\bm{\lambda}^{*}\circ\bm{S}^{0}+({\bm{1}}-\bm{\lambda}^{*})\circ\bm{r}$
for some $\bm{\lambda}^{*}\in\mathbb{R}^{n}$ with
${\bm{0}}\ll\bm{\lambda}^{*}\ll{\bm{1}}$, $\sum_{j\in\Omega}S^{*}_{j}=N$ and
$s(d_{I}\mathcal{L}+{\rm
diag}(\bm{\beta}\circ\bm{\lambda}^{*}\circ(\bm{S}^{0}-\bm{r})))\leq 0$, and
$\bm{I}(t)\to\bm{0}$ as $t\to\infty$.

  * (ii)

$\bm{S}^{*}=\bm{r}$ and $\|\bm{I}(t)-\bm{I}^{*}\|_{\infty}\to 0$ as
$t\to\infty$, where
$\bm{I}^{*}=(N-\sum_{j\in\Omega}r_{j})\bm{\alpha}\gg\bm{0}$.

Moreover, (i) holds if $N\leq\sum_{j\in\Omega}r_{j}$ while (ii) holds if
$N>N^{*}_{\bm{S}^{0},\bm{r}}$, where $N^{*}_{\bm{S}^{0},\bm{r}}$ is given by

$N^{*}_{\bm{S}^{0},\bm{r}}:=\sup\left\\{\bm{\lambda}^{*}\cdot\bm{S}^{0}+(1-\bm{\lambda}^{*})\cdot\bm{r}:\
{\bm{0}}\leq\lambda^{*}\leq{\bf 1}\quad\text{and}\quad s(d_{I}\mathcal{L}+{\rm
diag}(\bm{\beta}\circ\bm{\lambda}^{*}\circ(\bm{S}_{0}-\bm{r})))\leq
0\right\\}.$ (3.7)

###### Proof.

Let $J_{i}(t):=\int_{0}^{t}I_{i}(\tau)d\tau$ for $i\in\Omega$ and $t\geq 0$
and define

$\tilde{J}_{i}:=\int_{0}^{\infty}I_{i}(\tau)d\tau=\lim_{t\to\infty}J_{i}(t)\in(0,\infty].$

We consider two cases.
Case 1. $\tilde{J}_{i_{0}}<\infty$ for some $i_{0}\in\Omega$. By Lemma 2.7,
there exists $C_{1}>0$ such that

$\bm{I}_{M}(t)\leq C_{1}I_{i}(t),\quad\forall\;t\geq 1,i\in\Omega.$ (3.8)

Therefore,

$\tilde{J}_{i}=\tilde{J}_{i}(1)+\int_{1}^{\infty}I_{i}(\tau)d\tau\leq\tilde{J}_{i}(1)+C_{1}\int_{1}^{\infty}I_{i_{0}}(\tau)d\tau\leq{J}_{i}(1)+C_{1}\tilde{J}_{i_{0}},\quad\forall\;i\in\Omega.$

Hence, $|\tilde{J}_{i}|<\infty$ for all $i\in\Omega$ and
$\|\bm{\tilde{J}}\|_{\infty}<\infty$. Since $|I^{\prime}_{i}|$ is uniformly
bounded for $t\in(0,\infty)$ and $\tilde{J}_{i}<\infty$, it follow from Lemma
2.8 that $I_{i}(t)\to 0$ as $t\to\infty$ for each $i\in\Omega$.

By (3.8) again, we have

$\|\tilde{\bm{J}}-\bm{J}(t)\|_{\infty}\leq
C_{1}\int_{t}^{\infty}I_{i_{0}}(\tau)d\tau\to 0\ \ \text{as}\ \ t\to\infty.$

By the first equation of (3.1), we have

$\frac{d}{dt}(S_{i}-r_{i})=-\beta_{i}(S_{i}-r_{i})I_{i},\quad i\in\Omega,t>0.$

It follows that

$S_{i}-r_{i}=(S_{i}^{0}-r_{i})e^{-\beta_{i}J_{i}(t)},\quad i\in\Omega,t>0.$
(3.9)

This implies that

$S_{i}\to
S_{i}^{*}:=r_{i}+(S_{i}^{0}-r_{i})e^{-\beta_{i}\tilde{J}_{i}}=e^{-\beta_{i}\tilde{J}_{i}}S_{i}^{0}+(1-e^{-\beta_{i}\tilde{J}_{i}})r_{i}\quad\text{as}\quad
t\to\infty,\quad i\in\Omega.$

Define $\lambda_{i}^{*}=e^{-\beta_{i}\tilde{J}_{i}}$, $i\in\Omega$. Then,
$\bm{0}\ll\bm{\lambda}\ll\bm{1}$,
$\bm{S}^{*}=\bm{\lambda}^{*}\circ\bm{S}^{0}+({\bm{1}}-\bm{\lambda}^{*})\circ\bm{r}$,
and $\bm{S}(t)\to\bm{S}^{*}$ as $t\to\infty$. Noticing
$\sum_{i\in\Omega}(S_{i}+I_{i})=N$ for any $t\geq 0$ and $I_{i}(t)\to 0$ as
$t\to\infty$, we have $\sum_{i}S_{i}^{*}=N$.

Let $\bm{\psi}$ be the positive eigenvector corresponding with
$s:=s(d_{I}\mathcal{L}^{T}+{\rm
diag}(\bm{\beta}\circ\bm{\lambda}^{*}\circ(\bm{S}^{0}-\bm{r})))$ with
$\bm{\psi}_{M}=1$. Note that $s=s(d_{I}\mathcal{L}+{\rm
diag}(\bm{\beta}\circ\bm{\lambda}^{*}\circ(\bm{S}^{0}-\bm{r})))$. Then, we
have

$\displaystyle\frac{d}{dt}\sum_{i\in\Omega}\psi_{i}I_{i}$ $\displaystyle=$
$\displaystyle
d_{I}\sum_{i,j\in\Omega}L_{ij}\psi_{i}I_{j}+\sum_{i\in\Omega}(\beta_{i}S_{i}-\gamma_{i})\psi_{i}I_{i}$
$\displaystyle=$ $\displaystyle
d_{I}\sum_{i,j\in\Omega}L_{ji}\psi_{j}I_{i}+\sum_{i\in\Omega}(\beta_{i}S_{i}-\gamma_{i})\psi_{i}I_{i}$
$\displaystyle=$ $\displaystyle
s\sum_{i\in\Omega}\psi_{i}I_{i}+\sum_{i\in\Omega}\beta_{i}(S_{i}-S_{i}^{*})\psi_{i}I_{i}\geq\left(s-\bm{\beta}_{M}\|\bm{S}-\bm{S}^{*}\|_{\infty}\right)\sum_{i\in\Omega}\psi_{i}I_{i}.$

Therefore,

$N\geq\sum_{i\in\Omega}\psi_{i}I_{i}\geq
e^{\int_{0}^{t}(s-\bm{\beta}_{M}\|\bm{S}(\tau)-\bm{S}^{*}\|_{\infty})d\tau}\sum_{i\in\Omega}\psi_{i}I_{i}^{0},\quad
t\geq 0.$

Taking log for both sides and dividing by $t$, we obtain

$\frac{\ln
N-\ln\sum_{i\in\Omega}\psi_{i}I_{i}^{0}}{t}+\beta_{M}\frac{\int_{0}^{t}\|\bm{S}(\tau)-\bm{S}^{*}\|_{\infty}d\tau}{t}\geq
s.$

Letting $t\to\infty$ and noticing $\bm{S}(t)\to\bm{S}^{*}$ as $t\to\infty$, we
have $s\leq 0$.

Case 2. $\tilde{J}_{i}=\infty$ for all $i\in\Omega$. Fix $i_{1}\in\Omega$.
Then, $J_{i_{1}}(t)\to\infty$ as $t\to\infty$. By (3.8) and (3.9),

$\|\bm{S}(t)-\bm{r}\|_{\infty}\leq\|\bm{S}_{0}-\bm{r}\|_{\infty}e^{-\frac{\bm{\beta}_{m}}{C_{1}}J_{i_{1}}(t)}\to
0\quad\text{as}\ t\to\infty.$

Noticing $\sum_{i\in\Omega}(S_{i}+I_{i})=N$ for all $t\geq 0$,
$\sum_{i\in\Omega}I_{i}\to N-\sum_{i\in\Omega}r_{i}$ as $t\to\infty$. By Lemma
2.6 and the second equation of (3.1), we have

$\lim_{t\to\infty}\Big{\|}\bm{I}(t)-{\Big{(}N-\sum_{i\in\Omega}r_{i}\Big{)}}\bm{\alpha}\Big{\|}_{\infty}=0.$

This proves (ii).

Finally, it is easy to see that $N\leq\sum_{i\in\Omega}r_{i}$ implies (i). If
(i) holds, then
$N=\sum_{i\in\Omega}S^{*}_{i}=\sum_{i\in\Omega}(\lambda^{*}_{i}S_{i}^{0}+(1-\lambda^{*}_{i})r_{i})$
for some $\bm{\lambda}\in\mathbb{R}^{n}$ with
$\bm{0}\ll\bm{\lambda}^{*}\ll\bm{1}$ and $s\leq 0$. So, we have $N\leq
N^{*}_{\bm{S}^{0},\bm{r}}$. Therefore, if $N>N^{*}_{\bm{S}^{0},\bm{r}}$, then
(ii) holds. ∎

######  Remark 3.4.

By Theorem 3.3, whether limiting the movement of susceptible people can
eliminate the disease is depending on the total population $N$: if
$N\leq\sum_{j\in\Omega}r_{j}$, then disease can be controlled; if $N$ is
large, then it cannot be controlled by limiting the movement of susceptible
people. These predictions agree with those in [12, 16] where the predictions
are based on the asymptotic profiles of the EE as $d_{S}\to 0$.

######  Remark 3.5.

Alternatively, we may use a Lyapunov function to prove Theorem 3.3. Define

$V(\bm{S},\bm{I})=\sum_{i\in\Omega}\left(\frac{\beta_{i}S_{i}^{2}}{2\gamma_{i}}+I_{i}\right).$

Then it is easy to check that

$\dot{V}(\bm{S},\bm{I})=-\sum_{i\in\Omega}\frac{1}{\gamma_{i}}\left(\beta_{i}S_{i}-\gamma_{i}\right)^{2}I_{i}\leq
0.$

##  4 Model with the standard incidence mechanism

In the current section, we focus on the model with standard incidence
mechanism, that is

$\begin{cases}\bm{S}^{\prime}=d_{S}\mathcal{L}\bm{S}-(\bm{\beta}\circ\bm{S}\circ\bm{I})/(\bm{S}+\bm{I})+\bm{\gamma}\circ\bm{I}&t>0,\cr\bm{I}^{\prime}=d_{I}\mathcal{L}\bm{I}+(\bm{\beta}\circ\bm{S}\circ\bm{I})/(\bm{S}+\bm{I})-\bm{\gamma}\circ\bm{I}&t>0.\end{cases}$
(4.1)

We will study the global dynamics of the model when $d_{S}=0$ or $d_{I}=0$ as
in the previous section.

Here, we define the low-risk, moderate-risk, and high-risk patches ${H}^{+}$,
$H^{0}$ and ${H}^{-}$ as

$H^{-}:=\\{i\in\Omega:1<r_{i}\\},\quad
H^{0}:=\\{i\in\Omega:1=r_{i}\\},\quad\text{and}\quad
H^{+}:=\\{i\in\Omega:1>r_{i}\\}.$

The sets $H^{-}$, $H^{0}$, and $H^{+}$ partition $\Omega$ into disjoint
subsets, and either of these sets can be empty. Note that
$1/r_{i}=\beta_{i}/\gamma_{i}$ is the basic reproduction number of patch $i$
when the system is decoupled.

######  Theorem 4.1.

Suppose that (A1)-(A3) holds, $d_{I}=0$ and $d_{S}>0$. Let $(\bm{S},\bm{I})$
be the solution of (4.1). Then,

$\lim_{t\to\infty}\bm{S}(t)=k\bm{\alpha},\quad\lim_{t\to\infty}I_{i}(t)=0,\quad\forall\;i\in
H_{1},\quad\text{and}\quad\lim_{t\to\infty}I_{i}(t)=\frac{\beta_{i}-\gamma_{i}}{\gamma_{i}}k\alpha_{i},\quad\forall\;i\in
H_{2},$

where $H_{1}:=H^{0}\cup H^{-}\cup\Omega_{\bm{I}^{0}}^{0}$,
$H_{2}:=H^{+}\cap\Omega^{+}_{\bm{I}^{0}}$ and $k=N/({1+\sum_{i\in
H_{2}}\frac{\beta_{i}-\gamma_{i}}{\gamma_{i}}\alpha_{i}}).$

###### Proof.

First, we claim that there exists $\underline{S}>0$ such that
$\liminf_{t\to\infty}S_{i}(t)\geq\underline{S}$ for all $i\in\Omega$. To see
this, by the equation of $\bm{S}$,

$\displaystyle\frac{d}{dt}\sum_{i\in\Omega}S_{i}(t)$
$\displaystyle=\sum_{i\in\Omega}\gamma_{i}I_{i}-\sum_{i\in\Omega}\beta_{i}\frac{I_{i}}{I_{i}+S_{i}}S_{i}\vskip
6.0pt plus 2.0pt minus
2.0pt\displaystyle\geq\bm{\gamma}_{m}\sum_{i\in\Omega}I_{i}-\bm{\beta}_{M}\sum_{i\in\Omega}S_{i}\vskip
6.0pt plus 2.0pt minus 2.0pt$
$\displaystyle=\bm{\gamma}_{m}\left(N-\sum_{i\in\Omega}S_{i}\right)-\bm{\beta}_{M}\sum_{i\in\Omega}S_{i}\vskip
6.0pt plus 2.0pt minus
2.0pt=\bm{\gamma}_{m}N-(\bm{\beta}_{M}+\bm{\gamma}_{m})\sum_{i\in\Omega}S_{i}.$

So we can choose $t_{0}>0$ such that

$\sum_{i\in\Omega}S_{i}(t)\geq\frac{\bm{\gamma}_{m}N}{2(\bm{\beta}_{M}+\bm{\gamma}_{m})},\quad\forall\;t\geq
t_{0}.$ (4.2)

By Lemma 2.7, there exists $c_{0}>0$ such that

$e^{d_{S}\mathcal{L}}\bm{S}(t)\geq
c_{0}\|e^{d_{S}\mathcal{L}}\bm{S}(t)\|_{1}\bm{1},\quad\forall\;t\geq 0.$ (4.3)

By the equation $\bm{S}$ again,

$\bm{S}^{\prime}\geq d_{S}\mathcal{L}\bm{S}-\bm{\beta}_{M}\bm{S}.$

So by the comparison principle and (4.2)-(4.3), we have

$\displaystyle\bm{S}(t+1)$ $\displaystyle\geq
e^{-\bm{\beta}_{M}}e^{d_{S}\mathcal{L}}\bm{S}(t)\geq
c_{0}e^{-\bm{\beta}_{M}}\|e^{d_{S}\mathcal{L}}\bm{S}(t)\|_{1}\bm{1}$
$\displaystyle=c_{0}e^{-\bm{\beta}_{M}}\|\bm{S}(t)\|_{1}\bm{1}\geq
c_{0}e^{-\bm{\beta}_{M}}\frac{\bm{\gamma}_{m}N}{2(\bm{\beta}_{M}+\bm{\gamma}_{m})}\bm{1},\
\ \forall\;t\geq t_{0}.$

It follows that

$S_{i}(t)\geq\underline{S}:=c_{0}e^{-\bm{\beta}_{M}}\frac{\bm{\gamma}_{m}N}{2(\bm{\beta}_{M}+\bm{\gamma}_{m})},\quad\forall\;i\in\Omega,\
t\geq t_{0}+1.$ (4.4)

Define

$V(\bm{S},\bm{I})=\frac{1}{2}\sum_{i\in\Omega}\theta_{i}S_{i}^{2}+\frac{1}{2}\sum_{i\in
H_{2}}\theta_{i}\frac{\gamma_{i}}{\beta_{i}-\gamma_{i}}I_{i}^{2},$

where $\bm{\theta}=1/\bm{\alpha}$. Then, we have

$\dot{V}(\bm{S},\bm{I})=\sum_{i,j\in\Omega}\theta_{i}L_{ij}S_{i}S_{j}-\sum_{i\in
H_{1}}\theta_{i}\frac{(\beta_{i}-\gamma_{i})S_{i}-\gamma_{i}I_{i}}{S_{i}+I_{i}}S_{i}I_{i}-\sum_{i\in
H_{2}}\frac{\theta_{i}}{\beta_{i}-\gamma_{i}}\frac{[(\beta_{i}-\gamma_{i})S_{i}-\gamma_{i}I_{i}]^{2}I_{i}}{S_{i}+I_{i}}.$
(4.5)

By Theorem 2.4, $\sum_{i,j\in\Omega}\theta_{i}L_{ij}S_{i}S_{j}\leq 0$, where
the equality holds if and only if $\bm{S}$ is a multiple of $\bm{\alpha}$.

It is easy to see that if $I_{i}(0)=0$ for some $i\in\Omega$ then $I_{i}(t)=0$
for all $t>0$. If $i\in H^{-}\cup H^{0}$, then $\beta_{i}\leq\gamma_{i}$ and,
by the equation of $\bm{I}$,

$I^{\prime}_{i}(t)=\frac{((\beta_{i}-\gamma_{i})S_{i}-\gamma_{i}I_{i})I_{i}}{S_{i}+I_{i}}\leq\frac{-\gamma_{i}I_{i}^{2}}{S_{i}+I_{i}}\leq-\frac{\gamma_{i}}{N}I_{i}^{2}.$

It follows that $I_{i}(t)\to 0$ as $t\to\infty$ for all $i\in H_{1}$.
Moreover, integrating the equation of $I_{i}$ on $(0,\infty)$, we find that

$-I^{0}_{i}=\int_{0}^{\infty}I_{i}^{\prime}(t)dt=\int_{0}^{\infty}\frac{(\beta_{i}-\gamma_{i})S_{i}-\gamma_{i}I_{i}}{S_{i}+I_{i}}I_{i}dt,\quad\forall\;i\in
H^{-}\cup H^{0}.$

Noticing $(\beta_{i}-\gamma_{i})S_{i}-\gamma_{i}I_{i}\leq 0$ for $i\in
H^{-}\cup H^{0}$, we have

$0\leq-\int_{0}^{\infty}\frac{(\beta_{i}-\gamma_{i})S_{i}-\gamma_{i}I_{i}}{S_{i}+I_{i}}S_{i}I_{i}dt<\infty.$

So integrating (4.5) over $(0,\infty)$, we find that

$0\leq-\int_{0}^{\infty}\sum_{i,j\in\Omega}\theta_{i}L_{ij}S_{i}S_{j}dt<\infty\quad\text{and}\quad\int_{0}^{\infty}\frac{[(\beta_{i}-\gamma_{i})S_{i}-\gamma_{i}I_{i}]^{2}I_{i}}{S_{i}+I_{i}}dt<\infty,\quad\forall\;i\in
H_{2}.$ (4.6)

Since $\underline{S}/2\leq S_{i}(t)+I_{i}(t)<N$ for all $t\gg 1$, we have
$\sup_{t\geq 0}|S_{i}^{\prime}(t)|<\infty$ and $\sup_{t\geq
0}|I_{i}^{\prime}(t)|<\infty$. We then deduce from (4.6) and Lemma 2.8 that

$\lim_{t\to\infty}\sum_{i,j\in\Omega}\theta_{i}L_{ij}S_{i}S_{j}=0$ (4.7)

and

$\lim_{t\to\infty}\frac{[(\beta_{i}-\gamma_{i})S_{i}-\gamma_{i}I_{i}]^{2}I_{i}}{S_{i}+I_{i}}=0,\quad\forall\;i\in
H_{2}.$ (4.8)

Observe that

$I_{i}^{\prime}(t)=((\beta_{i}-\gamma_{i})S_{i}-\gamma_{i}I_{i})\frac{I_{i}}{S_{i}+I_{i}}\geq\Big{(}\frac{\underline{S}}{2}(\beta_{i}-\gamma_{i})-\gamma_{i}I_{i}\Big{)}\frac{I_{i}}{S_{i}+I_{i}},\quad
t\gg 1,\ i\in H_{2}.$

Therefore,
$\liminf_{t\to\infty}I_{i}(t)\geq\frac{\underline{S}(\beta_{i}-\gamma_{i})}{2\gamma_{i}}>0$
for $i\in H_{2}$. As a result, we obtain from (4.8) that

$\lim_{t\to\infty}[(\beta_{i}-\gamma_{i})S_{i}-\gamma_{i}I_{i}]^{2}=0,\quad\forall\;i\in
H_{2}.$ (4.9)

Since
$\bm{0}<\underline{S}\bm{1}\leq\liminf_{t\to\infty}\bm{S}(t)\leq\limsup_{t\to\infty}\bm{S}(t)\leq
N\bm{1}$, possibly after passing to a subsequence, we may suppose that
$\bm{S}(t)\to\bm{S}^{*}$ as $t\to\infty$. Thanks to (4.7), we have
$\sum_{i,j\Omega}\theta_{i}L_{ij}S_{i}^{*}S_{j}^{*}=0$. In view of Theorem
2.4, $\bm{S}^{*}=k^{*}\bm{\alpha}$ for some positive number $k^{*}$. By (4.9),
we have that $I_{i}(t)\to\frac{k^{*}(\beta_{i}-\gamma_{i})}{\gamma_{i}}$ as
$t\to\infty$, for each $i\in H_{2}$. Recalling that
$\|\bm{S}(t)+\bm{I}(t)\|_{1}=N$ for all $t\geq 0$, then $N=k^{*}(1+\sum_{j\in
H_{2}}(\beta_{j}-\gamma_{j})\alpha_{j}/\gamma_{j})$, which yields that
$k^{*}=N/((1+\sum_{j\in H_{2}}(\beta_{j}-\gamma_{j})\alpha_{j}/\gamma_{j}))$.
Since $k^{*}$ is independent of the chosen subsequence, the desired result
follows.

∎

######  Remark 4.2.

By Theorem 4.1, the infected people will live exactly on
$H^{+}\cap\Omega^{+}_{\bm{I}^{0}}$ (i.e. the high-risk patches with positive
initial infected cases) for the model with standard incidence mechanism when
the dispersal of infected people is limited. As before, these predictions
agree with previous studies for the model with standard incidence mechanism in
[11, 4], which are based on the asymptotic profiles of the EE as $d_{I}\to 0$.
The main difference is that our results are depending on the distribution of
initial data $\bm{I}^{0}$.

Next, we examine the scenario of a total restriction of the movement of the
susceptible population. In this direction, our result reads as follows.

######  Theorem 4.3.

Suppose that (A1)-(A3) holds, $d_{S}=0$ and $d_{I}>0$. Let
$(\bm{S}(t),\bm{I}(t))$ be the solution of (1.1). Then the following
conclusions hold:

  1. (i)

If $H^{-}\cup H^{0}$ is nonempty, then $\bm{I}(t)\to\bm{0}$ as $t\to\infty$
and $\liminf_{t\to\infty}S_{j}(t)>0$ for $j\in H^{-}\cup H^{0}$.

  2. (ii)

If $H^{+}=\Omega$, then
$(\bm{S}(t),\bm{I}(t))\to(I^{*}(\bm{r}\circ\bm{\alpha}/(\bm{1}-\bm{r})),I^{*}\bm{\alpha})$
as $t\to\infty$, where $I^{*}=N/(1+\|\bm{1}/(\bm{1}-\bm{r})\|_{1})$.

###### Proof.

Define

$V(\bm{S},\bm{I})=\frac{1}{2}\sum_{i\in\Omega}\theta_{i}\frac{(1-r_{i})}{r_{i}}S_{i}^{2}+\frac{1}{2}\sum_{i\in\Omega}\theta_{i}I_{i}^{2},$

where $\bm{\theta}=\bm{1}/\bm{\alpha}$. Then
$-\frac{N^{2}}{2}\sum_{i\in\Omega}\theta_{i}\leq V(\bm{S},\bm{I})\leq
N^{2}(1+\|\bm{1}/\bm{r}\|_{\infty})\sum_{i\in\Omega}\theta_{i}$ and

$\dot{V}(\bm{S},\bm{I})=\sum_{i,j\in\Omega}\theta_{i}L_{ij}I_{i}I_{j}-\sum_{i\in\Omega}\frac{\theta_{i}\beta_{i}^{2}}{\gamma_{i}}\big{(}(1-r_{i})S_{i}-r_{i}I_{i}\big{)}^{2}\frac{I_{i}}{S_{i}+I_{i}}.$

By Theorem 2.4, $\sum_{i,j\in\Omega}\theta_{i}L_{ij}I_{i}I_{j}\leq 0$, where
the equality holds if and only if $\bm{I}$ is a multiple of $\bm{\alpha}$.
Denote $\mathcal{N}:=\\{(\bm{S},\bm{I})\in\mathcal{E}:\
\dot{V}(\bm{S},\bm{I})=0\\}$. Then,

$\mathcal{N}=\left\\{(\bm{S},\bm{I})\in\mathcal{E}:\
\big{(}(1-r_{i})S_{i}-r_{i}I_{i}\big{)}^{2}I_{i}=0\ \ \forall i\in\Omega\ \
\text{and}\ \bm{I}=a\bm{\alpha}\text{ for some }a\in\mathbb{R}_{+}\right\\}.$

From this point, we distinguish two cases.

(i) Suppose that $H^{-}\cup H^{0}\neq\emptyset$. In this case, the maximal
connected invariant subset of $\mathcal{N}$ is

$\mathcal{N}_{1}:=\left\\{(\bm{S},\bm{I})\in\mathcal{E}:\
\bm{I}=\bm{0}\right\\}.$

Hence by the invariance principle [10, Theorem 4.3.4],
$(\bm{S}(t),\bm{I}(t))\to\mathcal{N}_{1}$ as $t\to\infty$. Thus
$\bm{I}(t)\to\bm{0}$ as $t\to\infty$.

Next, fix $i\in H^{-}\cup H^{0}$. Then

$S_{i}^{\prime}=(\gamma_{i}I_{i}-(\beta_{i}-\gamma_{i})S_{i})\frac{I_{i}}{S_{i}+I_{i}}\geq\frac{\gamma_{i}I_{i}^{2}}{S_{i}+I_{i}}\geq\frac{\bm{\gamma}_{m}}{N}I_{i}^{2}.$

Integrating both sides and taking $\liminf$ as $t\to\infty$, we obtain

$\liminf_{t\to\infty}S_{i}(t)\geq
S_{i}^{0}+\frac{\bm{\gamma}_{m}}{N}\int_{0}^{\infty}I_{i}^{2}dt>0.$

(ii) Suppose that $H^{-}\cup H^{0}=\emptyset$, that is $H^{+}=\Omega$. Then,
in addition to the set $\mathcal{N}_{1}$, the set
$\mathcal{N}_{2}:=\\{(I^{*}(\bm{r}\circ\bm{\alpha}/(\bm{1}-\bm{r})),I^{*}\bm{\alpha})\\}$
is also a maximal connected invariant subset of $\mathcal{N}$, where
$I^{*}=N/(1+\|\bm{1}/(\bm{1}-\bm{r})\|_{1})$. Hence by the invariance
principle [10, Theorem 4.3.4], $(\bm{S}(t),\bm{I}(t))\to\mathcal{N}_{i}$ for
some $i=1,2$ as $t\to\infty$. Observe that $\mathcal{N}_{2}$ contains a single
point. To complete the proof of the result, it remains to show that
$\bm{I}(t)\nrightarrow\bm{0}$ as $t\to\infty$. To this end, setting
$F_{i}(t)=I_{i}/(S_{i}+I_{i})$ for $t>0$ and $i\in\Omega$, then
$(\bm{S}(t),\bm{I}(t))$ solves the cooperative system

$\begin{cases}\bm{S}^{\prime}=(\bm{\gamma}\circ\bm{I}-(\bm{\beta}-\bm{\gamma})\circ\bm{S})\circ\bm{F},&t>0,\cr\bm{I}^{\prime}=d_{I}\mathcal{L}\bm{I}+((\bm{\beta}-\bm{\gamma})\circ\bm{S}-\bm{\gamma}\circ\bm{I})\circ\bm{F},&t>0.\end{cases}$
(4.10)

Taking

$(\underline{\bm{S}},\underline{\bm{I}})=(l(\bm{\gamma}/(\bm{\beta}-\bm{\gamma}))\circ\bm{\alpha},l\bm{\alpha})$

where
$l:=\min\\{(\bm{S}(1)\circ(\bm{\beta}-\bm{\gamma})/\bm{\gamma}\circ\bm{\alpha})_{m},(\bm{I}(1)/\bm{\alpha})_{m}\\}>0$,
we have that $(\underline{\bm{S}},\underline{\bm{I}})$ solves (4.10) and
$(\underline{\bm{S}},\underline{\bm{I}})\leq(\bm{S}(1),\bm{I}(1))$. Therefore,
by the comparison principle for cooperative systems, we have that
$(\underline{\bm{S}},\underline{\bm{I}})\leq(\bm{S}(t),\bm{I}(t))$ for all
$t\geq 1$. This shows that the omega limit set of $(\bm{S}(t),\bm{I}(t))$
doesn’t intersect with $\mathcal{N}_{1}$. Therefore,
$(\bm{S}(t),\bm{I}(t))\to\mathcal{N}_{2}$ as $t\to\infty$.

∎

######  Remark 4.4.

By Theorem 4.3, if $H^{-}\cup H^{0}$ is nonempty, the disease can be
controlled by limiting the movement of susceptible people. This agrees with
previous studies for the model with standard incidence mechanism in [1, 4],
which are based on the asymptotic profiles of the EE as $d_{I}\to 0$. If
$H^{-}\cup H^{0}$ is empty, then the disease cannot be controlled by limiting
the movement of susceptible people.

## References

  * [1] L. J. S. Allen, B. M. Bolker, Y. Lou, and A. L. Nevai.  Asymptotic profiles of the steady states for an SIS epidemic patch model.  SIAM Journal on Applied Mathematics, 67(5):1283–1309, 2007.
  * [2] L. J. S. Allen, B. M. Bolker, Y. Lou, and A. L. Nevai.  Asymptotic profiles of the steady states for an SIS epidemic reaction-diffusion model.  Discrete Contin. Dyn. Syst., 21(1):1–20, 2008.
  * [3] K. Castellano and R. B. Salako.  On the effect of lowering population’s movement to control the spread of an infectious disease.  Journal of Differential Equations, 316:1–27, 2022.
  * [4] S. Chen, J. Shi, Z. Shuai, and Y. Wu.  Asymptotic profiles of the steady states for an SIS epidemic patch model with asymmetric connectivity matrix.  Journal of Mathematical Biology, 80(7):2327–2361, 2020.
  * [5] J. T. Doumatè, T. B. Issa, and R. B. Salako.  Competition-exclusion and coexistence in a two-strain SIS epidemic model in patchy environments.  arXiv:2308.10348.
  * [6] Daozhou Gao.  How does dispersal affect the infection size?  SIAM Journal on Applied Mathematics, 80(5):2144–2169, 2020.
  * [7] Daozhou Gao and Yuan Lou.  Impact of state-dependent dispersal on disease prevalence.  Journal of Nonlinear Science, 31(5):73, 2021.
  * [8] Hongbin Guo, Michael Y Li, and Zhisheng Shuai.  Global stability of the endemic equilibrium of multigroup sir epidemic models.  Canadian applied mathematics quarterly, 14(3):259–284, 2006.
  * [9] Hongbin Guo, Michael Y Li, and Zhisheng Shuai.  Global dynamics of a general class of multistage models for infectious diseases.  SIAM Journal on applied mathematics, 72(1):261–279, 2012.
  * [10] D. Henry.  Geometric theory of semilinear parabolic equations.  1981\.
  * [11] H. Li and R. Peng.  Dynamics and asymptotic profiles of endemic equilibrium for SIS epidemic patch models.  Journal of mathematical biology, 79:1279–1317, 2019.
  * [12] H. Li and R. Peng.  An SIS epidemic model with mass action infection mechanism in a patchy environment.  Studies in Applied Mathematics, 150(3):650–704, 2023.
  * [13] H. Li, R. Peng, and Z. Wang.  On a diffusive susceptible-infected-susceptible epidemic model with mass action mechanism and birth-death effect: analysis, simulations, and comparison with other mechanisms.  SIAM J. Appl. Math., 78(4):2129–2153, 2018.
  * [14] M. Y. Li and Z. Shuai.  Global-stability problem for coupled systems of differential equations on networks.  J. Differential Equations, 248(1):1–20, 2010.
  * [15] R. Salako and Y. Wu.  On degenerate reaction-diffusion epidemic models with mass action or standard incidence mechanism.  Submitted.
  * [16] R. Salako and Y. Wu.  On the dynamics of an epidemic patch model with mass-action transmission mechanism and asymmetric dispersal patterns.  Submitted.
  * [17] Hongying Shu, Dejun Fan, and Junjie Wei.  Global stability of multi-group seir epidemic models with distributed delays and nonlinear transmission.  Nonlinear Analysis: Real World Applications, 13(4):1581–1592, 2012\.
  * [18] Zhisheng Shuai and Pauline van den Driessche.  Global stability of infectious disease models using lyapunov functions.  SIAM Journal on Applied Mathematics, 73(4):1513–1532, 2013.
  * [19] Y. Wu and X. Zou.  Asymptotic profiles of steady states for a diffusive SIS epidemic model with mass action infection mechanism.  J. Differential Equations, 261(8):4424–4447, 2016.
  * [20] Chunmei Zhang, Wenxue Li, and Ke Wang.  Graph-theoretic method on exponential synchronization of stochastic coupled networks with markovian switching.  Nonlinear Analysis: Hybrid Systems, 15:37–51, 2015.