11institutetext: Daewa Kim 22institutetext: Department of Mathematics, West
Virginia University, 94 Beechurst Ave, Morgantown, WV 26506 22email:
<EMAIL_ADDRESS>33institutetext: Annalisa Quaini 44institutetext:
Department of Mathematics, University of Houston, 3551 Cullen Blvd, Houston TX
77204 44email<EMAIL_ADDRESS>
# A 2D kinetic model for crowd dynamics with disease contagion

Daewa Kim and Annalisa Quaini

###### Abstract

We focus on the modeling and simulation of an infectious disease spreading in
a medium size population occupying a confined environment, such as an airport
terminal, for short periods of time. Because of the size of the crowd and
venue, we opt for a kinetic type model. The paper is divided into two parts.
In the first part, we adopt the simplifying assumption that people’s walking
speed and direction are given. The resulting kinetic model features a variable
that denotes the level of exposure to people spreading the disease, a
parameter describing the contagion interaction strength, and a kernel function
that is a decreasing function of the distance between a person and a spreading
individual. Such model is tested on problems involving a small crowd in a
square walkable domain. In the second part, ideas from the simplified model
are used to incorporate disease spreading in a kinetic theory approach for
crowd dynamics, i.e. the walking speed and direction result from interaction
with other people and the venue.

##  1 Introduction

The goal of this paper is to extend to 2D a 1D model presented in
kim_quaini2021  to study the early stage of an infectious disease spreading in
an intermediate size population occupying a confined environment, e.g. a ER
wait room, for a short period of time (minutes or hours). In this context,
classical epidemiological models fail as they rely on averaged large
population behaviors over a long time span (weeks or months). We focus on a
kinetic theory approach Bellomo2017_book ; Bellomo2021  for a disease that it
spreads with close proximity of individuals, like, e.g. measles, influenza or
COVID-19. In the first part of the paper, we present a model for the disease
spreading based on the simplifying assumption that people’s walking speed and
direction are given. In the second part, ideas from the simplified model are
used to incorporate disease spreading in a kinetic theory approach for crowd
dynamics.

The COVID-19 pandemic has motivated a large body of literature in mathematical
epidemiology. Several possible models have been investigated, including
individual-based Markov models (see, e.g., Flandoli2021 ), Susceptible,
Exposed, Infectious or Recovered (SEIR) models (see, e.g., Hethcote2000 ), and
networks of nodes where each mode dynamics has a SEIR structure (see, e.g.,
Boscheri2021 ; Gatto2020 ). Some works (e.g., Mpeshe2021 ) model the evolution
of the COVID-19 epidemics in combination with other phenomena, such as fear.
Other works focus a multiscale modeling approach. A notable example that
features an interdisciplinary framework including applied mathematics,
immunology, economics and virology is Bellomo_preprint . Similarly, modeling
and simulation of disease spreading in pedestrian crowds has recently become a
topic of increasing relevance. In Derjany2020 , contact evolution is combined
with a stochastic infection spread model to simulate disease spreading in
queuing pedestrians. Agent-based numerical simulations of pedestrian dynamics
are used in Harweg2021  to assess the behavior of pedestrians in public places
in the context of contact transmission of infectious diseases and gather
insights about exposure times and the overall effectiveness of distancing
measures. Finally, we mention the combination of a microscopic force-based
model with a contact tracking method in Rathinakumar_quaini2020  to simulate
the initial spreading of a highly infectious airborne disease in a confined
environment.

The model that we propose for disease spreading in a walking crowd takes
inspiration from the work on emotional contagion (i.e., spreading of fear or
panic) in Ref. Bertozzi2014ContagionSI ; Bertozzi2015 . In order to focus on
the spreading mechanism, we first introduce a simplifying assumption: the
walking speed and direction are given. The key features of the resulting
simplified model are:

  * -

a variable that denotes the level of exposure to people spreading the disease,
with the underlying idea that the more a person is exposed the more likely
they are to get infected;

  * -

a model parameter that describes the contagion interaction strength;

  * -

a kernel function that is a decreasing function of the distance between a
person and a spreading individual.

We show preliminary results for a problem involving a small crowd in a square
walkable domain with no obstacles or walls.

In the second part of the paper, the simplifying assumption is removed: the
walking speed and direction result from a kinetic model of crowd dynamics,
which incorporates a term to account for disease spreading. Kinetic models are
derived from observing the system at the mesoscale (see, e.g., Agnelli2015 ;
Bellomo2011383 ; Bellomo2017_book ; Bellomo2013_new ; Bellomo2015_new ;
Bellomo2016_new ; Bellomo2019_new ), i.e. the intermediate scale between the
macroscopic one (see, e.g., Cristiani2014_book ; Hughes2003 ; 5773492 ) and
the microscopic one (see, e.g., ASANO2010842 ; Chraibi2019 ; DAI20132202 ;
Helbing1995 ). The framework for mesoscale models comes from the kinetic
theory of gases, with the main difference that interactions of “active”
particles are irreversible, non-conservative and, in some cases, non-local and
nonlinearly additive Aristov_2019 . In general, one derives the kinetic model
from interactions at the microscopic scale and the same modeling principles
can lead to macroscopic (or hydrodynamic) models. See, e.g., Ref. Bellomo2020
. We remark that the approach to model crowd dynamics with disease contagion
presented here is different from the one in Ref. kim_quaini2020 , where two
models (one for pedestrian dynamics and one for disease spreading) are
coupled.

The particular kinetic model for crowd dynamics that we extend to account for
disease spreading was first presented in kim_quaini . Based on earlier works
Agnelli2015 ; Bellomo2015_new , the model’s main features are

  * -

discrete walking directions;

  * -

interactions modeled using tools of stochastic games;

  * -

heuristic, deterministic modeling of the walking speed corroborated by
experimental evidence Schadschneider2011545 .

In kim_quaini , the model has been shown to compare favorably against
experimental data related to egression from a facility of a medium-sized group
of people (40 to 138 pedestrians) Kemloh . In addition, we demonstrated that
realistic scenarios, such as passengers moving through one terminal of Hobby
Airport in Houston (USA), can be reproduced kim_quaini2021 .

The paper is organized as follows. We introduce our simplified contagion model
in Sec. 2 and its numerical discretization in Sec. 3. Numerical results are
shown in Sec. 4. Sec. 5 presents the kinetic theory approach for crowd
dynamics extended to incorporate disease spreading. Conclusions are drawn in
Sec. 6.

##  2 A simplified two-dimensional kinetic model

We start from an agent-based model. At the microscopic level we consider a
group of $N$ people divided between healthy or not spreading the disease yet
($N_{h}$) and actively spreading ($N_{s}$), with $N=N_{h}+N_{s}$. If person
$n$ belongs to the group of healthy (or not spreading) people, we denote with
$q_{n}\in[0,1)$ their level of exposure to people spreading the disease, with
the underlying idea that the more a person is exposed the more likely they are
to get infected. If person $n$ is actively spreading the disease, then
$q_{n}=1$ and this value stays constant throughout the entire simulation time.

Let $[t_{0},t_{f}]$ be a time interval of interest. Let
${\mbox{\boldmath{$x$}}}_{n}(t)=(x_{n}(t),y_{n}(t))^{T}$ and $v_{n}(t)$ denote
the position and speed of person $n$ for $t\in[t_{0},t_{f}]$, with initial
position ${\mbox{\boldmath{$x$}}}_{n}(t_{0})$ and initial speed $v_{n}(t_{0})$
being given. The microscopic model reads for $n=1,2,3,\dots,N$:

$\displaystyle\frac{d{\mbox{\boldmath{$x$}}}_{n}}{dt}=v_{n}(\cos\theta,\sin\theta)^{T},~{}\frac{dq_{n}}{dt}=\gamma\max\\{(q_{n}^{*}-q_{n}),0\\},~{}q_{n}^{*}=\frac{\sum_{m=1}^{N}\kappa_{n,m}q_{m}}{\sum_{m=1}^{N_{s}}\kappa_{n,m}},$
(1)

where the walking speed $v_{n}$ and walking direction $\theta$ are assumed to
be given. In model (1), $q_{n}^{\ast}$ corresponds to a weighted average
“level of sickness” surrounding person $n$, with $\kappa_{n,m}$ serving as the
weight in the average. We define $\kappa_{n,m}$ as follow

$\kappa_{n,m}=\left\\{\begin{split}&\kappa(|{\mbox{\boldmath{$x$}}}_{n}-{\mbox{\boldmath{$x$}}}_{m}|)=\dfrac{R}{(|{\mbox{\boldmath{$x$}}}_{n}-{\mbox{\boldmath{$x$}}}_{m}|^{2}+R^{2})\pi}\quad\text{if
}m\text{ is spreading}\\\ &0\hskip
176.407pt\text{otherwise}\end{split}\right.$ (2)

Notice that if person $m$ is spreading the disease, the interaction kernel is
a decreasing function of mutual distance between two people and is
parametrized by an interaction distance $R$. The value of $R$ is set so that
the value of $\kappa_{n,m}$ is “small” at about 6 ft or 2 m. More details are
given in Sec. 4. The meaning of parameter $\gamma$ in (1) is contagion
strength: for $\gamma=0$ there is no contagion, while for $\gamma\neq 0$ the
contagion occurs and it is faster for larger the values of $\gamma$. Note that
obviously the level of exposure can only increase over time. In addition, the
second equation in (1) ensures that the people spreading the disease will
constantly have $q_{n}=1$ in time.

In order to derive a kinetic model from the agent-based model (1), we
introduce the empirical distribution:

$h^{N}=\frac{1}{N}\sum_{n=1}^{N}\delta({\mbox{\boldmath{$x$}}}-{\mbox{\boldmath{$x$}}}_{n}(t))\delta(q-q_{n}(t)),\\\
$

where $\delta$ is the Dirac delta measure. We assume that all people remain in
a fixed compact domain $\Omega$ for the entire time interval under
consideration:
$({\mbox{\boldmath{$x$}}}_{n}(t),q_{n}(t))\in\Omega\subset\mathbb{R}^{3}$, for
$n=1,2,3,\dots,N$ and $t\in[t_{0},t_{f}]$. Let $\mathcal{P}(\mathbb{R}^{3}$)
be the space of probability measures on $\mathbb{R}^{3}$. The sequence
$\\{h^{N}\\}$ is relatively compact in the weak∗ sense (see, e.g.,
Billingsley1999 ). Therefore, there exists a subsequence $\\{h^{N_{k}}\\}_{k}$
such that $h^{N_{k}}$ converges to $h$ with weak∗-convergence in
$\mathcal{P}(\mathbb{R}^{3}$) and pointwise convergence in time as
$k\rightarrow\infty$.

Let $\psi\in C_{0}^{1}(\mathbb{R}^{3})$ be a test function. We have

$\displaystyle\frac{d}{dt}\langle
h^{N},\psi\rangle_{{\mbox{\boldmath{$x$}}},q}$
$\displaystyle=\frac{d}{dt}\biggl{<}\frac{1}{N}\sum_{n=1}^{N}\delta({\mbox{\boldmath{$x$}}}-{\mbox{\boldmath{$x$}}}_{n}(t))\delta(q-q_{n}(t)),\psi\biggr{>}_{{\mbox{\boldmath{$x$}}},q}$
$\displaystyle=\frac{d}{dt}\frac{1}{N}\sum_{n=1}^{N}\psi({\mbox{\boldmath{$x$}}}_{n}(t),q_{n}(t))$
$\displaystyle=\frac{1}{N}\sum_{n=1}^{N}\left(\psi_{x}v_{n}\cos\theta+\psi_{y}v_{n}\sin\theta+\psi_{q}\gamma\max\\{(q_{n}^{*}-q_{n}),0\\}\right)$
$\displaystyle=\langle
h^{N},\psi_{x}v\cos\theta\rangle_{{\mbox{\boldmath{$x$}}},q}+\langle
h^{N},\psi_{y}v\sin\theta\rangle_{{\mbox{\boldmath{$x$}}},q}$
$\displaystyle\quad\quad+\frac{\gamma}{N}\sum_{n=1}^{N}\psi_{q}\max\left\\{\left(\frac{\sum_{m=1}^{N}\kappa_{n,m}q_{n}}{\sum_{m=1}^{N}\kappa_{n,m}}-q_{n}\right),0\right\\},$
(3)

where $\langle\cdot\rangle_{{\mbox{\boldmath{$x$}}},q}$ means integration
against both $x$ and $q$.

Let us define

$\rho({\mbox{\boldmath{$x$}}})=\frac{1}{N}\sum_{n=1}^{N}\delta({\mbox{\boldmath{$x$}}}-{\mbox{\boldmath{$x$}}}_{n})\\\
$

and

$m({\mbox{\boldmath{$x$}}})=\biggl{<}q,\frac{1}{N}\sum_{m=1}^{N}\delta({\mbox{\boldmath{$x$}}}-{\mbox{\boldmath{$x$}}}_{m})\delta(q-q_{m})\biggr{>}_{{\mbox{\boldmath{$x$}}},q}=\frac{1}{N}\sum_{m=1}^{N}\delta({\mbox{\boldmath{$x$}}}-{\mbox{\boldmath{$x$}}}_{m})q_{m},\\\
$

We have

$\displaystyle\frac{1}{N}\sum_{m=1}^{N}\kappa(|{\mbox{\boldmath{$x$}}}_{n}-{\mbox{\boldmath{$x$}}}_{m}|)$
$\displaystyle=\biggl{<}\kappa(|{\mbox{\boldmath{$x$}}}_{n}-\tilde{{\mbox{\boldmath{$x$}}}}|),\frac{1}{N}\sum_{m=1}^{N}\delta(\tilde{{\mbox{\boldmath{$x$}}}}-{\mbox{\boldmath{$x$}}}_{m})\biggr{>}_{{\mbox{\boldmath{$x$}}}}=\kappa\ast\rho({\mbox{\boldmath{$x$}}}_{n}),$
$\displaystyle\frac{1}{N}\sum_{m=1}^{N}\kappa(|{\mbox{\boldmath{$x$}}}_{n}-{\mbox{\boldmath{$x$}}}_{m}|)q_{m}$
$\displaystyle=\biggl{<}\kappa(|{\mbox{\boldmath{$x$}}}_{n}-\tilde{{\mbox{\boldmath{$x$}}}}|),\frac{1}{N}\sum_{m=1}^{N}\delta(\tilde{{\mbox{\boldmath{$x$}}}}-{\mbox{\boldmath{$x$}}}_{m})q_{m}\biggr{>}_{{\mbox{\boldmath{$x$}}}}=\kappa\ast
m({\mbox{\boldmath{$x$}}}_{n}),$

where $\langle\cdot\rangle_{\mbox{\boldmath{$x$}}}$ means integration only in
$x$. Then, we can rewrite eq. (2) as

$\displaystyle\frac{d}{dt}\langle
h^{N},\psi\rangle_{{\mbox{\boldmath{$x$}}},q}$ $\displaystyle=\langle
h^{N},\psi_{x}v\cos\theta\rangle_{{\mbox{\boldmath{$x$}}},q}+\langle
h^{N},\psi_{y}v\sin\theta\rangle_{{\mbox{\boldmath{$x$}}},q}$
$\displaystyle\quad\quad+\gamma\biggl{<}h^{N},\psi_{q}\max\left\\{\frac{\kappa\ast
m_{s}}{\kappa\ast\rho_{s}}-q,0\right\\}\biggr{>}_{{\mbox{\boldmath{$x$}}},q}.$
(4)

Via integration by parts, eq. (4) leads to

$h_{t}^{N}+\nabla\cdot(v(\cos\theta,\sin\theta)^{T}h^{N})+\gamma(\max\\{(q^{\ast}-q),0\\}h^{N})_{q}=0,$
(5)

where $q^{*}$ is the local _average_ sickness level weighted by (2):

$q^{\ast}(t,{\mbox{\boldmath{$x$}}})=\frac{\iint\kappa(|{\mbox{\boldmath{$x$}}}-\overline{{\mbox{\boldmath{$x$}}}}|)h(t,\overline{{\mbox{\boldmath{$x$}}}},q)qdqd\overline{{\mbox{\boldmath{$x$}}}}}{\iint\kappa(|{\mbox{\boldmath{$x$}}}-\overline{{\mbox{\boldmath{$x$}}}}|)h(t,\overline{{\mbox{\boldmath{$x$}}}},q)dqd\overline{{\mbox{\boldmath{$x$}}}}}.$
(6)

Letting $k\rightarrow\infty$, we obtain the limiting kinetic model from the
subsequence $h^{N_{k}}$:

$h_{t}+\nabla\cdot(v(\cos\theta,\sin\theta)^{T}h)+\gamma(\max\\{(q^{\ast}-q),0\\}h)_{q}=0,$
(7)

where $h(t,{\mbox{\boldmath{$x$}}},q)$ is the probability of finding at time
$t$ and position $x$ a person with level of exposure $q$ if $q\in[0,1)$ or a
person spreading the disease if $q=1$.

It is convenient to switch to non-dimensional variables. Let $D$ be the
largest distance a pedestrian can cover in domain $\Omega$ and $V_{M}$ the
largest speed a person can reach. Then, we can define characteristic time
$T=D/V_{M}$. Moreover, let $\rho_{M}$ be the maximum people density per square
meter. The non-dimensional quantities are: position
$\hat{{\mbox{\boldmath{$x$}}}}={\mbox{\boldmath{$x$}}}/D$, walking speed
$\hat{v}=v/V_{M}$, time $\hat{t}=t/T$, and distribution function
$\hat{h}=h/\rho_{M}$. For ease of notation though, we will drop the hat with
the understanding that all quantities are non-dimensional unless otherwise
specified.

##  3 Discretization in space and time

In this section, we present the numerical discretization of model (7).

Let us start from the space discretization. Divide the spatial domain into a
number of cells
$[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}]\times[y_{j-\frac{1}{2}},y_{j+\frac{1}{2}}]$
of length $\Delta x$ and $\Delta y$, respectively. The discrete mesh points
$x_{i}$ and $y_{j}$ are given by

$\displaystyle x_{i}=i\Delta x,\quad x_{i+1/2}=x_{i}+\frac{\Delta
x}{2}=\Big{(}i+\frac{1}{2}\Big{)}\Delta x,$ $\displaystyle y_{j}=j\Delta
y,\quad y_{j+1/2}=y_{j}+\frac{\Delta y}{2}=\Big{(}j+\frac{1}{2}\Big{)}\Delta
y,$

for $i=0,1,\dots,N_{x}$ and $j=0,1,\dots,N_{y}$, The contagion level domain is
partitioned into subdomains $[q_{k-\frac{1}{2}},q_{k+\frac{1}{2}}]$ of length
of $\Delta q$ with $k\in 1,2,\dots,N_{q}$.

$\displaystyle q_{k}=k\Delta q,\quad q_{k+1/2}=q_{k}+\frac{\Delta
q}{2}=\Big{(}k+\frac{1}{2}\Big{)}\Delta q.$

Here, $N_{x}$, $N_{y}$ and $N_{q}$ are the total number of points in $x-$,
$y-$ and $q-$ directions, respectively.

Let us denote $h_{i,j,k}=h(t,x_{i},y_{j},q_{k})$. The average level of
sickness is computed using a midpoint rule for the integrals in (6):

$\displaystyle q^{\ast}(t,x_{i},y_{j})\approx q^{\ast}_{i,j}$
$\displaystyle=\frac{\sum_{k}\sum_{\overline{j}}\sum_{\overline{i}}\kappa_{\overline{i},\overline{j}}h_{\overline{i},\overline{j},k}q_{k}\Delta
q\Delta x\Delta
y}{\sum_{k}\sum_{\overline{j}}\sum_{\overline{i}}\kappa_{\overline{i},\overline{j}}h_{\overline{i},\overline{j},k}\Delta
q\Delta x\Delta y},$ (8)

where if at $(x_{\overline{i}},y_{\overline{j}})$ there is nonzero probability
of finding a person that is spreading the disease:

$\kappa_{\overline{i},\overline{j}}=\frac{R}{((x_{i}-x_{\overline{i}})^{2}+(y_{j}-y_{\overline{j}})^{2}+R^{2})\pi},$

otherwise $\kappa_{\overline{i},\overline{j}}=0$.

We start with a first-order semi-discrete upwind scheme for problem (7), which
reads:

$\displaystyle\partial_{t}h_{i,j,k}$
$\displaystyle+\frac{{\eta}_{i,j,k}-{\eta}_{i-1,j,k}}{\Delta
x}+\frac{{\phi}_{i,j,k}-{\phi}_{i,j-1,k}}{\Delta
y}+\gamma\frac{\xi_{i,j,k+\frac{1}{2}}-\xi_{i,j,k-\frac{1}{2}}}{\Delta q}=0,$
(9)

where

$\displaystyle\eta_{i,j,k}$
$\displaystyle=v\cos\theta\,h_{i,j,k},\quad\phi_{i,j,k}=v\sin\theta\,h_{i,j,k},$
(10) $\displaystyle\xi_{i,j,k+\frac{1}{2}}$
$\displaystyle=\max\left\\{q_{i,j}^{\ast}-q_{k+\frac{1}{2}},0\right\\}h_{i,j,k}.$
(11)

Next, let us discretize in time. Let $\Delta t\in\mathbb{R}$,
$t^{n}=t_{0}+n\Delta t$, with $n=0,...,N_{T}$ and $t_{f}=t_{0}+N_{T}\Delta t$,
where $t_{f}$ is the end of the time interval under consideration. Moreover,
we denote by $y^{n}$ the approximation of a generic quantity $y$ at the time
$t^{n}$. For the time discretization of problem (9), we use the forward Euler
scheme:

$\displaystyle h_{i,j,k}^{l+1}=h_{i,j,k}^{l}-\Delta
t\left(\dfrac{{\eta}^{l}_{i,j,k}-{\eta}^{l}_{i-1,j,k}}{\Delta
x}+\dfrac{{\phi}^{l}_{i,j,k}-{\phi}^{l}_{i,j-1,k}}{\Delta
y}+\gamma\dfrac{\xi^{l}_{i,j,k+\frac{1}{2}}-\xi^{l}_{i,j,k-\frac{1}{2}}}{\Delta
q}\right).$ (12)

To construct a scheme that is of second order, we add a flux limiter. Let
$\varphi$ be a slope limiter function. for example, one could choose the Van
Leer function:

$\varphi(\theta)=\frac{|\theta|+\theta}{1+|\theta|}.$ (13)

The space discretized eq. (7) now reads:

$\displaystyle\partial_{t}h_{i,j,k}$
$\displaystyle+\frac{\overline{\eta}_{i,j,k}-\overline{\eta}_{i-1,j,k}}{\Delta
x}+\frac{\overline{\phi}_{i,j,k}-\overline{\phi}_{i,j-1,k}}{\Delta y}$
$\displaystyle+\gamma\frac{\xi_{i,j,k+\frac{1}{2}}-\xi_{i,j,k-\frac{1}{2}}}{\Delta
q}+\gamma\frac{C_{i,j,k+\frac{1}{2}}-C_{i,j,k-\frac{1}{2}}}{\Delta q}=0,$ (14)

where

$\displaystyle\overline{\eta}_{i,j,k}$
$\displaystyle=\eta_{i,j,k}+\frac{\eta_{i+1,j,k}-\eta_{i,j,k}}{2}\varphi\left(\frac{\eta_{i,j,k}-\eta_{i-1,j,k}}{\eta_{i+1,j,k}-\eta_{i,j,k}}\right),$
(15) $\displaystyle\overline{\phi}_{i,j,k}$
$\displaystyle=\phi_{i,j,k}+\frac{\phi_{i,j+1,k}-\phi_{i,j,k}}{2}\varphi\left(\frac{\phi_{i,j,k}-\phi_{i,j-1,k}}{\phi_{i,j+1,k}-\phi_{i,j,k}}\right),$
(16) $\displaystyle q_{k+\frac{1}{2}}$
$\displaystyle=\frac{q_{k}+q_{k+1}}{2},~{}s_{i,j,k+\frac{1}{2}}=q^{\ast}_{i,j}-q_{k+\frac{1}{2}}$
$\displaystyle W_{i,j,k-\frac{1}{2}}$
$\displaystyle=h_{i,j,k}-h_{i,j,k-\frac{1}{2}}$ $\displaystyle
C_{i,j,k+\frac{1}{2}}$
$\displaystyle=\dfrac{1}{2}\left|s_{i,j,k+\frac{1}{2}}\right|\left(1-\dfrac{\Delta
t}{\Delta
q}\left|s_{i,j,k+\frac{1}{2}}\right|\right)W_{i,j,k-\frac{1}{2}}\varphi\left(\dfrac{W_{i,j,\textbf{k}-\frac{1}{2}}}{W_{i,j,k-\frac{1}{2}}}\right).$
(17)

In eq. (17), the subscript k is k-1 if $s_{i,j,k-\frac{1}{2}}>0$ and k+1 if
$s_{i,j,k-\frac{1}{2}}<0$. Recall that $\eta_{i,j,k}$ in (15) and
$\phi_{i,j,k}$ in (16) are defined in (10), while $\xi_{i,j,k+\frac{1}{2}}$ is
defined in (11). Thanks to the flux limiter, scheme (3) is of second-order
scheme in velocity.

For the time discretization of problem (3), we adopt again the forward Euler
scheme

$\displaystyle h_{i,j,k}^{l+1}=h_{i,j,k}^{l}-\Delta t\Big{(}$
$\displaystyle\dfrac{\overline{\eta}^{l}_{i,j,k}-\overline{\eta}^{l}_{i-1,j,k}}{\Delta
x}+\dfrac{\overline{\phi}^{l}_{i,j,k}-\overline{\phi}^{l}_{i,j-1,k}}{\Delta
y}$ $\displaystyle+$
$\displaystyle\gamma\dfrac{\xi^{l}_{i,j,k+\frac{1}{2}}-\xi^{l}_{i,j,k-\frac{1}{2}}}{\Delta
q}+\gamma\dfrac{C^{l}_{i,j,k+\frac{1}{2}}-C^{l}_{i,j,k-\frac{1}{2}}}{\Delta
q}\Big{)}$ (18)

The time step $\Delta t$ for scheme (12) or (3) is set as

$\displaystyle\Delta t=\frac{1}{2}\min\Bigg{\\{}\frac{\Delta
x}{\max_{k}q_{k}},\frac{\Delta y}{\max_{k}q_{k}},\frac{\Delta
q}{2\gamma\max_{k}q_{k}}\Bigg{\\}}$ (19)

to satisfy the Courant-Friedrichs-Lewy (CFL) condition.

##  4 Numerical results

We assess the approach presented in Sec. 3 through a series of tests. For all
the problems, the computational domain in the $xy$-plane is
$[0,10]\times[0,10]$ m2 before non-dimensionalization, while $q\in[0,1]$. The
dimensionless quantities are obtained by using the following reference
quantities: $D=10\sqrt{2}$ m, $V_{M}=1$ m/s, $T=10\sqrt{2}$ s, $\rho_{M}=7$
people/m2. In all the tests, we take the initial density to be constant in
space and equal to $\rho_{M}$. So, the group occupying the walkable domain
consists of 700 people.

We set the interaction distance $R=1$ m since this choice makes the value of
the kernel function relatively small at a distance of 2 m (or about 6 ft). See
Fig. 1.

\begin{overpic}[width=238.49231pt,grid=false]{./kernel_function_0621.pdf}
\end{overpic} Figure 1: Kernel function vs the distance between people for
interaction radius $R=1,1.5,2$.

To discretize the walkable domain, we take $\Delta x=\Delta y=0.1$ m, while we
set $\Delta q=0.01$. We will consider two values for the contagion strength:
$\gamma=50$ and $\gamma=10$. The time step is set to $\Delta t=5\cdot 10^{-5}$
s for $\gamma=50$ and to $\Delta t=25\cdot 10^{-5}$ s for $\gamma=10$
according to (19). First, we keep the group of people still (i.e., $v=0$) to
make sure that the level of exposure evolves as expected. Then, in a second
set of tests, we let people walk and observe how the motion affects the
spreading of the disease. We run each simulation for $t\in(0,10]$ s.

###  4.1 Tests with $\boldsymbol{v=0}$

We adapt two initial conditions from kim_quaini2021 :

  * -

IC1: people that are certainly spreading (i.e., $q=1$) are located outside the
circle of radius 1 m centered at $(5,5)$ m, while within that circle we place
people that have certainly not been exposed (i.e, $q=0$). See Fig. 2 (left).

  * -

IC2: people that are certainly spreading (i.e., $q=1$) are located outside the
circle of radius 2 m centered at $(5,5)$ m, while the rest of the people
located within that circle have certainly not been exposed (i.e, $q=0$). See
Fig. 2 (right).

The difference between the two boundary conditions is that all the healthy
people in IC1 are exposed to spreading people, while some healthy people in
IC2 are not exposed to spreading people. Thus, we expect that the level of
exposure for the people that are centrally located rises faster for IC1 than
for IC2.

\begin{overpic}[width=182.1196pt,grid=false,tics=10]{./test1_contagion_time_0.pdf}
\end{overpic}\begin{overpic}[width=182.1196pt,grid=false,tics=10]{./test2_contagion_time_0.pdf}
\end{overpic}

Figure 2: Initial conditions IC1 (left) and IC2 (right).

Fig. 3 shows the evolution of the distribution density $h$ for initial
conditions IC1 and IC2 with $\gamma=50$, together with a view of the results
on a section of the 3D domain. We observe that the level of exposure of the
central group of healthy people in IC1 increases quickly. In particular, we
see that the the shape from $q^{*}$ quickly moves away from being close to a
sharp discontinuity and gets closer to a paraboloid, since the level of
exposure increases faster for the people closer to the circle separating
healthy people from spreading people. As expected, the rise in the level of
exposure is much slower for the simulation with initial condition IC2. This is
evident when comparing the second and fourth row in Fig. 3. Moreover, notice
that the healthy people close to the center of the circle in IC2 get very
little exposure to spreading people.

\begin{overpic}[width=134.42113pt,grid=false,tics=10]{./test1_contagion_time_30000}
\put(-18.0,45.0){IC1}. \end{overpic}

\begin{overpic}[width=134.42113pt,grid=false,tics=10]{./test1_contagion_time_150000}
\end{overpic}

\begin{overpic}[width=134.42113pt,grid=false,tics=10]{./test1_contagion_time_300000}
\end{overpic}

\begin{overpic}[width=134.42113pt,grid=false,tics=10]{./test1_slice_time_30000}
\put(-18.0,45.0){IC1} \end{overpic}

\begin{overpic}[width=134.42113pt,grid=false,tics=10]{./test1_slice_time_150000}
\end{overpic}

\begin{overpic}[width=134.42113pt,grid=false,tics=10]{./test1_slice_time_300000}
\end{overpic}

\begin{overpic}[width=134.42113pt,grid=false,tics=10]{./test2_contagion_time_30000}
\put(-18.0,45.0){IC2} \end{overpic}

\begin{overpic}[width=134.42113pt,grid=false,tics=10]{./test2_contagion_time_150000}
\end{overpic}

\begin{overpic}[width=134.42113pt,grid=false,tics=10]{./test2_contagion_time_300000}
\end{overpic}

\begin{overpic}[width=134.42113pt,grid=false,tics=10]{./test2_slice_time_30000}
\put(-18.0,45.0){IC2} \end{overpic}

\begin{overpic}[width=134.42113pt,grid=false,tics=10]{./test2_slice_time_150000}
\end{overpic}

\begin{overpic}[width=134.42113pt,grid=false,tics=10]{./test2_slice_time_300000}
\end{overpic}

Figure 3: Tests with $v=0$: evolution of the distribution density $h$ for
initial condition IC1 (first row) and corresponding results on section $y=0.5$
(second row), evolution of $h$ for initial condition IC2 (third row) and
corresponding results on section $y=0.5$ (fourth row). Note that in all the
subfigures time is dimensional while space is non-dimensional. In both cases,
we set $\gamma=50$. The white dashed line in the images on row two and four
represents $q^{*}$. The legend for all the images is the same as in Fig. 2.

Next, we consider initial condition IC1 and compare the evolution of the
distribution density $h$ for $\gamma=50$ and $\gamma=10$. See Fig. 4. We see
that the level of exposure of the central group of healthy people increases
faster for larger values of $\gamma$, as expected. Parameter $\gamma$ plays a
key role in the evolution of the disease spreading according to model (7). If
one was interested in modeling a realistic scenario, such as the spreading of
COVID-19 in a ER waiting room, $\gamma$ would have to be carefully tuned.

\begin{overpic}[width=134.42113pt,grid=false,tics=10]{./test1_contagion_time_30000}
\put(-23.0,45.0){$\gamma=50$} \end{overpic}

\begin{overpic}[width=134.42113pt,grid=false,tics=10]{./test1_contagion_time_150000}
\end{overpic}

\begin{overpic}[width=134.42113pt,grid=false,tics=10]{./test1_contagion_time_300000}
\end{overpic}

\begin{overpic}[width=134.42113pt,grid=false,tics=10]{./test1_slice_time_30000}
\put(-23.0,45.0){$\gamma=50$} \end{overpic}

\begin{overpic}[width=134.42113pt,grid=false,tics=10]{./test1_slice_time_150000}
\end{overpic}

\begin{overpic}[width=134.42113pt,grid=false,tics=10]{./test1_slice_time_300000}
\end{overpic}

\begin{overpic}[width=134.42113pt,grid=false,tics=10]{./test0_contagion_time_6000}
\put(-23.0,45.0){$\gamma=10$} \end{overpic}

\begin{overpic}[width=134.42113pt,grid=false,tics=10]{./test0_contagion_time_30000}
\end{overpic}

\begin{overpic}[width=134.42113pt,grid=false,tics=10]{./test0_contagion_time_60000}
\end{overpic}

\begin{overpic}[width=134.42113pt,grid=false,tics=10]{./test0_slice_time_6000}
\put(-23.0,45.0){$\gamma=10$} \end{overpic}

\begin{overpic}[width=134.42113pt,grid=false,tics=10]{./test0_slice_time_30000}
\end{overpic}

\begin{overpic}[width=134.42113pt,grid=false,tics=10]{./test0_slice_time_60000}
\end{overpic}

Figure 4:  Tests with $v=0$: evolution of the distribution density $h$ for
$\gamma=50$ (first row) and corresponding results on section $y=0.5$ (second
row), evolution of $h$ for $\gamma=10$ (third row) and corresponding results
on section $y=0.5$ (fourth row). In both cases, the initial condition is IC1.
Note that in all the subfigures time is dimensional while space is non-
dimensional. The white dashed line in the images on row two and four
represents $q^{*}$. The legend for all the images is the same as in Fig. 2.

This first set of tests was meant to verify our implementation of method
described in Sec. 3 and to check that the disease spreading term in eq. (7)
(i.e., the third term on the left-hand side) produced the expected outcomes.
Next, we will set people in motion.

###  4.2 Tests with prescribed walking velocity

We assign to all people walking speed $v=1$ m/s and walking direction
$\theta=\pi/4$, as if they were headed to the upper right corner of the domain
in the $xy$-plane. Once the people in the spreading phase of the disease have
left the domain, we assume they cannot spread to the people in the walkable
domain anymore.

We consider again IC1 and IC2 with $\gamma=50$. Fig. 5 shows the evolution of
the distribution density $h$ for both initial conditions. As one would expect,
the motion contributes to lowering the exposure level in the both cases, since
some of the spreading people leave the domain as soon as the motion starts.
This difference is obvious when one compares the second (resp., fourth) row in
Fig. 3 with the second (resp., fourth) row in Fig. 5.

\begin{overpic}[width=138.76157pt,grid=false,tics=10]{./test3_contagion_time_30000}
\put(57.4,46.0){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(0,1){48.5}}
\put(57.4,46.0){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(-1,-4){9.8}}
\put(47.6,6.8){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(0,1){49.7}}
\put(57.1,94.2){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(-1,-4){9.5}}
\put(-14.0,45.0){IC1} \end{overpic}

\begin{overpic}[width=138.76157pt,grid=false,tics=10]{./test3_contagion_time_150000}
\put(57.4,46.0){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(0,1){48.5}}
\put(57.4,46.0){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(-1,-4){9.8}}
\put(47.6,6.8){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(0,1){49.7}}
\put(57.1,94.2){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(-1,-4){9.5}}
\end{overpic}

\begin{overpic}[width=138.76157pt,grid=false,tics=10]{./test3_contagion_time_240000}
\put(57.4,46.0){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(0,1){48.5}}
\put(57.4,46.0){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(-1,-4){9.8}}
\put(47.6,6.8){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(0,1){49.7}}
\put(57.1,94.2){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(-1,-4){9.5}}
\end{overpic}

\begin{overpic}[width=134.42113pt,grid=false,tics=10]{./test3_slice_time_30000.jpg}
\put(-18.0,45.0){IC1} \end{overpic}

\begin{overpic}[width=134.42113pt,grid=false,tics=10]{./test3_slice_time_150000.jpg}
\end{overpic}

\begin{overpic}[width=134.42113pt,grid=false,tics=10]{./test3_slice_time_240000.jpg}
\end{overpic}

\begin{overpic}[width=138.76157pt,grid=false,tics=10]{./test4_contagion_time_30000}
\put(57.4,46.0){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(0,1){48.5}}
\put(57.4,46.0){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(-1,-4){9.8}}
\put(47.6,6.8){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(0,1){49.7}}
\put(57.1,94.2){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(-1,-4){9.5}}
\put(-14.0,45.0){IC2} \end{overpic}

\begin{overpic}[width=138.76157pt,grid=false,tics=10]{./test4_contagion_time_150000}
\put(57.4,46.0){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(0,1){48.5}}
\put(57.4,46.0){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(-1,-4){9.8}}
\put(47.6,6.8){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(0,1){49.7}}
\put(57.1,94.2){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(-1,-4){9.5}}
\end{overpic}

\begin{overpic}[width=138.76157pt,grid=false,tics=10]{./test4_contagion_time_240000}
\put(57.4,46.0){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(0,1){48.5}}
\put(57.4,46.0){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(-1,-4){9.8}}
\put(47.6,6.8){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(0,1){49.7}}
\put(57.1,94.2){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(-1,-4){9.5}}
\end{overpic}

\begin{overpic}[width=134.42113pt,grid=false,tics=10]{./test4_slice_time_30000.jpg}
\put(-18.0,45.0){IC2} \end{overpic}

\begin{overpic}[width=134.42113pt,grid=false,tics=10]{./test4_slice_time_150000.jpg}
\end{overpic}

\begin{overpic}[width=134.42113pt,grid=false,tics=10]{./test4_slice_time_240000.jpg}
\end{overpic}

Figure 5:  Tests with prescribed walking velocity: the first row shows the
evolution of the distribution density $h$ for initial condition IC1 and
highlighted in green is the section whose results are reported in the second
row. The third row shows the evolution of $h$ for initial condition IC2 and
highlighted in green is the section whose results are reported in the fourth
row. Note that in all the subfigures time is dimensional while space is non-
dimensional. For both cases, $\gamma=50$. The white dashed line in the images
on row two and four represents $q^{*}$. The legend for all the images is the
same as in Fig. 2.

Finally, we modify the initial conditions to show that our model can handle
also more complex scenarios featuring uncertainties. IC1 and IC2 are changed
to:

  * -

IC1-bis: people are positioned like in IC1 but the probabilities of finding
people with $q=1$ and $q=0$ is reduced and another value of $q$ is assigned at
a given $(x,y)$, i.e. at every location in the walkable domain the initial
distribution has two bumps in $q$. See Fig. 6 (left).

  * -

IC2-bis: people are positioned like in IC2 but the probabilities of finding
people with $q=1$ and $q=0$ is reduced and, similarly to IC1-bis, another
value of $q$ is assigned at a given $(x,y)$. See Fig. 6 (right).

Fig. 7 shows the evolution of the distribution density $h$ for initial
conditions IC1-bis and IC2-bis with $\gamma=50$.

\begin{overpic}[width=182.1196pt,grid=false,tics=10]{./test5_contagion_time_0.pdf}
\put(48.6,35.6){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(0,1){33.7}}
\put(48.5,69.2){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(-1,-4){6.9}}
\put(41.6,8.3){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(0,1){33.5}}
\put(48.6,35.6){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(-1,-4){6.8}}
\end{overpic}\begin{overpic}[width=182.1196pt,grid=false,tics=10]{./test6_contagion_time_0.pdf}
\put(48.6,35.6){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(0,1){33.7}}
\put(48.5,69.2){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(-1,-4){6.9}}
\put(41.6,8.3){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(0,1){33.5}}
\put(48.6,35.6){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(-1,-4){6.8}}
\end{overpic}

\begin{overpic}[width=182.1196pt,grid=false,tics=10]{./test5_slice_time_0.jpg}
\end{overpic}

\begin{overpic}[width=182.1196pt,grid=false,tics=10]{./test6_slice_time_0.jpg}
\end{overpic}

Figure 6: Tests with prescribed walking velocity: initial conditions IC1-bis
(top left) and IC2-bis (top right). The initial data on the section
highlighted in green in the top panels is reported in the corresponding bottom
panel.

\begin{overpic}[width=138.76157pt,grid=false,tics=10]{./test5_contagion_time_30000}
\put(57.4,46.0){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(0,1){48.5}}
\put(57.4,46.0){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(-1,-4){9.8}}
\put(47.6,6.8){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(0,1){49.7}}
\put(57.1,94.2){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(-1,-4){9.5}}
\put(-18.0,45.0){IC1-bis} \end{overpic}

\begin{overpic}[width=138.76157pt,grid=false,tics=10]{./test5_contagion_time_150000}
\put(57.4,46.0){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(0,1){48.5}}
\put(57.4,46.0){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(-1,-4){9.8}}
\put(47.6,6.8){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(0,1){49.7}}
\put(57.1,94.2){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(-1,-4){9.5}}
\end{overpic}

\begin{overpic}[width=138.76157pt,grid=false,tics=10]{./test5_contagion_time_240000}
\put(57.4,46.0){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(0,1){48.5}}
\put(57.4,46.0){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(-1,-4){9.8}}
\put(47.6,6.8){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(0,1){49.7}}
\put(57.1,94.2){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(-1,-4){9.5}}
\end{overpic}

\begin{overpic}[width=134.42113pt,grid=false,tics=10]{./test5_slice_time_30000.jpg}
\put(-22.0,43.0){IC1-bis} \end{overpic}

\begin{overpic}[width=134.42113pt,grid=false,tics=10]{./test5_slice_time_150000.jpg}
\end{overpic}

\begin{overpic}[width=134.42113pt,grid=false,tics=10]{./test5_slice_time_240000.jpg}
\end{overpic}

\begin{overpic}[width=138.76157pt,grid=false,tics=10]{./test6_contagion_time_30000}
\put(57.4,46.0){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(0,1){48.5}}
\put(57.4,46.0){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(-1,-4){9.8}}
\put(47.6,6.8){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(0,1){49.7}}
\put(57.1,94.2){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(-1,-4){9.5}}
\put(-18.0,45.0){IC2-bis} \end{overpic}

\begin{overpic}[width=138.76157pt,grid=false,tics=10]{./test6_contagion_time_150000}
\put(57.4,46.0){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(0,1){48.5}}
\put(57.4,46.0){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(-1,-4){9.8}}
\put(47.6,6.8){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(0,1){49.7}}
\put(57.1,94.2){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(-1,-4){9.5}}
\end{overpic}

\begin{overpic}[width=138.76157pt,grid=false,tics=10]{./test6_contagion_time_240000}
\put(57.4,46.0){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(0,1){48.5}}
\put(57.4,46.0){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(-1,-4){9.8}}
\put(47.6,6.8){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(0,1){49.7}}
\put(57.1,94.2){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(-1,-4){9.5}}
\end{overpic}

\begin{overpic}[width=134.42113pt,grid=false,tics=10]{./test6_slice_time_30000.jpg}
\put(-22.0,43.0){IC2-bis} \end{overpic}

\begin{overpic}[width=134.42113pt,grid=false,tics=10]{./test6_slice_time_150000.jpg}
\end{overpic}

\begin{overpic}[width=134.42113pt,grid=false,tics=10]{./test6_slice_time_240000.jpg}
\end{overpic}

Figure 7:  Tests with prescribed walking velocity: the first row shows the
evolution of the distribution density $h$ for initial condition IC1-bis and
highlighted in green the section whose results are reported in the second row.
The third row shows the evolution of $h$ for initial condition IC2-bis and
highlighted in green the section whose results are reported in the fourth row.
Note that in all the subfigures time is dimensional while space is non-
dimensional. In both cases, we set $\gamma=50$. The white dashed line in the
images on row two and four represents $q^{*}$. The legend for all the images
is the same as in Fig. 6.

##  5 A more complex 2D kinetic model

In this section, we remove the simplifying assumption used in Sec. 2, i.e.
walking speed $v$ and direction $\theta$ are no longer assumed to be given.
Instead, they result from the interaction with the surrounding people.
Following Agnelli2015 ; kim_quaini ; kim_quaini2021 , we account for the
granular feature of the system (i.e., the fact that the distance between
pedestrians can range from small to large) by discretizing $\theta$. Instead
of being continuous, variable $\theta$ can only take values in the set:

$I_{\theta}=\left\\{\theta_{i}=\frac{i-1}{N_{d}}2\pi:i=1,\dots,N_{d}\right\\},$

where $N_{d}$ is the maximum number of possible directions. As for the walking
speed $v$, we assume it depends on the local level of congestion as
demonstrated by experimental studies. Thus, variable $v$ is treated as a
deterministic variable.

Let ${\mbox{\boldmath{$x$}}}=(x,y)$ denote position. For the above mentioned
assumptions on $\theta$ and $v$, the distribution function can be written as

$f(t,{\mbox{\boldmath{$x$}}},\theta)=\sum_{i=1}^{N_{d}}f^{i}(t,{\mbox{\boldmath{$x$}}})\delta(\theta-\theta_{i}),$
(20)

where $f^{i}(t,{\mbox{\boldmath{$x$}}})$ represents the people that, at time
$t$ and position $x$, move with direction $\theta_{i}$, and $\delta$ denotes
the Dirac delta function. Like for the model in Sec. 2, it is convenient to
switch to non-dimensional variables: position
$\hat{{\mbox{\boldmath{$x$}}}}={\mbox{\boldmath{$x$}}}/D$, walking speed
$\hat{v}=v/V_{M}$, time $\hat{t}=t/T$, and distribution function
$\hat{f}=f/\rho_{M}$. Once again, we will drop the hat for ease of notation
though, with the understanding that all variables from now on are non-
dimensional. Due to the normalization of $f$ and the $f_{i}$ functions, the
dimensionless local density is obtained by summing the distribution functions
over the set of directions:

$\displaystyle\rho(t,{\mbox{\boldmath{$x$}}})=\sum_{i=1}^{N_{d}}f^{i}(t,{\mbox{\boldmath{$x$}}}).$
(21)

The model in Agnelli2015 ; kim_quaini ; kim_quaini2021  features two key
parameters:

  * -

$\alpha\in[0,1]$, which represents the quality of the walkable domain: if
$\alpha=1$ the domain is clear, while if $\alpha=0$ an obstruction is present.
Note that $\alpha$ could be dependent on space (and possibly time) in a
prescribed way.

  * -

$\varepsilon\in[0,1]$, which represents the overall level of stress:
$\varepsilon=0$ indicates no stress, while $\varepsilon=1$ indicates a high
stress stress situation.

Parameter $\alpha$ plays a role in determining the walking speed: where
$\alpha=1$ people can walk at the desired speed (i.e., $V_{M}$), while where
$\alpha=0$ people are forced to slow down or stop. See, e.g., kim_quaini  on
how to define the walking speed $v=v[\rho](t,{\mbox{\boldmath{$x$}}},\alpha)$,
with the square brackets denoting that $v$ depends on $\rho$ in a functional
way. Parameter $\varepsilon$ plays a role in the selection of the walking
direction when interacting with other pedestrians. Indeed, the choice of the
waking direction among the possible directions in $I_{\theta}$ results from
four competing factors:

  1. 1.

The goal to reach a target, like, e.g., the gate of an airport terminal.

  2. 2.

The desire to avoid collisions with the walls and other obstacles.

  3. 3.

The tendency to look for less crowded areas.

  4. 4.

The tendency to follow the stream or “herd”.

Parameter $\varepsilon$ weights between 3 (prevailing when stress is low) and
4 (prevailing when stress is high).

In order to explain how the four factors are modeled, we need to introduce
some terminology. Interactions involve three types of active particles (or
people):

  * -

test particle with state $({\mbox{\boldmath{$x$}}},\theta_{i})$: they are
representative of the whole system;

  * -

candidate particle with state $({\mbox{\boldmath{$x$}}},\theta_{h})$: they can
reach in probability the state of the test particle after individual-based
interactions with the venue (walls, obstacles, or target) or with field
particles;

  * -

field particle with state $({\mbox{\boldmath{$x$}}},\theta_{k})$: their
interactions with candidate particles determines a possible change of state.

Note that if a test particle changes their state (in probability) into that of
the test of the test particle as a result of the interactions with field
people, the test particle loses their state.

We call $\mathcal{A}_{h}(i)$ the probability that a candidate particle with
direction $\theta_{h}$ adjusts their direction to $\theta_{i}$ (the direction
of the test particle) due to the interaction with the venue. The transition
probability that the candidate particle changes their direction to
$\theta_{i}$ (direction of the test particle) in the search for less crowded
areas if their stress level is low or to $\theta_{k}$ (direction of the field
particle) if their stress level is high is denoted by $\mathcal{B}_{hk}(i)$.
The sets of all transition probabilities
$\mathcal{A}=\\{\mathcal{A}_{h}(i)\\}_{h,i=1,\dots,N_{d}}$ and
$\mathcal{B}=\\{\mathcal{B}_{hk}(i)\\}_{h,k,i=1,\dots,N_{d}}$ form the so-
called tables of games that model the game played by active particle
interacting with the venue and other particles, respectively. See, e.g.,
kim_quaini  for more details.

In order to state the mathematical model, we need two more ingredients:

  * -

the interaction rate with the venue $\mu[\rho]$ : it models the frequency of
interactions between candidate particles and the venue and it is assume to
have a functional dependence on $\rho$. In fact, it is easier for pedestrians
to see the walls, possible obstacles, and their target if density is low.

  * -

the interaction rate with other pedestrians $\eta[\rho]$: it defines the
number of binary encounters per unit time. This rate depends on $\rho$ in a
functional way too since if the local density increases the interaction rate
also increases.

For possible ways to set $\mu[\rho]$ and $\eta[\rho]$, see, e.g., kim_quaini .

In order to derive the mathematical model, one has to take the balance of
particles in an elementary volume of the space of microscopic states,
considering the net flow into such volume due to transport and interactions.
Let $\mathcal{J}^{i}[f]$ be the net balance of particles that move with
direction $\theta_{i}$:

$\displaystyle\mathcal{J}^{i}[f](t,{\mbox{\boldmath{$x$}}},q)=\mu[\rho]\left(\sum_{h=1}^{N_{d}}\int_{0}^{1}\mathcal{A}_{h}(i)f^{h}(t,{\mbox{\boldmath{$x$}}},q_{h})dq_{h}-f^{i}(t,{\mbox{\boldmath{$x$}}},q)\right)$
$\displaystyle\quad+\eta[\rho]\left(\sum_{h,k=1}^{N_{d}}\int_{0}^{1}\int_{0}^{1}\mathcal{B}_{hk}(i)f^{h}(t,{\mbox{\boldmath{$x$}}},q_{h})f^{k}(t,{\mbox{\boldmath{$x$}}},q_{k})dq_{h}dq_{k}-f^{i}(t,{\mbox{\boldmath{$x$}}},q)\rho(t,{\mbox{\boldmath{$x$}}})\right)$
(22)

due to interactions with the venue (first term at the left-hand side of (22))
and with the surrounding particles (second term at the left-hand side of
(22)). Then, the model is given by:

$\displaystyle\frac{\partial f^{i}}{\partial t}$
$\displaystyle+\nabla\cdot\left(v(\cos\theta_{i},\sin\theta_{i})^{T}f^{i}(t,{\mbox{\boldmath{$x$}}},q)\right)$
$\displaystyle=\mathcal{J}^{i}[f](t,{\mbox{\boldmath{$x$}}},q)+\gamma\frac{\partial(\max\\{(q^{\ast}-q),0\\}f^{i})}{\partial
q}$ (23)

for $i=1,2,\dots,N_{d}$. Problem (22)-(23) is completed with equation (21) for
the density and a suitable law $v=v[\rho](t,{\mbox{\boldmath{$x$}}},\alpha)$
that relates the walking speed to the density.

While modeling motion in 2 dimensions ($x$ and $y$), eq. (23) is a 3D problem
in variables $x,y$, and $q$. Its numerical approximation requires a carefully
designed scheme to contain the computational cost. Here, we are only going to
present ideas on how to design such a scheme. Numerical results will be
presented in a follow-up paper. Following what we have done in kim_quaini2020b
for emotional contagion, we will split model (23) into subproblems that are
easier to solve using operator splitting (see, e.g., glowinski2003finite ).
For example, with the Lie splitting algorithm one could break problem (23)
into: (i) a pure advection problem that features the contagion term, i.e. the
last term in eq. (23), and (ii) a problem involving the interaction with the
venue and other pedestrians. Then, suitable finite difference scheme will be
applied for the space and time discretization of problems (i) and (ii).

##  6 Conclusions

We presented a way to design a kinetic type model to simulate the early stage
of an infectious disease spreading in an intermediate size population
occupying a confined environment for a short period of time.

In order to focus on the spreading mechanism, in the first part of the paper
we adopted a strong simplifying assumption: people’s walking speed and
direction are given. The disease spreading is modeled using three main
ingredients: an additional variable for the level of exposure to people
spreading the disease, a parameter for the contagion interaction strength, and
a kernel function that is a decreasing function of the distance between a
person and a spreading individual. For such simplified model, we propose a
first-order and a second-order finite difference scheme. We tested our
numerical approach on simple 2D problems and verified that the level of
exposure evolved as expected, even in scenarios featuring uncertainties.

In the second part of the paper, we used ideas from the simplified model to
incorporate disease spreading in a kinetic theory approach for crowd dynamics.
The appealing feature of kinetic (or mesoscopic) models for crowd dynamics is
the flexibility in accounting for multiple interactions (hard to achieve in
microscopic models) and heterogeneous behavior in people (hard to achieve in
macroscopic models). Modeling crowd dynamics and disease contagion in two
spatial dimensions requires the solution of a 3D problem, the additional
variable being the contagion level. The design of an efficient numerical
scheme for such problem will be object of a follow-up paper.

###### Acknowledgements.

This work has been partially supported by NSF through grant DMS-1620384.

## References

  * [1] J. P. Agnelli, F. Colasuonno, and D. Knopoff.  A kinetic theory approach to the dynamics of crowd evacuation from bounded domains.  Mathematical Models and Methods in Applied Sciences, 25(01):109–129, 2015.
  * [2] V.V. Aristov.  Biological systems as nonequilibrium structures described by kinetic methods.  Results in Physics, 13:102232, 2019.
  * [3] M. Asano, T Iryo, and M. Kuwahara.  Microscopic pedestrian simulation model combined with a tactical model for route choice behaviour.  Transportation Research Part C: Emerging Technologies, 18(6):842–855, 2010.
  * [4] B. Aylaj, N. Bellomo, L. Gibelli, and A. Reali.  A unified multiscale vision of behavioral crowds.  Mathematical Models and Methods in Applied Sciences, 30(1):1–22, 2020.
  * [5] N. Bellomo and A. Bellouquid.  On the modeling of crowd dynamics: Looking at the beautiful shapes of swarms.  Networks and Heterogeneous Media, 6(3):383–399, 2011.
  * [6] N. Bellomo, A. Bellouquid, L. Gibelli, and N Outada.  A quest towards a mathematical theory of living systems.  Modeling and Simulation in Science, Engineering and Technology. (Birkhäuser, 2017), 2017.
  * [7] N. Bellomo, A. Bellouquid, and D. Knopoff.  From the microscale to collective crowd dynamics.  SIAM Multiscale Modeling & Simulation, 11(3):943–963, 2013.
  * [8] N. Bellomo, R. Bingham, M.K. Chaplain, G. Dosi, G. Forni, D.A. Knopoff, J. Lowengrub, R. Twarock, and M.E. Virgillito.  A multiscale model of virus pandemic: Heterogeneous interactive entities in a globally connected world.  Mathematical Models and Methods in Applied Sciences, 30(8):1591–1651, 2020.
  * [9] N. Bellomo, D. Burini, G. Dosi, L. Gibelli, D.A. Knopoff, N. Outada, P. Terna, and M.E Virgillito.  What is life? a perspective of the mathematical kinetic theory of active particles.  Mathematical Models and Methods in Applied Sciences, 31(8), 2021\.
  * [10] N. Bellomo and L. Gibelli.  Toward a mathematical theory of behavioral-social dynamics for pedestrian crowds.  Mathematical Models and Methods in Applied Sciences, 25(13):2417–2437, 2015.
  * [11] N. Bellomo and L. Gibelli.  Behavioral crowds: Modeling and Monte Carlo simulations toward validation.  Computers & Fluids, 141:13–21, 2016.
  * [12] N. Bellomo, L. Gibelli, and N. Outada.  On the interplay between behavioral dynamics and social interactions in human crowds.  Kinetic and Related Models, 12(2):397–409, 2019.
  * [13] Andrea L. Bertozzi, Jesús Rosado, Martin B. Short, and Li Wang.  Contagion shocks in one dimension.  Journal Statistical Physics, 158:647–664, 2015.
  * [14] P. Billingsley.  Convergence of probability measures.  Wiley-Interscience publication, 1999.
  * [15] W. Boscheri, G. Dimarco, and L. Pareschi.  Modeling and simulating the spatial spread of an epidemic through multiscale kinetic transport equations.  Mathematical Models and Methods in Applied Sciences, 31:1–39, 2021\.
  * [16] M. Chraibi, A. Tordeux, A. Schadschneider, and A. Seyfried.  Modelling of Pedestrian and Evacuation Dynamics, pages 649–669.  Springer-New York, 2019.
  * [17] E. Cristiani, B. Piccoli, and A. Tosin.  Multiscale Modeling of Pedestrian Dynamics.  Springer-Cham, 2014.
  * [18] J. Dai, X. Li, and L. Liu.  Simulation of pedestrian counter flow through bottlenecks by using an agent-based model.  Physica A: Statistical Mechanics and its Applications, 392(9):2202–2211, 2013.
  * [19] P. Derjany, S. Namilae, D. Liu, and A. Srinivasan.  Multiscale model for the optimal design of pedestrian queues to mitigate infectious disease spread.  PLOS ONE, 15:1–21, 2020.
  * [20] F. Flandoli, E. La Fauci, and M. Riva.  Individual-based markov model of virus diffusion: Comparison with covid-19 incubation period, serial interval and regional time series.  Mathematical Models and Methods in Applied Sciences, 31, 2021.
  * [21] M. Gatto, E. Bertuzzo, L. Mari, S. Miccoli, L. Carraro, R. Casagrandi, and A. Rinaldo.  Spread and dynamics of the covid-19 epidemic in italy: Effects of emergency containment measures.  Proceedings of the National Academy of Sciences, 117(19):10484–10491, 2020.
  * [22] R. Glowinski.  Finite element methods for incompressible viscous flow, in Handbook of numerical analysis, P. G. Ciarlet, J. L. Lions (Eds), volume 9.  North-Holland, 2003.
  * [23] T. Harweg, D. Bachmann, and F. Weichert.  Agent-based simulation of pedestrian dynamics for exposure time estimation in epidemic risk assessment.  Journal of Public Health, 2021.
  * [24] D. Helbing and P. Molnár.  Social force model for pedestrian dynamics.  Physical review. E, 51:4282–4286, 1998.
  * [25] H.W. Hethcote.  The mathematics of infectious diseases.  SIAM Review, 42(4):599–653, 2000.
  * [26] R. L. Hughes.  The flow of human crowds.  Annual Review of Fluid Mechanics, 35(1):169–182, 2003.
  * [27] D. Kim, K. O’Connell, W. Ott, and A. Quaini.  A kinetic theory approach for 2D crowd dynamics with emotional contagion.  Mathematical Models and Methods in Applied Sciences, Accepted., 2020\.  https://arxiv.org/abs/2012.08108.
  * [28] D. Kim and A. Quaini.  A kinetic theory approach to model pedestrian dynamics in bounded domains with obstacles.  Kinetic & Related Models, 12(6):1273–1296, 2019.
  * [29] D. Kim and A. Quaini.  Coupling kinetic theory approaches for pedestrian dynamics and disease contagion in a confined environment.  Mathematical Models and Methods in Applied Sciences, 30(10):1893–1915, 2020.
  * [30] D. Kim and A. Quaini.  A kinetic theory approach to model crowd dynamics with disease contagion.  arXiv:2103.15151, 2021.
  * [31] S.C. Mpeshe and N. Nyerere.  Modeling the dynamics of coronavirus disease pandemic coupled with fear epidemics.  Computational and Mathematical Methods in Medicine, 2021, 2021.
  * [32] K. Rathinakumar and A. Quaini.  A microscopic approach to study the onset of a highly infectious disease spreading.  Mathematical Biosciences, 329:108475, 2020.
  * [33] A. Schadschneider and A. Seyfried.  Empirical results for pedestrian dynamics and their implications for modeling.  Networks and Heterogeneous Media, 6(3):545–560, 2011.
  * [34] A. Shende, M. P. Singh, and P. Kachroo.  Optimization-based feedback control for pedestrian evacuation from an exit corridor.  IEEE Transactions on Intelligent Transportation Systems, 12(4):1167–1176, 2011.
  * [35] A. U. K Wagoum, A. Tordeux, and W. Liao.  Understanding human queuing behaviour at exits: An empirical study.  Royal Society Open Science, 4(1):160896, 2017.
  * [36] L. Wang, M. B. Short, and A. L. Bertozzi.  Efficient numerical methods for multiscale crowd dynamics with emotional contagion.  Mathematical Models and Methods in Applied Sciences, 27(1):205–230, 2017.


######  Lemma 5.6.

There exists $\delta\in(0,1)$ and an open neighborhood
$\mathcal{O}\subseteq\mathbb{T}^{2}$ of $0$ for which

$R(\xi)\leq\frac{1-\delta}{100}P(\xi)$

for $\xi\in\mathcal{O}$.

###### Proof.

In (39), we have written $R(\xi)=R_{1}(\xi)+R_{2}(\xi)$ and noted, on account
of Proposition 4.3, that $R_{2}(\xi)=o(P(\xi)/100)$ as $\xi\to 0$. Thus, to
prove the statement at hand, it suffices to find an open neighborhood
$\mathcal{O}\subseteq\mathbb{T}^{2}$ of $0$ and $\rho\in(0,1)$ for which

$\mathcal{E}(\xi):=\frac{179}{1200}\zeta^{8}+\frac{1}{6}\eta\zeta^{4}-\frac{23}{900}\eta\zeta^{6}-\frac{3}{100}\eta^{2}\zeta^{4}+\frac{47}{150}\eta^{3}\zeta^{2}=100R_{1}(\xi)\leq\rho
P(\xi)$

for all $\xi=(\eta,\zeta)\in\mathcal{O}$. As a first approximation, we set
$\mathcal{O}=\\{\xi=(\eta,\zeta):P(\xi)<1\\}$ and decompose it as
$\mathcal{O}=\mathcal{R}_{1}\cup\mathcal{R}_{2}$ where $\mathcal{R}_{1}$ and
$\mathcal{R}_{2}$ are given by (22) and (23), respectively. In view of our
analysis in Example 12, we observe that, for
$\xi=(\eta,\zeta)\in\mathcal{R}_{1}$, $0\leq\zeta^{2}\leq-\sqrt{2}\eta$ and,
in particular, $\eta$ is non-positive. Therefore,

$\mathcal{E}(\xi)\leq\frac{179}{1200}4\eta^{4}+0+\frac{23}{900}(\sqrt{2})^{3}\eta^{4}+0+0=\left(\frac{179}{300}+\frac{23\sqrt{2}}{450}\right)\eta^{4}\leq\frac{7}{9}P(\xi)$

for $\xi=(\eta,\zeta)\in\mathcal{R}_{1}$. On $\mathcal{R}_{2}$, we have

$\max\left\\{\eta^{4},\left(1-\frac{1}{\sqrt{2}}\right)^{2}\zeta^{4}\right\\}\leq
P(\xi)$

and therefore

$\displaystyle\mathcal{E}(\xi)\leq\frac{179}{1200}(6+4\sqrt{2})^{2}P(\xi)^{2}+\frac{1}{6}(6+4\sqrt{2})P(\xi)^{5/4}$
$\displaystyle\hskip
56.9055pt+\frac{23}{900}(6+4\sqrt{2})^{3/2}P(\xi)^{7/4}+0+\frac{47}{150}(2+\sqrt{2})P(\xi)^{5/4}\leq
150P(\xi)^{5/4}$

where we have used the fact that $P(\xi)<1$ for $\xi\in\mathcal{O}$. Thus, by
further restricting the open set $\mathcal{O}$ so that $150P(\xi)^{1/4}<7/9$,
the preceding estimates ensure that $\mathcal{E}(\xi)\leq\rho P(\xi)$ for all
$\xi\in\mathcal{O}$ where $\rho=7/9$. ∎

###### Proof of Proposition 5.5..

By an appeal to the preceding lemma, let $\mathcal{O}\subseteq\mathbb{T}^{2}$
be an open neighborhood of $0$ for which

$\absolutevalue{\widehat{\psi}(\xi)}=e^{-P(\xi)/100+R(\xi)}\leq e^{-\gamma
P(\xi)}$ (41)

for $\xi\in\mathcal{O}$ where $\gamma=\delta/100>0$. Using similar arguments
to those which appear in the proof of Theorem 5.1, we find that

$\|\varphi^{(n)}\|_{\infty}=\sup_{x\in\mathbb{Z}^{2}}\absolutevalue{\psi^{(n)}(x)}\leq\rho^{n}+\int_{\mathcal{O}}\absolutevalue{\widehat{\psi}(\xi)}^{n}\,d\xi$

for all $n\in\mathbb{N}_{+}$ where $0<\rho<1$. Focusing on the integral above,
we make the change of variables $\xi\mapsto T(n^{-G}\xi)$ to see that

$n^{5/8}\int_{\mathcal{O}}\absolutevalue{\widehat{\psi}(\xi)}^{n}\,d\xi=\int_{n^{G}(\mathcal{U})}\exp\left(-n\left(\frac{\widetilde{P}(n^{-G}\xi)}{100}-\widetilde{R}(n^{-G}\xi)\right)\right)\,d\xi$

for each $n\in\mathbb{N}_{+}$ where $\mathcal{U}=T^{-1}(\mathcal{O})$. Upon
noting that $T(n^{-G}\xi)\in\mathcal{O}$ whenever $\xi\in n^{G}(\mathcal{U})$,
it follows from Lemma 3.3 and (41) that

$\exp\left(-n\left(\frac{\widetilde{P}(n^{-G}\xi)}{100}-\widetilde{R}(n^{-G}\xi)\right)\right)\chi_{{}_{n^{G}(\mathcal{U})}}(\xi)\leq
e^{-n\gamma\widetilde{P}(n^{-G}\xi)}\chi_{{}_{n^{G}(\mathcal{U})}}(\xi)\leq
e^{M}e^{-C(\eta^{2}+\zeta^{4})}$

for all $n\in\mathbb{N}_{+}$ and $\xi=(\eta,\zeta)\in\mathbb{R}^{2}$; here $M$
and $C$ are positive constants and, for a measurable set
$\mathcal{A}\subseteq\mathbb{R}^{2}$, $\chi_{{}_{\mathcal{A}}}$ denotes the
indicator function of $\mathcal{A}$. Also, thanks to the computation (40) and
the fact that $\widetilde{R}_{2}(\xi)=o(\widetilde{P}(\xi))$ as $\xi\to 0$, we
have

$\lim_{n\to\infty}\exp\left(-n\left(\frac{\widetilde{P}(n^{-G}\xi)}{100}-\widetilde{R}(n^{-G}\xi)\right)\right)\chi_{{}_{n^{G}(\mathcal{U})}}(\xi)=0,$

for almost every $\xi\in\mathbb{R}^{d}$. Given that $\xi\mapsto
e^{M}e^{-C(\eta^{2}+\zeta^{4})}\in L^{1}(\mathbb{R}^{2})$, an appeal to the
dominated convergence theorem is justified and we conclude that

$\displaystyle\lim_{n\to\infty}n^{5/8}\|\psi^{(n)}\|_{\infty}$
$\displaystyle=$
$\displaystyle\lim_{n\to\infty}n^{5/8}\rho^{n}+\lim_{n\to\infty}n^{-5/8}\int_{\mathcal{O}}\absolutevalue{\widehat{\psi}(\xi)}^{n}\,d\xi$
$\displaystyle=$ $\displaystyle
0+\lim_{n\to\infty}\int_{n^{G}(\mathcal{U})}\exp\left(-n\left(\frac{\widetilde{P}(n^{-G}\xi)}{100}-\widetilde{R}(n^{-G}\xi)\right)\right)\,d\xi$
$\displaystyle=$ $\displaystyle 0.$

∎

To our knowledge, there is no known theory that is able to treat this example
and, in particular, the asymptotic behavior of $\|\psi^{(n)}\|_{\infty}$
remains unknown. ∎

##  6 Discussion

###  6.1 Another perspective

As we noted in the introduction and demonstrated in the proof of Theorem 3.1,
key to the analysis in this paper is the observation that

$\varphi(t)=H_{P}^{t}(0)=H_{\widetilde{P}}^{t}(0)=\widetilde{\varphi}(t)$

for all $t>0$ where $\widetilde{P}=P\circ T$ can be seen, in some sense, as a
“dual” symbol to $P$. In terms of our motivating example in which
$P(\xi)=(\eta+\zeta^{2})^{2}+\eta^{4}$ and

$\Lambda=\partial_{x_{1}}^{4}+\partial_{x_{2}}^{4}+2i\partial_{x_{1}}\partial_{x_{2}}^{2}-\partial_{x_{1}}^{2},$

$\widetilde{P}(\xi)=\eta^{2}+(\eta-\zeta^{2})^{4}$ for
$\xi=(\eta,\zeta)\in\mathbb{R}^{2}$ and corresponds to the constant
coefficient operator

$\widetilde{\Lambda}=-\partial_{x_{1}}^{2}+\partial_{x_{1}}^{4}+4i\partial_{x_{1}}^{3}\partial_{x_{2}}^{2}-6\partial_{x_{1}}^{2}\partial_{x_{2}}^{4}-4i\partial_{x_{1}}\partial_{x_{2}}^{6}+\partial_{x_{2}}^{8}.$

Though the operator $\Lambda$ is a self-adjoint fourth-order elliptic operator
and is therefore well understood, $\widetilde{\Lambda}$ is mysterious and
somewhat poorly behaved as the following proposition shows.

######  Proposition 6.1.

The operator $\widetilde{\Lambda}$ with symbol
$\widetilde{P}=\eta^{2}+(\eta-\zeta^{2})^{4}$ is not hypoelliptic. Still,
$-\widetilde{\Lambda}$ generates a semigroup
$\\{e^{-t\widetilde{\Lambda}}\\}_{t>0}$ with heat kernel given by

$H_{\widetilde{P}}^{t}(x)=\frac{1}{(2\pi)^{2}}\int_{\mathbb{R}^{2}}e^{-t\widetilde{P}(\xi)}e^{-ix\cdot\xi}\,d\xi$

defined for $t>0$ and $x\in\mathbb{R}^{2}$. This heat kernel has

$\widetilde{\varphi}(t)=H_{\widetilde{P}}^{t}(0)\asymp\begin{cases}t^{-1/2}&0<t\leq
1\\\ t^{-5/8}&t\geq 1\end{cases}$

for $t>0$ and, further,

$\lim_{t\to
0}t^{1/2}\widetilde{\varphi}(t)=\frac{1}{\pi^{2}}\Gamma(5/4)^{2}\hskip
28.45274pt\mbox{and}\hskip
28.45274pt\lim_{t\to\infty}t^{5/8}\widetilde{\varphi}(t)=\frac{1}{2\pi^{3/2}}\Gamma(9/8).$

###### Proof.

The statements concerning $H_{\widetilde{P}}$ and $\widetilde{\varphi}$ follow
directly from our analysis in Examples 7 and 9 through the identity
$\varphi=\widetilde{\varphi}$. Thus, it remains to prove that
$\widetilde{\Lambda}$ is not hypoelliptic. To see this, we compute

$\partial^{4}_{\zeta}\widetilde{P}(\eta,\zeta)=384\zeta^{4}-1152\zeta^{2}(\eta-\zeta^{2})+144(\zeta-\zeta^{2})^{2}$

to see that

$\lim_{\zeta\to\infty}\frac{\partial^{4}_{\zeta}\widetilde{P}(\zeta^{2},\zeta)}{\widetilde{P}(\zeta^{2},\zeta)}=384.$

Thus, for the multi-index $\alpha=(0,4)$,

$\lim_{|\xi|\to\infty}\frac{\widetilde{P}^{(\alpha)}(\xi)}{\widetilde{P}(\xi)}\neq
0;$

here, for a multi-index $\alpha=(\alpha_{1},\alpha_{2})$,
$\widetilde{P}^{(\alpha)}=\partial_{\eta}^{\alpha_{1}}\partial_{\zeta}^{\alpha_{2}}\widetilde{P}$
in the notation of L. Hörmander [Hö83]. By virtue of Theorem 11.1.1 of [Hö83],
we conclude that $\widetilde{\Lambda}$ is not hypoelliptic. ∎

Of course, as we discussed in the introduction, the theory formulated in this
article has a version in which one begins with the operator
$\widetilde{\Lambda}$ whose heat kernel $H_{\widetilde{P}}$, along the
diagonal, is well behaved in large time but whose small time behavior is
elusive and is deduced more easily from $H_{P}$ via the correspondence
$H_{P}^{t}(0)=H_{\widetilde{P}}^{t}(0)$.

###  6.2 Future Directions

This article has focused on on-diagonal large-time asymptotics for the heat
kernels $H_{P}$ of certain inhomogeneous operators $\Lambda$ and symbols $P$
and well-behaved associated perturbations. At present, we do not have a good
understanding of the off-diagonal behavior of these kernels in large time.
Knowledge of this behavior could greatly improve our understanding of the
theory presented in this article and it would also inform our study of
convolution powers. In particular, having a good handle on this off-diagonal
behavior would inform on the stability of convolution powers (seen as
numerical difference schemes) as it has in [Th65], [CF22] (see also [Th69] and
Section 6 of [RSC17]).

The symbols $P(\xi)$ treated throughout this article all take the special form
$P_{1}(\eta+Q(\zeta))+P_{2}(\eta)$ resemblant of our motivating example in
which $P(\eta,\zeta)=(\eta+\zeta^{2})^{2}+\eta^{4}$. Our methods appear to be
more broadly applicable, however. Consider, for example, the symbol

$P^{\prime}(\xi)=(\eta+\zeta^{2})+\eta^{2}\zeta^{2}$

defined for $\xi=(\eta,\zeta)\in\mathbb{R}^{2}$. Akin to our motivating
example, $P^{\prime}$ lacks a tractable scaling from which large-time
asymptotics for $\varphi^{\prime}(t)=H_{P^{\prime}}^{t}(0)$ may be computed
directly. Still, by composing $P^{\prime}$ with the measure-preserving
transformation $T(\eta,\zeta)=(\eta-\zeta^{2},\zeta)$, we find that

$\lim_{t\to\infty}t\widetilde{P^{\prime}}(t^{-1/2}\eta,t^{-1/6}\zeta)=\eta^{2}+\zeta^{6}$

for $(\eta,\zeta)\in\mathbb{R}^{2}$ where
$\widetilde{P^{\prime}}=P^{\prime}\circ T$. Accompanied by the fact that
$t\widetilde{P^{\prime}}(t^{-1/2}\eta,t^{-1/6}\zeta)\geq(\eta^{2}+\zeta^{4}-2)/8$
for $t\geq 1$ and $(\eta,\zeta)\in\mathbb{R}^{2}$, the change of variables
$(\eta,\zeta)\mapsto(t^{-1/2}\eta,t^{-1/6}\zeta)$ followed by an application
of the dominated convergence theorem shows that

$\lim_{t\to\infty}t^{2/3}\varphi^{\prime}(t)=\frac{1}{(2\pi)^{2}}\int_{\mathbb{R}^{2}}e^{-\eta^{2}-\zeta^{6}}\,d\eta\,d\zeta=\frac{1}{2\pi^{3/2}}\Gamma(7/6).$

Another class of example can be produced easily by considering powers of
symbols $P(\xi)=P_{1}(\eta+Q(\zeta))+P_{2}(\eta)$ on $\mathbb{R}^{d}$
satisfying the hypotheses of Theorem 3.1. Specifically, given such a symbol
$P(\xi)$ and $\kappa>0$, our methods show that the heat kernel
$H_{P^{\kappa}}$ associated to $P(\xi)^{\kappa}$ satisfies the on-diagonal
asymptotics

$H_{P^{\kappa}}^{t}(0)\asymp\begin{cases}t^{-\mu_{0}/\kappa}&0<t\leq 1\\\
t^{-\mu_{\infty}/\kappa}&t\geq 1\end{cases}$

for $t>0$ where $\mu_{0}$ and $\mu_{\infty}$ are those given by Theorem 3.1.
At this time, describing precisely the set of examples to which our methods
apply is an open question. These questions will be explored in a forthcoming
article.

##  Appendix A Technical Estimates

######  Lemma A.1.

Let $P$ be a positive homogeneous function and let
$F\in\operatorname{End}(\mathbb{R}^{a})$ be such that $\\{t^{F}\\}$ is non-
expanding. If there exits $E\in\operatorname{Exp}(P)$ such that $[E,F]=0$,
then, for any $\delta>0$, there exists $\rho>0$ for which

$\rho P(t^{F}\eta)\leq\delta P(\eta)$

for all $\eta\in\mathbb{R}^{a}$ and $0<t\leq 1$.

###### Proof.

Because $P(0)=0$, it suffices to prove that, for some $C>0$,

$P(t^{F}\eta)\leq CP(\eta)$

for all $0<t\leq 1$ and non-zero $\eta\in\mathbb{R}^{a}$. Denote by $S$ the
unital level set of $P$ and observe that

$K=\\{t^{F}\eta:0<t\leq 1,\,\eta\in S\\}$

is a bounded set because $S$ is compact and $\\{t^{F}\\}_{0<t\leq 1}$ is
uniformly bounded. By virtue of the continuity of $P$, it follows that, for
every $0<t\leq 1$ and $\eta\in S$,

$P(t^{F}\eta)\leq\sup_{\xi\in K}P(\xi)=:C<\infty.$

Given any non-zero $\eta\in\mathbb{R}^{a}$, Proposition 4.1 of [BR22]
guarantees that $\eta=r^{E}\eta_{0}$ where $r=P(\eta)>0$, $\eta_{0}\in S$, and
$E$ is as in the statement of the lemma. Using the fact that $E$ and $F$
commute, we have

$P(t^{F}\eta)=P(r^{E}t^{F}\eta_{0})=rP(t^{F}\eta_{0})\leq rC=CP(\eta),$

as desired. ∎

######  Lemma A.2.

Let $P$ be a positive homogeneous function on $\mathbb{R}^{a}$. Then, for any
positive constants $\rho_{1},\rho_{2}$, there exists $\epsilon>0$ for which

$\epsilon P(\xi)\leq\rho_{1}P(\zeta+\xi)+\rho_{2}P(\zeta)$

for all $\xi,\zeta\in\mathbb{R}^{a}$.

###### Proof.

Let $\rho_{1},\rho_{2}>0$ and define $R:\mathbb{R}^{2a}\to\mathbb{R}$ by

$R(\xi,\zeta)=\rho_{1}P(\zeta+\xi)+\rho_{2}P(\zeta)$

for $(\xi,\zeta)\in\mathbb{R}^{a}\times\mathbb{R}^{a}=\mathbb{R}^{2a}$. It is
clear that $R$ is a continuous and positive-definite function on
$\mathbb{R}^{2a}$. For any $E\in\operatorname{Exp}(P)$, observe that

$R(t^{E\oplus
E}(\xi,\zeta))=R(t^{E}\xi,t^{E}\zeta)=\rho_{1}P(t^{E}\zeta+t^{E}\xi)+\rho_{2}P(t^{E}\zeta)=tR(\xi,\zeta)$

for all $t>0$ and $(\xi,\zeta)\in\mathbb{R}^{2a}$. Thus, $R$ is homogeneous
with respect to $E\oplus E$ and, because $\\{t^{E}\\}$ is a contracting group
on $\mathbb{R}^{a}$, it is evident that $\\{t^{E\oplus E}\\}$ is a contracting
group on $\mathbb{R}^{2a}$. In view of Proposition 1.1 of [BR22], we conclude
that $R$ is a positive homogeneous function. Denote by $S$ the unital level
set of $R$ and, because $S$ is compact and $(\xi,\zeta)\mapsto P(\xi)$ is
continuous and does not identically vanish on $S$, we have

$\epsilon:=\left(\sup_{(\xi,\zeta)\in S}P(\xi)\right)^{-1}>0.$

Given a non-zero $(\xi,\zeta)\in\mathbb{R}^{2a}$, Proposition 4.1 of [BR22]
ensures that $(\xi,\zeta)=r^{E\oplus
E}(\xi_{0},\zeta_{0})=(r^{E}\xi_{0},r^{E}\zeta_{0})$ where $r=R(\xi,\zeta)>0$,
$(\xi_{0},\zeta_{0})\in S$ and $E\in\operatorname{Exp}(P)$. With this, we
observe that

$\epsilon P(\xi)=\epsilon rP(\xi_{0})\leq
r=R(\xi,\zeta)=\rho_{1}P(\zeta+\xi)+\rho_{2}P(\zeta).$

Since this inequality holds trivially when $(\xi,\zeta)=(0,0)$, the proof is
complete. ∎

By making analogous arguments to those in the proof above, we easily obtain
the following lemma.

######  Lemma A.3.

Let $P$ be a positive homogeneous function on $\mathbb{R}^{a}$. Then, for any
positive constants $\rho_{1},\rho_{2}$, there exists $\epsilon>0$ for which

$\epsilon P(\zeta+\xi)\leq\rho_{1}P(\xi)+\rho_{2}P(\zeta)$

for all $\xi,\zeta\in\mathbb{R}^{a}$.

######  Lemma A.4.

Let $P_{1}$ and $P_{2}$ be positive homogeneous functions on $\mathbb{R}^{a}$
and let $E_{1}\in\operatorname{Exp}(P_{1})$ and
$E_{2}\in\operatorname{Exp}(P_{2})$. If $[E_{1},E_{2}]=0$ and
$\\{t^{E_{1}-E_{2}}\\}$ is non-expanding, then, for any $\delta>0$, there are
constants $\rho,M>0$ for which

$\rho P_{1}(t^{E_{1}-E_{2}}\eta)\leq M+\delta P_{2}(\eta)$

for all $\eta\in\mathbb{R}^{a}$ and $0<t\leq 1$.

###### Proof.

Let $\delta>0$. We shall first prove that, there are constants $\rho,M>0$ for
which

$\rho P_{1}(\eta)\leq M+\delta P_{2}(\eta)$

for all $\eta\in\mathbb{R}^{a}$. We note that this is the desired inequality
evaluated at $t=1$; we shall obtain the full inequality by scaling. Now,
denote by $S$ and $B$ the unital level set and unit ball associated to
$P=\delta P_{2}$. Because $B$ is relatively compact and
$\\{t^{E_{1}-E_{2}}\\}_{0<t\leq 1}$ is uniformly bounded,

$K=\\{r^{E_{1}-E_{2}}\eta:0<r\leq 1,\,\eta\in B\\}$

is bounded and hence relatively compact. Since $P_{2}$ is positive and does
not vanish identically on $K$,

$\rho:=\left(\sup_{\xi\in K}P_{1}(\xi)\right)^{-1}>0.$

Now, for any $\eta\in\mathbb{R}^{a}\setminus B$, it follows from Proposition
4.1 of [BR22] that $\eta=r^{-E_{2}}\eta_{0}$ for $r=1/P(\eta)=1/(\delta
P_{2}(\eta))\leq 1$ and $\eta_{0}\in S$. Consequently,

$\rho P_{1}(\eta)=\rho
P_{1}(r^{-E_{2}}\eta_{0})=\frac{\rho}{r}P_{1}(r^{E_{1}-E_{2}}\eta_{0})\leq\frac{1}{r}=\delta
P_{2}(\eta).$

Of course, $P_{1}(\eta)$ is bounded on $B$ and so it follows that, for some
constant $M>0$,

$\rho P_{1}(\eta)\leq M+\delta P_{2}(\eta)$

for all $\eta\in\mathbb{R}^{a}$, as claimed. Finally, for $0<t\leq 1$ and
$\eta\in\mathbb{R}^{a}$, we apply the preceding inequality to $t^{-E_{2}}\eta$
to see that

$\rho P_{1}(t^{E_{1}-E_{2}}\eta)=t\rho P_{1}(t^{-E_{2}}\eta)\leq
t\left(M+\delta P_{2}(t^{-E_{2}}\eta)\right)\leq tM+\delta P_{2}(\eta)\leq
M+\delta P_{2}(\eta)$

as desired. ∎

###### Proof of Lemma 3.3.

We shall first prove the estimate (12). For the lower estimate in (12),
simultaneous appeals to Lemmas A.1 and A.4 give positive constants $\rho_{1}$,
$\rho_{2}$, and $M$ for which

$\rho_{1}P_{1}(\zeta)\leq M+P_{2}(\zeta)\hskip 28.45274pt\mbox{and}\hskip
28.45274pt\rho_{2}P_{1}(t^{E_{1}-E_{2}}\eta)\leq P_{1}(\eta)/2$

for $\zeta,\eta\in\mathbb{R}^{a}$ and $0<t\leq 1$. With the constants
$\rho_{1}$ and $\rho_{2}$ in hand, we appeal to Lemma A.2 to find a positive
constant $\epsilon>0$ for which

$\epsilon P_{1}(\xi)\leq\rho_{1}P_{1}(\zeta+\xi)+\rho_{2}P_{1}(\zeta)$

for all $\xi,\zeta\in\mathbb{R}^{a}$. Consequently, for
$\xi,\eta\in\mathbb{R}^{a}$ and $0<t\leq 1$,

$\epsilon
P_{1}(\xi)\leq\rho_{1}P_{1}(t^{E_{1}-E_{2}}\eta+\xi)+\rho_{2}P_{1}(t^{E_{1}-E_{2}}\eta)\leq
M+P_{2}(t^{E_{1}-E_{2}}\eta+\xi)+P_{1}(\eta)/2.$

Upon setting $C=\min\\{\epsilon,1/2\\}>0$, it follows that

$C(P_{1}(\eta)+P_{1}(\xi))-M\leq P_{1}(\eta)/2+\epsilon P_{1}(\xi)-M\leq
P_{1}(\eta)+P_{2}(t^{E_{1}-E_{2}}\eta+\xi)$

for all $\eta,\xi\in\mathbb{R}^{a}$ and $0<t\leq 1$. For the upper estimate in
(12), we first make an appeal to Lemma A.4 to obtain positive constants
$M^{\prime}$ and $C_{1}$ for which

$P_{1}(\eta)\leq M^{\prime}+C_{1}P_{2}(\eta)$

for all $\eta\in\mathbb{R}^{a}$. By virtue of Lemma A.3, there is a positive
constant $C_{2}$ for which

$P_{2}(\zeta+\xi)\leq C_{2}(P_{2}(\zeta)+P_{2}(\xi))$

for all $\xi,\zeta\in\mathbb{R}^{a}$. Finally, an appeal to Lemma A.1
guarantees a positive constant $C_{3}$ for which

$P_{2}(t^{E_{1}-E_{2}}\eta)\leq C_{3}P_{2}(\eta)$

for all $\eta\in\mathbb{R}^{a}$ and $0<t\leq 1$. Upon combining these
estimates, we obtain

$\displaystyle P_{1}(\eta)+P_{2}(t^{E_{1}-E_{2}}\eta+\xi)$ $\displaystyle\leq$
$\displaystyle
M^{\prime}+C_{1}P_{2}(\eta)+C_{2}P_{2}(t^{E_{1}-E_{2}}\eta)+C_{2}P_{2}(\xi)$
$\displaystyle\leq$ $\displaystyle
M^{\prime}+(C_{1}+C_{2}C_{3})P_{2}(\eta)+C_{2}P_{2}(\xi)$

for all $\eta,\xi\in\mathbb{R}^{a}$ and $0<t\leq 1$. Our desired (upper)
estimate in (12) follows immediately by taking
$C^{\prime}=\max\\{C_{2},C_{1}+C_{2}C_{3}\\}$ and so the proof of (12) is
complete.

To prove (11), we first make an appeal to Lemma A.4 to obtain positive
constants $\rho$ and $M$ for which

$\rho P_{1}(t^{E_{1}-E_{2}}\eta)\leq M+P_{2}(\eta)/2$

for all $\eta\in\mathbb{R}^{a}$ and $0<t\leq 1$. With this $\rho$, an appeal
to Lemma A.2 yields $\epsilon>0$ for which

$\epsilon P_{1}(\xi)\leq P_{1}(\zeta+\xi)+\rho P_{1}(\zeta)$

for all $\xi,\zeta\in\mathbb{R}^{a}$. Upon setting $C=\min\\{\epsilon,1/2\\}$,
we obtain

$\displaystyle CP_{1}(\xi)+P_{2}(\eta))-M$ $\displaystyle\leq$
$\displaystyle\epsilon P_{1}(\xi)+P_{2}(\eta)/2-M$ $\displaystyle\leq$
$\displaystyle P_{1}(t^{E_{1}-E_{2}}\eta+\xi)+\rho
P_{1}(t^{E_{1}-E_{2}}\eta)+P_{2}(\eta)/2-M$ $\displaystyle\leq$ $\displaystyle
P_{1}(t^{E_{1}-E_{2}}\eta+\xi)+P_{2}(\eta)$

for all $\eta,\xi\in\mathbb{R}^{a}$ and $0<t\leq 1$ and this is precisely the
lower estimate in (12). Making use of Lemmas A.4 and A.3, the upper estimate
in (12) is established in a similar way to that for (11); we leave the
remaining details to the reader. ∎

Our final goal in this appendix is to prove Proposition 4.3. Before the proof,
we present two lemmas.

######  Lemma A.5.

Assume the notation and hypotheses of Proposition 4.3 and define

$\widetilde{P}_{\infty}(\xi)=P_{1}(\eta)+P_{2}(-Q(\zeta))\hskip
28.45274pt\mbox{and}\hskip
28.45274pt\widetilde{\mathcal{E}}_{\infty}(\xi)=\widetilde{P}(\xi)-\widetilde{P}_{\infty}(\xi)$

for $\xi=(\eta,\zeta)\in\mathbb{R}^{d}$. Then
$\widetilde{\mathcal{E}}_{\infty}(\xi)$ is subhomogeneous with respect to
$G=E_{1}\oplus F_{2}$.

###### Proof.

Let $\epsilon>0$ and $K\subseteq\mathbb{R}^{d}$ be a compact set. Since
$\\{t^{E_{1}-E_{2}}\\}$ is contracting and $P_{2}$ and $Q$ are continuous, we
can find $0<t_{0}\leq 1$ for which

$\absolutevalue{(P_{1}(\eta)+P_{2}(t^{E_{1}-E_{2}}\eta-Q(\zeta))-\widetilde{P}_{\infty}(\xi)}=\absolutevalue{P_{2}(t^{E_{1}-E_{2}}\eta-Q(\zeta))-P_{2}(-Q(\zeta))}<\epsilon$

for all $0<t\leq t_{0}$ and $\xi=(\eta,\zeta)\in K$. Thus, for every $0<t\leq
t_{0}$ and $\xi=(\eta,\zeta)\in K$, we have

$\displaystyle\absolutevalue{\widetilde{\mathcal{E}}_{\infty}(t^{G}\xi)}$
$\displaystyle=$
$\displaystyle\absolutevalue{P_{1}(t^{E_{1}}\eta)+P_{2}(t^{E_{1}}\eta-Q(t^{F_{2}}\zeta))-P_{1}(t^{E_{1}}\eta)-P_{2}(-Q(t^{F_{2}}\zeta))}$
$\displaystyle=$
$\displaystyle\absolutevalue{P_{2}(t^{E_{1}}\eta-t^{E_{2}}Q(\zeta))-P_{2}(-t^{E_{2}}Q(\zeta))}$
$\displaystyle=$ $\displaystyle
t\absolutevalue{P_{2}(t^{E_{1}-E_{2}}\eta-Q(\zeta)-P_{2}(-Q(\zeta))}$
$\displaystyle\leq$ $\displaystyle\epsilon t,$

as desired. ∎

The following lemma asserts loosely that $\widetilde{P}$ is approximately
homogeneous with respect to $G$ in small time whenever $\\{t^{E_{1}-E_{2}}\\}$
is contracting.

######  Lemma A.6.

Assume the notation and hypotheses of Proposition 4.3. Then, there exists a
compact set $K\subseteq\mathbb{R}^{d}$ such that, for any $\tau>0$ and
$0<\epsilon<1$, there exists $0<t_{0}\leq\tau$ for which

$\mathcal{O}:=\\{\xi\in\mathbb{R}^{d}:\xi=0\,\,\mbox{or}\,\,\xi=t^{G}\xi^{\prime}\,\,\mbox{for}\,\,0<t<t_{0}\,\,\mbox{and}\,\,\xi^{\prime}\in
K\\}$

is an open neighborhood of $0$ and, for each non-zero
$\xi=t^{G}\xi^{\prime}\in\mathcal{O}$,

$(1-\epsilon)t\leq\widetilde{P}(\xi)=\widetilde{P}(t^{G}\xi^{\prime})\leq(1+\epsilon)t$

###### Proof.

Let $K=S_{\widetilde{P}_{\infty}}$ be the unital level set of
$\widetilde{P}_{\infty}$ and, by an appeal to the preceding lemma, let
$0<t_{0}\leq\tau$ be such that

$\absolutevalue{\widetilde{\mathcal{E}}_{\infty}(t^{G}\xi^{\prime})}\leq\epsilon
t$

whenever $0<t\leq t_{0}$ and $\xi^{\prime}\in K$. In this notation, we observe
that $\mathcal{O}$, as defined in the statement, coincides with the
$\widetilde{P}_{\infty}$-adapted open ball
$B_{t_{0}}=\\{\xi\in\mathbb{R}^{d}:\widetilde{P}_{\infty}(\xi)<t_{0}\\}$ since
$G\in\operatorname{Exp}(\widetilde{P}_{\infty})$; in particular,
$\mathcal{O}=B_{t_{0}}$ is necessarily an neighborhood of $0$ (see Proposition
4.1 of [BR22]). Further, for each non-zero
$\xi=t^{G}\xi^{\prime}\in\mathcal{O}=B_{t_{0}}$, we have
$\widetilde{P}_{\infty}(\xi)=t$ and therefore

$\absolutevalue{\widetilde{P}(\xi)-t}=\absolutevalue{\widetilde{\mathcal{E}}_{\infty}(\xi)}=\absolutevalue{\widetilde{\mathcal{E}}_{\infty}(t^{G}\xi^{\prime})}\leq\epsilon
t$

so that
$(1-\epsilon)t\leq\widetilde{P}(\xi)=\widetilde{P}(t^{G}\xi^{\prime})\leq(1+\epsilon)t$.
∎

###### Proof of Proposition 4.3.

In view of Lemma 4.2, we need only to prove the sufficiency of the
subhomogeneity condition. Specifically, we prove that
$\widetilde{R}(\xi)=o(\widetilde{P}(\xi))$ as $\xi\to 0$ when
$\\{t^{E_{1}-E_{2}}\\}$ is contracting and $\widetilde{R}$ is subhomogeneous
with respect to $G$. To this end, let $\epsilon>0$ and take $K$ as in the
preceding lemma. With our assumption that $\widetilde{R}$ is subhomogeneous
with respect to $G$, let $\tau>0$ be such that
$\absolutevalue{\widetilde{R}(t^{G}\xi^{\prime})}\leq\epsilon t/2$ whenever
$\xi^{\prime}\in K$ and $0<t\leq\tau$. By an appeal to the preceding lemma,
let $0<t_{0}\leq\tau$ for which

$\widetilde{P}(\xi)=\widetilde{P}(t^{G}\xi^{\prime})\geq t/2$

whenever $\xi=t^{G}\xi^{\prime}$ is a non-zero member of the open neighborhood

$\mathcal{O}=\\{\xi\in\mathbb{R}^{d}:\xi=0\,\,\mbox{or}\,\,\xi=t^{G}\xi^{\prime}\,\,\mbox{for}\,\,0<t<t_{0}\,\,\mbox{and}\,\,\xi^{\prime}\in
K\\}$

of $0$ in $\mathbb{R}^{d}$. Thus, for any non-zero
$\xi=t^{G}\xi^{\prime}\in\mathcal{O}$,

$\absolutevalue{\widetilde{R}(\xi)}=\absolutevalue{\widetilde{R}(t^{G}\xi^{\prime})}\leq\epsilon
t/2\leq\epsilon\widetilde{P}(t^{G}\xi^{\prime})=\epsilon\widetilde{P}(\xi).$

By the continuity of $\widetilde{R}$ and $\widetilde{P}$, this estimate
clearly holds when $\xi=0$ and so we have shown that
$\widetilde{R}(\xi)=o(\widetilde{P}(\xi))$ as $\xi\to 0$ as was asserted. ∎

Acknowlegements: Evan Randles would like to thank Cornell University and Colby
College for financial support and Cornell University for hosting him during
his 2021-2022 sabbatical leave from Colby College. Laurent Saloff-Coste is
partially supported by NSF grant DMS-1707589 and DMS-2054593.

## References

  * [Ba01] G. Barbatis  Explicit estimates on the fundamental solution of higher-order parabolic equations with measurable coefficients.  Journal of Differential Equations, 174(2):442-463, 2001.
  * [Ba98] G. Barbatis.  Sharp heat kernel bounds and Finsler-type metrics.  The Quarterly Journal of Mathematics, 49(3):261-277, 1998.
  * [BD96] G. Barbatis and E. B. Davies  Sharp bounds on heat kernels of higher order uniformly elliptic operators.  Journal of Operator Theory, 36(1):197-198, 1996.
  * [BR22] H. Q. Bui and E. Randles.  A Generalized Polar-Coordinate Integration Formula with Applications to the Study of Convolution Powers of Complex-Valued Functions on $\mathbb{Z}^{d}$.  Journal of Fourier Analysis and Applications, 28(2):1-74, 2022.
  * [CKS87] E. A. Carlen, S. Kusuoka, and D. Stroock.  Upper bounds for symmetric Markov transition functions.  In Annales de l’IHP Probabilités et statistiques, 23(S2):245-287.
  * [Coe22] L. Coeuret.  Local limit theorem for complex valued sequences.  arXiv:2201.01514
  * [Co96] T. Coulhon.  Ultracontractivity and Nash type inequalities.  Journal of functional analysis, 141(2):510-539, 1996.
  * [Da89] E. B. Davies.  Heat kernels and spectral theory. Cambridge university press; 1989.
  * [Da95] E. B. Davies.  Uniformly Elliptic Operators With Measurable Coefficients.  J. Funct. Anal., 132:141-169, 1995.
  * [Da95a] E. B. Davies.  Long-time Asymptotics of Fourth Order Parabolic Equations.  J. Math. Anal., 67:323–345, 1995.
  * [DSC14] P. Diaconis and L. Saloff-Coste.  Convolution powers of complex functions on Z.  Math. Nachrichten, 287(10):1106–1130, 2014.
  * [Ei69] S. D. Eidelman.  Parabolic Systems North-Holland/Wolters-Noordhoff, Amsterdam, 1969.
  * [EN00] K.-J. Engel and R. Nagel.  One-Parameter Semigroups for Linear Evolution Equations, volume 194 of Graduate Texts in Mathematics.  Springer New York, 2000.
  * [EP70] M. A. Evgrafov and M. M. Postnikov.  Asymptotic behaviour of Green’s functions for parabolic and elliptic equations with constant coefficients’  Math. USSR Sbomik 11:1-24, 1970.
  * [CF22] J.-F. Coulombel and G. Faye.  Generalized Gaussian bounds for discrete convolutions powers.  Rev. Mat. Iberoam. doi: 10.4171/RMI/1338 (to appear in print, 2022)
  * [Hö83] L. Hörmander.  The Analysis of Linear Partial Differential Operators II.  Springer-Verlag Berlin Heidelberg, Berlin, 1983.
  * [Ra22] E. Randles.  Local Limit Theorems for Complex Functions on $\mathbb{Z}^{d}$.  arXiv:2201.01319 [math.CA], 2022.
  * [RSC15] E. Randles and L. Saloff-Coste.  On the convolution powers of complex functions on $\mathbb{Z}$.  J. Fourier Anal. Appl., 21(4):754–798, 2015.
  * [RSC17] E. Randles and L. Saloff-Coste.  Convolution powers of complex functions on $\mathbb{Z}^{d}$.  Rev. Matemática Iberoam., 33(3):1045–1121, 2017.
  * [RSC17a] E. Randles and L. Saloff-Coste.  Positive-homogeneous operators, heat kernel estimates and the Legendre-Fenchel transform.  In Baudoin F., Peterson, J. (eds) Stochastic Analysis and Related Topics: A Festschrift in Honor of Rodrigo Bañuelos., volume 72 of Progress in Probability, pages 1–47. Springer International Publishing, 2017.
  * [RSC20] E. Randles and L. Saloff-Coste.  Davies’ method for heat-kernel estimates: An extension to the semi-elliptic setting.  Transactions of the American Mathematical Society, 373(4):2525–2565, 2020.
  * [SC09] L. Saloff-Coste.  Sobolev Inequalities in Familiar and Unfamiliar Settings.  In: Sobolev Spaces In Mathematics I. International Mathematical Series, vol 8. Springer, New York, NY, 2009.
  * [Sp64] F. Spitzer.  Principles of Random Walk, volume 34 of Graduate Texts in Mathematics.  Springer New York, New York, NY, 1964.
  * [Th65] V. Thomée.  Stability of difference schemes in the maximum-norm.  J. Differ. Equ., 1(3):273–292, July 1965.
  * [Th69] V. Thomée.  Stability theory for partial difference operators.  SIAM Rev., 11(2):152–195, 1969.
  * [Ti82] K. Tintarev  Short time asymptotics for fundamental solutions of higher order parabolic equations,  Comm. Partial Differential Equations, 7:371-391, 1982.
  * [Va85] N. Varopoulos  Hardy-Littlewood theory for semigroups.  Journal of functional analysis, 63(2):240-260, 1985.