# Rapid Integrator for a Class of Multi-Contact Systems

Marion Anderson Department of Electrical and Computer Engineering, University
of Michigan, Ann Arbor, MI, USA<EMAIL_ADDRESS>Shai Revzen Department of
Electrical and Computer Engineering, University of Michigan, Ann Arbor, MI,
USA<EMAIL_ADDRESS>
###### Abstract

Many problems in robotics involve creating or breaking multiple contacts
nearly simultaneously or in an indeterminate order. We present a novel general
purpose numerical integrator based on the theory of Event Selected Systems
(ESS). Many multicontact models are ESS, which has recently been shown to
imply that despite a discontinuous vector field, the flow of these systems is
continuous, piecewise smooth, and has a well defined orbital derivative for
all trajectories, which can be rapidly computed. We provide an elementary
proof that our integrator is first-order accurate and verify numerically that
it is in fact second-order accurate as its construction anticipated. We also
compare our integrator, implemented in NumPy, to a MuJoCo simulation on models
with 2 to 100 contacts, and confirm that the increase in simulation time per
contact is nearly identical. The results suggest that this novel integrator
can be invaluable for modelling and control in many robotics applications.

##  1 Introduction

A defining property of robots is the ability to physically interact with their
environment. To control these interactions, roboticists usually rely on models
in the form of equations of motion, which allow a computer to anticipate the
trajectory robots and objects will follow. When this trajectory’s first-order
dependence on parameters and initial conditions (its “orbital derivative”) is
computable, the models become much more amenable to planning and control
optimization [40, 39, 18, 46]

These equations of motion are typically represented using ordinary
differential equations (ODEs). The ODEs which describe robot tasks where the
robot makes and breaks contact with its environment, such as legged locomotion
or dexterous gripping, commonly use “hard contacts” to describe the process of
making or breaking contact [40]. These models produce ODEs with a jump
discontinuity at the moment of contact (a “hybrid transition” at the crossing
of a “guard” surface) and are often viewed as a type of “hybrid systems” or
“non-smooth systems”, which require special attention to integrate accurately
(see [18] for terminology). One approach for dealing with this complication
has been to replace hard contact models with so-called “soft contact”
approximations [33, 45, 43], whose smoothness allows first-order dependencies
to be numerically computed. Unfortunately, whether hard or soft, it is common
for several contacts to be made or broken nearly simultaneously. Examples
include the tripod walking gaits of RHex [1], iSprawl [29], and Dante II [4],
humanoid grasping such as the Utah/MIT Dextrous Hand [26] and the Karlsruhe
dextrous hand [48], and the multitude of three-fingered dextrous manipulators
[5, 27, 38].

Although there have been numerous results recovering first-order
approximations for trajectories undergoing a single isolated hybrid transition
[17, 25], or even a pair of near simultaneous hybrid transitions [25], first
order approximations for trajectories undergoing three or more near-
simultaneous transitions have only recently been published [8, 13]. The
complexity of the problem with three or more transitions stems from the
combinatorial growth in the ordering of transition events [19]. For even the
simplest cases, it can be shown that the entire factorial number of transition
orderings is reachable by infinitesimally perturbing an initial condition.
Recent work has shown that common near-simultaneous collisions which arise in
robotics, under mild technical conditions, admit a first-order approximation
[8] called the Bouligand Derivative (herein referred to as the “B-Deriv”). For
these systems, the B-Deriv is continuous, piecewise linear, and
computationally equivalent to the one-sided difference quotient limit
(“directional derivative”) [8]. [13] developed an efficient algorithm for
computing the B-Deriv from this result. This algorithm does not require the
root finding numerical machinery commonly found in contact solvers [13, 32,
12, 44, 45].

Figure 1: A three-legged hopping robot we developed, which motivated much of
this work [36, 2]. It “preloads” its elastic legs prior to ground contact
using a cable assembly. The cable assembly can only exert force upward on the
legs. The legs immediately exert a force on the ground at the moment of
contact. As a result, the hopper contact force model is affine, and cannot be
readily solved using the machinery for linear complementarity problems [6].

Root finding is known to be a significantly expensive computation in solving
the mechanical forces that arise in multi-contact, and has drawn much
attention from the robotics community [6, 5]. Despite its slow speed, root
finding is used to achieve sufficient accuracy for simulation and control. The
speed-accuracy tradeoff is a common concept in numerical methods, and is often
used when evaluating numerical methods against one another [15]. Removing root
finding machinery without sacrificing (too much) accuracy could greatly reduce
some of the need for computing power in robotics, enabling faster simulation
and real-time control with slower (cheaper) on-board computers.

###  1.1 Related Work

Regardless of the contact model used, hybrid systems describing contact events
have historically been treated by roboticists in three ways: with impact maps
(sometimes called “reset maps”) [35, 47], as linear complimentarity problems
(LCPs) [3, 16], and as multi-valued maps, often called “differential
inclusions” [37, 41]. As each of these formalisms seeks to model similar
physical phenomena, there has been some work to find equivalence relations
between them [7].

LCPs are a broad class of matrix optimization problems, with much work done by
Cottle starting in the late 1960s [10, 9, 11]. As such, LCPs benefit from a
multitude of efficient solvers [11, Chps. 4,5]. They have much attention in
robotics due to their natural extension to resolving the forces that occur in
multi-contact [3, 16, 34]. Although LCPs have such a natural extension to
sticking contact, they can only approximate the forces for sliding contacts
with Coloumb friction [3], and cannot represent nonlinear contact dynamics
[6]. Regardless, they have achieved prolific use in both research and
industrial settings [32, 12, 44]. Differential inclusions, first introduced to
engineering by Moreau [37] can treat friction directly, and are considered to
most accurately describe the hybrid phenomena arising in robotics [30, 41, 6].
However, differential inclusions have not received much attention for robot
simulation. Modern rigid body software simulation suites either use soft
contact models (MuJoCo [45]), LCP solvers (DartSim [32], PyBullet [12]), or a
combination of both (Open Dynamics Engine [44]). Impact maps apply
discontinuous (impulsive) updates to the vector field, and have primarily been
used to model the single-contact phenomena that occur in bipedal locomotion
[41, 47, 20]. Recent work [8, 13] suggests that modelling multi-contact using
impacts maps may be more tractable than previously thought. Impact maps are
the most similar formalism to the work presented in this manuscript.

###  1.2 Contributions

We present, to our knowledge, the first integration algorithm to combine
conventional ODE integration within a single “hybrid domain” with rapid
resolution of multi-contact events that (at the limit) reproduces the correct
orbital derivative. We describe its practical use in Sec. 2 and its
implementation alongside an example code in Sec. 3. We numerically assess its
accuracy in Sec. 4, compare its performance to the robotic simulation software
MuJoCo in Sec. 5, and provide a simple proof of correctness to first-order in
Sec. 6.

##  2 Integrator Usage

This section is intended as a “manual” for the practitioner who wishes to use
our integrator with their hybrid system and to provide a precise definition of
requirements. See Sec. 5 for modelling examples.

###  2.1 Mathematical Preliminaries

Let $\mathfrak{X}\subseteq\mathbb{R}^{n}$ be a compact domain. Let
$\mathcal{H}\subseteq\mathbb{R}^{m}$ be an open neighborhood of
$0\in\mathbb{R}^{m}$.

(hybrid) dynamics vector field $\displaystyle f$
$\displaystyle:\mathfrak{X}\times\mathcal{H}\to\mathbb{R}^{n},$ event
function(s) $\displaystyle h$ $\displaystyle:\mathfrak{X}\to\mathcal{H},$

As $f$ and $h$ fully define the dynamics of interest, we will refer to the
pair $(f,h)$ as the “system.” We define the “guard” $\mathfrak{g}_{k}$ to be
the zero-level set of $h_{k}(\cdot)$, the $k$-th component of $h$, i.e.
$\mathfrak{g}_{k}:=\\{x\in\mathfrak{X}\;|\;h_{k}(x)=0\\}$. We will consider
various (compact) “$\varepsilon$-neighborhoods” of $\mathfrak{g}_{k}$ defined
using $h_{k}(\cdot)$ by
$N^{\varepsilon}:=\\{x\in\mathfrak{X}\;|\;-\varepsilon\leq
h(x)\leq\varepsilon\\}$ for $\varepsilon\geq 0$.

We further require that:

  1. 1.

$f(x,h)$ is $\mathcal{C}^{p}$ smooth in its first parameter.

  2. 2.

$f(x,y)=f(x,y^{\prime})$ if $y_{k}y^{\prime}_{k}>0$ for all $k$. In other
words, the $y\in\mathcal{H}$ dependence of $f(\cdot,y)$ is only on the signs
of the components of $y$.

  3. 3.

$h(\cdot)$ is $\mathcal{C}^{p}$.

  4. 4.

There exists some $c>0$ such that for all $y\in\mathcal{H}$, $k$, and for all
$x\in N^{\nu_{k}}$ we have $\nabla h_{k}(x)\cdot f(x,y)>c$.111This is known as
the “liveness condition” of ESS [8, Defn. 1].

  5. 5.

For every $k$ we require a $\nu_{k}>0$ such that for all (fixed)
$b\in\mathcal{H}$, and $x_{0}\in N^{\nu_{k}}$, the unique solution for the ODE
$\dot{x}=f(x,b)$ that passes through $x_{0}$ exists for the interval
$[-h_{k}(x_{0})/c,h_{k}(x_{0})/c]$. This condition ensures that all
trajectories starting in $N^{\nu_{k}}$ cross $\mathfrak{g}_{k}$.

By abuse of notation we define $f(x):=f(x,h(x))$, and note that $f(x)$ is
smooth except perhaps at the guards $\mathfrak{g}_{k}$. Since those are proper
co-dimension 1 manifolds (liveness ensures $\nabla h_{k}\neq 0$), the
dynamical system $\dot{x}=f(x)$ is a hybrid dynamical system in the sense of
[46] with guards $\mathfrak{g}_{k}$.

###  2.2 Practical Use

Our integrator requires definitions of $f$, $h$, the Jacobian of $h$,
$Dh:\mathfrak{X}\to L(\mathbb{R}^{m},\mathbb{R}^{n})$, and a precision
parameter, $\varepsilon>0$. Note that the Jacobian need not be closed-form.

After choosing appropriate event functions, $\varepsilon$ provides the primary
means to trade off speed and accuracy. Increasing $\varepsilon$ reduces the
number of integration steps performed at the risk of greater integration error
and vice-versa. For any state $x_{p}\in N^{\nu_{\varepsilon}}$, the user
should work to ensure $\|Dh(x_{p})\cdot f(x_{p})\|_{\infty}\leq\varepsilon$.
This is always achievable because given some $\bar{h}$, there exists
$\mu_{k}:=\sup_{\mathfrak{X}}\nabla{\bar{h}}_{k}(x)\cdot f(x)<\infty$ because
$\mathfrak{X}$ is compact. Define $h_{k}:=\varepsilon/\mu_{k}{\bar{h}}_{k}$ to
obtain the result.

##  3 Algorithm

###  3.1 Implementation

To our knowledge, the integration algorithm presented in this section is the
first to combine conventional ODE integration within a single hybrid domain
with rapid resolution of multi-contact events using the B-Deriv. In the case
of ESS, the B-Deriv is the first-order approximation to the flow, even in
those locations where it is not classically (Fréchet) differentiable [8]. We
use the B-Deriv to resolve multi-contact events using tangent line
approximations, effectively “projecting” the current state through the guards.
Our algorithm performs conventional integration and B-Deriv integration in an
alternating fashion. As such, we will refer to conventional integration as the
“integration step” and B-Deriv integration as the “projection step.” Due to
this alternating behavior, our algorithm utilizes a loop structure in
software, shown in Listing 1. At the start of each loop, the trajectory has
been integrated to time $t_{q}$ and state $x_{q}$.

Figure 2: Our algorithm integrating a piecewise-constant ESS (purple dotted
line). The system has two guards, $\mathfrak{g}_{1}$ (thick red line) and
$\mathfrak{g}_{2}$ (dashed thick blue line) with non-zero slope and intersect
at a point away from the origin. The integrator projects (purple line) when
the trajectory enters $N^{\nu_{k}}$, otherwise it performs conventional
integration (purple dots).

By requirement (4) when hybrid transitions occur they are (locally)
irrevocable. This implies that we need only consider each guard once, at most.
Our algorithm therefore maintains a pool of candidate guards which have not
been crossed. In the context of the projection step, we will refer to
candidate guards simply as guards for brevity. If $x_{q}$ lies in
$N^{\nu_{k}}$ of some guards, we determine (to first order) which guard will
be crossed first, and project to that guard with the first order approximation
of the impact map as in [47], removing it from the candidate pool. When
$x_{q}$ lies in $N^{\nu_{k}}$ for multiple guards, the loop structure of our
algorithm naturally performs the projection step multiple times consecutively,
equivalent to projecting across the guards in one computation. The projection
operation is the exact impact map for a constant vector field dynamic crossing
a co-dimension 1 hyper-plane guard. [13] showed that a suitably chosen
piecewise-constant vector field (the “corner vector field”) has the same
orbital B-Deriv as the original ESS did, and our algorithm employs a similar
intuition. If $x_{q}$ does not lie in any $N^{\nu_{k}}$ we integrate forward
with a conventional ODE integrator until the trajectory enters any
$N^{\nu_{k}}$ (we used dopri5 [21, II.4, p.178]; most other integration
schemes could also be used). To detect when the trajectory enters some
$N^{\nu_{k}}$, the conventional integrator must provide a means to detect such
an event and terminate integration. We used a termination condition detecting
if the trajectory had entered $N^{\nu_{k}}$ of at least one guard.

In common rhythmic hybrid systems, such as a 1D hopping robot, the same guard
(the ground) can be crossed multiple times in alternating directions (stance-
flight and flight-stance transitions) over a single trajectory. To integrate
such a system through multiple hops (guard crossings), our algorithm needs to
re-add guards to the candidate pool under certain conditions. In our current
implementation, when the candidate pool is empty, and the integration step has
just completed, we re-add any guards with negative-valued event functions. In
future implementations, we would like to generalize this behavior to a “guard
reset matrix,” which would allow us to specify how crossing one guard re-
enables other guards to become candidates. For the 1D hopper, crossing the
flight-stance guard indicates the trajectory may soon cross the stance-flight
guard and that it should be reset.

⬇ import numpy as np import favorite_integrator as dopri5 def integrate(f, h,
Dh, eps, x0, t0, tf): # Termination Condition def ntr_Nvk(x1, x2): return
((h(x1) < -eps) & (h(x2) >= -eps)).any() dopri5.term_early_func = ntr_Nvk #
Begin Integration x0 = np.asarray(y0) tq, xq = t0, x0.copy() tlog, xlog =
[tq], [xq.copy()] crossmsk = h(x0) >= 0 while tq < tf or not all(crossmsk): hx
= h(xq) if not any((hx > -eps)[$\sim$crossmsk]): tq, xq = dopri5(xq, tq, tf)
tq, xq = tq[1:], xq[1:] if all(crossmsk): crossmsk = h(xq[-1]) >= 0 else: #
Projection Step G, e = f(xq, hx), Dh(xq) eG = np.dot(e, G) dt = -hx / eG
dt[crossmsk] = inf j = argmin(dt) dt = dt[j] dx = G * dt tq = [tq + dt] xq =
[xq + dx] crossmsk[j] = True # End Loop / Next Iteration tlog.append(tq)
xlog.append(xq) tq, xq = tq[-1], xq[-1] return np.hstack(tlog),
np.vstack(xlog) Listing 1: Python implementation of our algorithm.
---

###  3.2 Differences from Council et al [13]

We derived our projection step from the algorithm presented by Council et al
in [13, Fig. 2.1]. Our projection step differs in where we sample the vector
field and event functions, and the assumptions we make about the trajectory’s
path through the event functions. A collection of $m$ event functions
$\\{h_{k}\\}$ for $k\in\\{1,\dots,m\\}$ which all take the value $0$ at a
point $\rho$ are called an “($m$-dimensional) ESS complex (crossing at
$\rho$).” Council’s algorithm assumes that the system has only one ESS
complex, that the B-Deriv and $Dh$ at $\rho$ are already known, and that the
coordinates have been transformed such that $\rho$ is the origin [13]. We
sought to relax these restrictions to support systems with multiple ESS
complexes and to minimize the additional modelling requirements beyond
defining a system’s event functions. To that end, we do not assume knowledge
of $\rho$ and instead sample $h$, $Dh$, and $f$ at a point in $N^{\nu_{k}}$
lying along the trajectory.

##  4 Order Test

Prior to proving properties of a numerical method, it is often useful to
evaluate the method on a toy problem with a known exact solution. Such an
evaluation is doubly useful both to build intuition for the method’s behavior
and to verify that the implementation lacks obvious errors. Numerical
integrators are commonly evaluated on the relationship between their step size
and normed integration error as bounded by an integer-power polynomial [21,
Chp II]. The degree of this polynomial minus one is called the integrator’s
“order” [21, Chp II]. Given an exact solution $y$, initial condition $y_{0}$,
initial time $t_{0}$, and step size $\delta t$, an integrator approximates
approximates the exact value $y(t_{0}+\delta t)$ as $y_{1}$. For an integrator
of order $p$ and some constant $K$, the error bound takes the form

$||y(t_{0}+\delta t)-y_{1}||\leq K\delta t^{p+1}$

[21, Eq 1.10, p.134].

Figure 3: Log-log plot of our algorithm’s RMS integration error versus
$1/\varepsilon$ (dotted red line) for an initial value problem. The data
interval used to determine the order (gray rectangle) and the fit extended
over the entire interval considered (dashed blue line) are overlaid. The
integrator achieved second-order accuracy, indicated by the slope of the blue
line. The flat tails for small and large $1/\varepsilon$ are expected.

We sought to empirically assess the order of our algorithm’s projection step
as a function of $\varepsilon$. Our integration step utilized the 5th-order
dopri5 integrator [21, II.4, p.178], and the related B-Deriv projection
algorithm developed by Council et al has been proven to be 2nd-order in
$\varepsilon$ [13]. As such, we expected our projection step to be at most
2nd-order. From this, we expected that any trajectory integration error should
be dominated by the error of our projection step, and not by our integration
step. For such an assessment, we needed a differential equation of sufficient
complexity which still had a known closed-form flow. Linear and affine
differential equations possess readily-computable exact solutions through the
matrix exponential, making them a natural choice. For simplicity, we chose our
guards to be the coordinate planes. We constructed a 3D system of affine
ordinary differential equations, described in Tab. 1

| $x<0$ | $x\geq 0$
---|---|---
$y<0$ | $-y+1$ | $-2y+1$
$y\geq 0$ | $y+1$ | $10x+1$
(a) $\dot{x}$

$x<0$ | $x\geq 0$
---|---
$x+1$ | $x/2+2$
$-x+1$ | $y+1$
(b) $\dot{y}$

$z\leq 0$ | $z>0$
---|---
$3z-1$ | $-z-1$
$3z-1$ | $-z-1$
(c) $\dot{z}$

Table 1: The 3D system of affine ODEs used for our order test. The $x$
velocity (a), $y$ velocity (b), and $z$ velocity (c) are grouped separately.

We integrated the trajectory with initial condition
$(x,y,z)=(-0.4,-0.15,0.3)$. This trajectory crossed all guards during the
interval $t\in(0,0.5)$. We computed integration error as the root-mean squared
error (RMS) of the integrated trajectory from the exact trajectory. We found
the RMS to be bounded by a polynomial of order at least $2.1$. From this
experiment, we expect our algorithm to be 2nd-order, as
$\mathrm{RMS}(\varepsilon)\leq K\varepsilon^{2.1}<K\varepsilon^{3}$, see Fig.
3.

##  5 Performance Comparison with MuJoCo

###  5.1 Calibration

To assess the practical merit of our algorithm, we sought to compare it with a
state-of-the-art simulation framework which explicitly treats collision
phenomena. We elected to use MuJoCo, short for “Multi-Joint dynamics with
Contact,” for our performance baseline [14]. We selected MuJoCo for the
simplicity of use with Python, and the ease of reading and writing modelling
parameters from its XML modelling framework. MuJoCo is a simulation package
which has received much attention for robot simulation [15, 49, 31, 45]. As
such, MuJoCo has benefitted from many years of dedicated optimization in C and
C++ [14], while we implemented our integrator in Python. It is well-known that
Python is less performant than C [22, 42], and it would be reasonable to not
only expect MuJoCo to outperform our algorithm in a head-on comparison of
execution time, but also that our algorithm’s execution time would be
dominated by the runtime performance of the Python virtual machine itself.

To establish a meaningful comparison between our algorithm and MuJoCo
uncomplicated by the details of software implementation, we needed an event-
selected toy robot model with a closed-form exact solution to calibrate both
the runtime performance and the accuracy of each solver. Inspired by the
hopping robot in Fig. 1, we wrote a basic model for a 1D hopping robot with no
damping:

$\dot{z}=\begin{cases}-g,&z>L_{0}\\\ -g+k(L_{0}-z)/m,&z\leq L_{0}\end{cases}$

with event functions $h_{1}(z)=z-z_{0},\;\;h_{2}(z)=-z+z_{0}$. $h_{2}$ may
seem redundant, as it has an identical guard to $h_{1}$, but it is easy to
verify that $h_{1}$ only satisfies requirement (4) when $\dot{z}<0$ and that
$h_{2}$ is needed to satisfy it when $\dot{z}>0$. We simulated an equivalent
hopper model in MuJoCo starting at rest with $z_{0}=2L_{0}=2$ from $t=0$ to
$t=2$. We configured the integration step in our integrator to have the same
time step ($\delta t=0.002$s) as MuJoCo’s integrator and then ran a bisection
search to find the $\varepsilon^{\prime}$ where our algorithm computed a
trajectory with equivalent RMS error from the exact solution as MuJoCo, shown
in Fig. 4

Figure 4: Time series output of our calibration: integrating a the trajectory of a 1D hopping robot. The top row of plots show the exact trajectory (dotted gray), MuJoCo’s result (dashed blue), and our result (red) overlaid. The bottom row shows the residual from the exact trajectory for both MuJoCo and our solver. The left plots show height above the ground (contact height $=1$) and the right plots show vertical velocity.  $\varepsilon^{\prime}$ | MuJoCo (s) | ours (s)
---|---|---
0.13 | 0.002436 | 0.010251
Table 2: Calibration parameters determined by simulating the same 1D hopper
initial value problem with MuJoCo and our algorithm. Our algorithm achieved
identical RMS integration error to MuJoCo at
$\varepsilon=\varepsilon^{\prime}$. The runtime in seconds of each method is
tabulated in the right two columns. Both solvers had maximum a time step of
0.002 seconds.

###  5.2 Experiment

Our algorithm’s potential gains are in its treatment of multiple guards. We
wanted to assess how its performance scaled versus MuJoCo with increasing
number of contacts. We developed a model of a plate in SE(2) falling due to
gravity onto a bed of evenly-spaced vertical springs with frictionless rolling
contacts. A model of this form allowed us to automatically generate models and
MuJoCo compatible XML of the plate falling onto arbitrarily many spring
contacts. We assumed the plate was long enough that it was above each spring
at all times. For the $i$th spring with x-position $x_{s_{i}}$, our flight-
stance event functions and Jacobians took the form:

$\displaystyle h_{i}$ $\displaystyle=L_{0}+\tan\theta\cdot(x-x_{s_{i}})-z$
$\displaystyle Dh_{i}$
$\displaystyle=\begin{bmatrix}\tan\theta,&-1,&sec^{2}\theta\cdot(x-x_{s_{i}})\end{bmatrix}$

for $i\in\\{1,\dots,m/2\\}$, with the stance-flight event functions and
Jacobians being their negative as with 1D hopper.

We generated 120 random initial conditions starting from rest centered over
the spring bed with $z_{0}\in(2,3)$ and $\theta_{0}\in(-0.1,0.1)$. For each
initial condition, we simulated the plate falling onto spring beds of 2 to 100
evenly-spaced identical spring-damper contacts with the precision parameter
determined from our calibration. For a large number of contacts (roughly 150
contacts), both MuJoCo and our algorithm started to exhibit irregular
behavior, leading us to believe their performance above 100 contacts was an
artifact of implementation, and not reflective of their algorithmic
performance. We evaluated each algorithm on its marginal runtime cost per
contact. We computed marginal cost as the derivative of the median runtime
with respect to the number of contacts. We found that the marginal runtime
cost per contact for our Python implementation was roughly the same as
MuJoCo’s, see Fig. 5.

Figure 5: The marginal cost of our algorithm (red) and MuJoCo (dashed blue)
for increasing number of spring-damper contacts. We integrated the same
$N=120$ random initial conditions for each number of contacts. Our algorithm
implemented in Python has roughly the same marginal cost per contact as
MuJoCo.

##  6 Proof of First-Order Correctness

This “Event Selected System” [8, Defn. 2] has a unique maximal flow [8, Cor.
5] $\Phi:\mathcal{F}\to\mathfrak{X}$ on a maximal flow domain
$\mathcal{F}\subseteq\mathbb{R}\times\mathfrak{X}$, and $(0,x)\in\mathcal{F}$
for all $x$.

We denote the first parameter of $\Phi$ with a superscript, and take
$\Phi^{t}$ to be an operator acting on a state to its right. Note the “flow
property”: $\Phi^{t}\Phi^{s}x=\Phi^{t+s}x$ whenever
$\\{(t+s,x),(s,x),(t,\Phi^{s}x)\\}\in\mathcal{F}$.

Let $\tau_{k}:N^{\nu_{k}}\to[-\nu_{k}/c,\nu_{k}/c]$ be the unique arrival time
(in future or past) $\tau_{k}(x)$ of initial condition $x\in N^{\nu_{k}}$ at
the guard $\mathfrak{g}_{k}$. This implies the existence of the impact point
map $\pi_{k}(x):=\Phi^{\tau_{k}(x)}x$,
$\pi_{k}:N^{\nu_{k}}\to\mathfrak{g}_{k}$, which is a non-linear projection,
i.e. $\pi_{k}\circ\pi_{k}=\pi_{k}$, $\pi_{k}|_{\mathfrak{g}_{k}}=\text{Id}$.

The correctness of our integration scheme builds on two parts: (i) a condition
for the validity of projection onto the guards, which can be used for speeding
up the computation of impacts on co-dimension one guard manifolds in any
smooth dynamical system; and (ii) the guarantee of ESS, arising from liveness,
that all transitions resolve into a unique, continuous and piecewise smooth
flow [8, Thm. 4].

###  6.1 Efficient projection onto guards

First, we derive a bound on the impact time for reaching all guards. Consider
an execution starting with $x\in N^{\varepsilon}$,
$\varepsilon\leq\min_{k}\nu_{k}$ and the arrival time $T:=\tau_{k}(x)$ for
reaching $\mathfrak{g}_{k}$ from a starting $x$ with $h_{k}(x)<0$. Observe
that

$\displaystyle-h_{k}(x)$
$\displaystyle=h_{k}(\Phi^{T}x)-h_{k}(x)=\int_{0}^{T}\frac{d}{dt}h_{k}(\Phi^{t}x)dt$
$\displaystyle=\int_{0}^{T}\nabla h_{k}(\Phi^{t}x)\cdot
f(\Phi^{t}x,h(\Phi^{t}x))dt\geq cT,$

i.e. the arrival time $\tau_{k}(x)\leq-h_{k}(x)/c\leq\varepsilon/c$, despite
any possible domain switching induces by guard crossings.

Consider the execution starting at $x$ with some component $h_{k}(x)<0$, and
follow it until its earliest guard crossing, where it crosses
$\mathfrak{g}_{j}$. The smoothness of each of the (at most) $2^{m}$ pieces of
$f(\cdot,\cdot)$ together the compactness of $\mathfrak{X}$ guarantee that $f$
is Lipschitz on $\mathfrak{X}$ with a Lipschitz constant $L$, and bounded by
some bound $M$. Because no guards are being crossed, we can define
$f_{b}(x):=f(x,h(x))$ the smooth vector field that applies until that first
guard is crossed.

From the commonly known [28, Thm. 3.4 p.96] we have:

$\displaystyle\|f_{b}(x)t+x-\Phi^{t}x\|\leq$
$\displaystyle\|x-\Phi^{0}x\|\mathrm{e}^{Lt}$
$\displaystyle+\frac{\mathrm{e}^{Lt}-1}{L}\|f_{b}(y)-f_{b}(x)\|_{\infty}$
$\displaystyle\leq$
$\displaystyle\frac{2M}{L}\left(\mathrm{e}^{-h_{j}(x)L/c}-1\right).$

Since the bound approaches zero at least linearly fast when
$\varepsilon>-h_{j}(x)$ does, it follows that for any $\delta>0$, we select
$\varepsilon>0$ such that $\|f_{b}(x)t+x-\Phi^{t}x\|<\delta$ for all
$t\leq\varepsilon/c$. Define the first order approximate impact time
$\hat{\tau}(x):=-h_{j}(x)/(\nabla h(x)\cdot f_{b}(x))$, and the first order
approximate impact map $\hat{\pi}_{j}(x):=f_{b}(x)\hat{\tau}(x)+x$. We wish to
bound the error $\|\hat{\pi}_{j}(x)-\pi_{j}(x)\|$.

Note $\hat{\tau}(x)\leq-h_{j}(x)/c\leq\varepsilon/c$, and therefore
$\|\hat{\pi}(x)-\Phi^{\hat{\tau}(x)}x\|<\delta$. For the true impact time we
have
$\|f_{b}(x)\tau(x)+x-\Phi^{\tau(x)}x\|=\|f_{b}(x)\tau(x)+x-\pi_{j}(x)x\|$. But
also
$\|(f_{b}(x)\tau(x)+x)-\hat{\pi}_{j}(x)\|=\|f_{b}(x)\||\hat{\tau}(x)-\tau(x)|<2M\varepsilon/c$.
Thus from the triangle inequality
$\|\hat{\pi}_{j}(x)-\pi_{j}(x)\|<\delta+2M\varepsilon/c$, which is order
$\varepsilon$ small.

From a practical standpoint it is more convenient to give a $\varepsilon_{k}$
for each component $h_{k}(\cdot)$, such that the worst case impact point error
$\|\hat{\pi}_{k}-\pi_{k}\|_{\infty}$ is sufficiently small for the simulation
needed.

The elementary calculations above provide nothing essentially novel about
computing first order impact approximations compared to what is already known
to the applied math community, but the approach of employing guard functions
to also provide an algorithmic hint as to whether projections are appropriate
to employ was not one we saw elsewhere. The value of our work here is that it
provided a general purpose integrator that applies to ESS, which model a broad
class of multi-contact problems which appear in robotics.

##  7 Discussion

In this work we have been primarily interested in presenting our algorithm and
its applicability to treating multiple contacts in robotics. We have provided
an example implementation, a simple proof of correctness to first-order,
empirically demonstrated that it can achieve second-order accuracy, and found
that, with increasing number of contacts, its decrease in performance (as
implemented in Python), is no worse than MuJoCo’s decrease in performance. It
would be interesting to see how our algorithm performs on more realistic robot
models, and separately against additional simulation suites like PyBullet [12]
or the Open Dynamics Engine [44]. MuJoCo is known to simulate robots faster
than most other simulation suites, so we would expect our algorithm to out-
perform them as well [15]. This makes our algorithm particularly well-suited
to real-time applications, such as online optimization and model-predictive
control of robots with 6 legs or more. It may now be tractable to design robot
manipulators with more than 3 fingers, with real-time trajectory optimization
enabled by our method. As state estimators are designed to be robust to noise,
the speed-accuracy tradeoff of our algorithm could be more heavily exploited
in contact-aware state estimation [23, 39].

It is worth noting that ESS is more general than the types of hybrid problems
MuJoCo and PyBullet treat [45, 12, 8], and our algorithm could find uses
outside of robotics, such as in stability analysis of power grids [24],
economics [6], or in fields as far removed from engineering as ecology and
neurobiology [8, 6].

## References

  * [1] R. Altendorfer, N. Moore, H. Komsuoglu, M. Buehler, H. B. Brown, D. McMordie, U. Saranli, R. Full, and D. E. Koditschek, Rhex: A biologically inspired hexapod runner, Autonomous Robots, 11 (2001), pp. 207–213, https://doi.org/10.1023/A:1012426720699.
  * [2] M. Anderson, T. McLaughlin, and S. Revzen, Serial elastic versus parallel elastic actuators in hopping, in Dynamic Walking, no. 513, Internet, June 2020, Dynamic Walking.  URL.
  * [3] M. Anitescu and F. A. Potra, Formulating dynamic multi-rigid-body contact problems with friction as solvable linear complementarity problems, Nonlinear Dynamics, 14 (1997), pp. 231–247, https://doi.org/10.1023/A:1008292328909.
  * [4] J. E. Bares and D. S. Wettergreen, Dante ii: Technical description, results, and lessons learned, The International Journal of Robotics Research, 18 (1999), pp. 621–649, https://doi.org/10.1177/02783649922066475.
  * [5] A. Bicchi, Hands for dexterous manipulation and robust grasping: a difficult road toward simplicity, IEEE Transactions on Robotics and Automation, 16 (2000), pp. 652–662, https://doi.org/10.1109/70.897777.
  * [6] B. Brogliato, Some perspectives on the analysis and control of complementarity systems, IEEE Transactions on Automatic Control, 48 (2003), pp. 918–935, https://doi.org/10.1109/TAC.2003.812777.
  * [7] B. Brogliato, A. Daniilidis, C. Lemaréchal, and V. Acary, On the equivalence between complementarity systems, projected systems and differential inclusions, Systems & Control Letters, 55 (2006), pp. 45–51, https://doi.org/10.1016/j.sysconle.2005.04.015.
  * [8] S. A. Burden, S. S. Sastry, D. E. Koditschek, and S. Revzen, Event-selected vector field discontinuities yield piecewise-differentiable flows, SIAM Journal of Applied Dynamical Systems, 15 (2016), pp. 1227–1267, https://doi.org/10.1137/15M1016588.
  * [9] R. Cottle, G. Habetler, and C. Lemke, On classes of copositive matrices, Linear Algebra and its Applications, 3 (1970), pp. 295–310, https://doi.org/10.1016/0024-3795(70)90002-9.
  * [10] R. W. Cottle, Numerical methods for complementarity problems in engineering and applied science, in Computing Methods in Applied Sciences and Engineering, 1977, I, R. Glowinski, J. L. Lions, and I. Laboria, eds., Berlin, Heidelberg, 1979, Springer Berlin Heidelberg, pp. 37–52.
  * [11] R. W. Cottle, J.-S. Pang, and R. E. Stone, The Linear Complementarity Problem, Society for Industrial and Applied Mathematics, 2009, https://doi.org/10.1137/1.9780898719000.
  * [12] E. Coumans and Y. Bai, Pybullet, a python module for physics simulation for games, robotics and machine learning.  http://pybullet.org, 2016–2021.
  * [13] G. Council, S. Revzen, and S. A. Burden, Representing and Computing the B-Derivative of the Piecewise-Differentiable Flow of a Class of Nonsmooth Vector Fields, Journal of Computational and Nonlinear Dynamics, 17 (2022), p. 091004, https://doi.org/10.1115/1.4054481.
  * [14] Deepmind Technologies, Mujoco documentation, https://mujoco.readthedocs.io/en/stable/overview.html.
  * [15] T. Erez, Y. Tassa, and E. Todorov, Simulation tools for model-based robotics: Comparison of bullet, havok, mujoco, ode and physx, in 2015 IEEE International Conference on Robotics and Automation (ICRA), 2015, pp. 4397–4404, https://doi.org/10.1109/ICRA.2015.7139807.
  * [16] N. Fazeli, R. Kolbert, R. Tedrake, and A. Rodriguez, Parameter and contact force estimation of planar rigid-bodies undergoing frictional contact, The International Journal of Robotics Research, 36 (2017), pp. 1437–1454, https://doi.org/10.1177/0278364917698749.
  * [17] A. F. Filippov, Existence and General Properties of Solutions of Discontinuous Systems, Springer Netherlands, Dordrecht, 1988, pp. 48–122, https://doi.org/10.1007/978-94-015-7793-9_3.
  * [18] R. Goebel, R. G. Sanfelice, and A. R. Teel, Hybrid dynamical systems, IEEE Control Systems Magazine, 29 (2009), pp. 28–93, https://doi.org/10.1109/MCS.2008.931718.
  * [19] R. Goebel, R. G. Sanfelice, and A. R. Teel, Hybrid dynamical systems, IEEE Control Systems Magazine, 29 (2009), pp. 28–93, https://doi.org/10.1109/MCS.2008.931718.
  * [20] R. D. Gregg and M. W. Spong, Bringing the compass-gait bipedal walker to three dimensions, in 2009 IEEE/RSJ International Conference on Intelligent Robots and Systems, 2009, pp. 4469–4474, https://doi.org/10.1109/IROS.2009.5354583.
  * [21] E. Hairer, S. P. Nørsett, and G. Wanner, Solving ordinary differential equations I: nonstiff problems, vol. 8 of Springer Series in Computational Mathematics, Springer-Verlag, Berlin, 2 ed., 1993, https://doi.org/10.1007/978-3-540-78862-1.
  * [22] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant, Array programming with NumPy, Nature, 585 (2020), pp. 357–362, https://doi.org/10.1038/s41586-020-2649-2.
  * [23] R. Hartley, M. Ghaffari, R. M. Eustice, and J. W. Grizzle, Contact-aided invariant extended kalman filtering for robot state estimation, The International Journal of Robotics Research, 39 (2020), pp. 402–430, https://doi.org/10.1177/0278364919894385.
  * [24] I. Hiskens and M. Pai, Trajectory sensitivity analysis of hybrid systems, IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications, 47 (2000), pp. 204–220, https://doi.org/10.1109/81.828574.
  * [25] A. Ivanov, The stability of periodic solutions of discontinuous systems that intersect several surfaces of discontinuity, Journal of Applied Mathematics and Mechanics, 62 (1998), pp. 677–685, https://doi.org/10.1016/S0021-8928(98)00087-2.
  * [26] S. Jacobsen, J. Wood, D. Knutti, and K. Biggers, The utah/m.i.t. dextrous hand: Work in progress, The International Journal of Robotics Research, 3 (1984), pp. 21–50, https://doi.org/10.1177/027836498400300402.
  * [27] W. Jongkind, Dextrous gripping in a hazardous environment: guidelines, fault tolerance and control, in Proceedings of IEEE Systems Man and Cybernetics Conference - SMC, vol. 1, 1993, pp. 509–514 vol.1, https://doi.org/10.1109/ICSMC.1993.384794.
  * [28] H. K. Khalil, Nonlinear systems, Prentice Hall, 3rd ed., 2002.
  * [29] S. Kim, J. E. Clark, and M. R. Cutkosky, isprawl: Autonomy, and the effects of power transmission, in Climbing and Walking Robots, Berlin, Heidelberg, 2005, Springer Berlin Heidelberg, pp. 859–867.
  * [30] M. Kiseleva, N. Kuznetsov, and G. Leonov, Theory of Differential Inclusions and Its Application in Mechanics, Springer International Publishing, Cham, 2018, pp. 219–239, https://doi.org/10.1007/978-3-319-62464-8_9.
  * [31] J. Koenemann, A. Del Prete, Y. Tassa, E. Todorov, O. Stasse, M. Bennewitz, and N. Mansard, Whole-body model-predictive control applied to the hrp-2 humanoid, in 2015 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), 2015, pp. 3346–3351, https://doi.org/10.1109/IROS.2015.7353843.
  * [32] J. Lee, M. X. Grey, S. Ha, T. Kunz, S. Jain, Y. Ye, S. S. Srinivasa, M. Stilman, and C. K. Liu, Dart: Dynamic animation and robotics toolkit, Journal of Open Source Software, 3 (2018), p. 500, https://doi.org/10.21105/joss.00500.
  * [33] Y. Li and I. Kao, A review of modeling of soft-contact fingers and stiffness control for dextrous manipulation in robotics, in Proceedings 2001 ICRA. IEEE International Conference on Robotics and Automation (Cat. No.01CH37164), vol. 3, 2001, pp. 3055–3060 vol.3, https://doi.org/10.1109/ROBOT.2001.933086.
  * [34] P. Lötstedt, Mechanical systems of rigid bodies subject to unilateral constraints, SIAM Journal on Applied Mathematics, 42 (1982), pp. 281–296, https://doi.org/10.1137/0142022.
  * [35] J. Lygeros, K. Johansson, S. Simic, J. Zhang, and S. Sastry, Continuity and invariance in hybrid automata, in Proceedings of the 40th IEEE Conference on Decision and Control (Cat. No.01CH37228), vol. 1, 2001, pp. 340–345 vol.1, https://doi.org/10.1109/CDC.2001.980123.
  * [36] T. McLaughlin, M. Anderson, and S. Revzen, Hopping with elastic restitution is more difficult than it seems, American Physical Society, 2020, http://meetings.aps.org/Meeting/MAR20/Session/S22.11.
  * [37] J. J. Moreau, Unilateral Contact and Dry Friction in Finite Freedom Dynamics, Springer Vienna, Vienna, 1988, pp. 1–82, https://doi.org/10.1007/978-3-7091-2624-0_1.
  * [38] T. Okada, Object-handling system for manual industry, IEEE Transactions on Systems, Man, and Cybernetics, 9 (1979), pp. 79–89, https://doi.org/10.1109/TSMC.1979.4310152.
  * [39] J. J. Payne, N. J. Kong, and A. M. Johnson, The uncertainty aware salted kalman filter: State estimation for hybrid systems with uncertain guards, in 2022 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), 2022, pp. 8821–8828, https://doi.org/10.1109/IROS47612.2022.9981218.
  * [40] M. Posa, C. Cantu, and R. Tedrake, A direct method for trajectory optimization of rigid bodies through contact, The International Journal of Robotics Research, 33 (2014), pp. 69–81, https://doi.org/10.1177/0278364913506757.
  * [41] M. Posa, M. Tobenkin, and R. Tedrake, Stability analysis and control of rigid-body systems with impacts and friction, IEEE Transactions on Automatic Control, 61 (2016), pp. 1423–1437, https://doi.org/10.1109/TAC.2015.2459151.
  * [42] L. Prechelt, An empirical comparison of seven programming languages, Computer, 33 (2000), pp. 23–29, https://doi.org/10.1109/2.876288.
  * [43] M. Raff and C. D. Remy, Hard vs. soft contacts in trajectory optimization of legged robots: An example with a compliant rubber foot, 2023, https://api.semanticscholar.org/CorpusID:259305317.
  * [44] R. L. Smith, Open dynamics engine.  http://ode.org, 2001–2004.
  * [45] E. Todorov, T. Erez, and Y. Tassa, Mujoco: A physics engine for model-based control, in 2012 IEEE/RSJ International Conference on Intelligent Robots and Systems, 2012, pp. 5026–5033, https://doi.org/10.1109/IROS.2012.6386109.
  * [46] E. Wendel and A. D. Ames, Rank deficiency and superstability of hybrid systems, Nonlinear Analysis: Hybrid Systems, 6 (2012), pp. 787–805, https://doi.org/10.1016/j.nahs.2011.09.002.
  * [47] E. Westervelt, J. Grizzle, and D. Koditschek, Hybrid zero dynamics of planar biped walkers, IEEE Trans Automat Contr, 48 (2003), pp. 42–56, https://doi.org/10.1109/TAC.2002.806653.
  * [48] G. Wohlke, The karlsruhe dextrous hand: grasp planning, programming and real-time control, in Proceedings of IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS’94), vol. 1, 1994, pp. 352–359 vol.1, https://doi.org/10.1109/IROS.1994.407370.
  * [49] Z. Zhang, J. Chen, Z. Chen, and W. Li, Asynchronous episodic deep deterministic policy gradient: Toward continuous control in computationally complex environments, IEEE Transactions on Cybernetics, 51 (2021), pp. 604–613, https://doi.org/10.1109/TCYB.2019.2939174.


# Constraint on the minimally extended varying speed of light using time
dilations in type Ia supernovae

Seokcheon Lee 1
1Department of Physics, Institute of Basic Science, Sungkyunkwan University,
Suwon 16419, Korea E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The Friedmann-Lemaître-Robertson-Walker model establishes the correlation
between redshifts and distances. It has a metric expansion of space. As a
result, the wavelength of photons propagating through the expanding space is
stretched, creating the cosmological redshift, $z$. It also relates the
frequency of light detected by a local observer to that emitted from a distant
source. In standard cosmology (i.e., a constant speed light model), this
relation is given by a factor $1/(1+z)$. However, this ratio is modified in
the minimally extended varying speed of light model (meVSL, $c=c_{0}a^{b/4}$)
as $1/(1+z)^{1-b/4}$. This time dilation effect is detected as the observed
rate of the time variation in the intensity of emitted radiation. The spectra
of type Ia supernovae (SNe Ia) provide a reliable way to measure the apparent
aging rate of distant objects. We use data on 13 high-redshift ($0.28\leq
z\leq 0.62$) SNe Ia to obtain $b=0.198\pm 0.415$ at the $1$-$\sigma$
confidence interval. The current data is too sparse to give meaningful
constrain on the meVSL and cannot distinguish the meVSL model from the
standard model.

###### keywords:

Varying speed of light – Time dilation – Type Ia Supernovae

††pubyear: 2023††pagerange: Constraint on the minimally extended varying speed
of light using time dilations in type Ia supernovae–Constraint on the
minimally extended varying speed of light using time dilations in type Ia
supernovae

##  1 Introduction

The spectrum of light coming from a distant source shows features in it, such
as absorption lines, emission lines, or other variations in light intensity,
and one can use them to determine the (cosmological) redshift Harrison (2000).
Redshift characterizes the relative difference between observed and emitted
wavelengths of an object. Friedmann’s solutions to Einstein’s equations of
general relativity (GR) by adopting the Robertson-Walker (RW) metric show that
the wavelength of photons propagating through the expanding space is
stretched, creating the redshift. Thus, the redshift, $z$, is the fundamental
observable in the Friedmann-Lemaître-Robertson-Walker (FLRW) model.

The stretching of the wavelength of photons in the expanding universe can be
related to the shrinking of the frequency by using the dispersion relation of
photons. Consequently, this effect should be detected as the time dilation
(TD) in the intensity of emitted radiation Ishak (2019).

The type Ia supernovae (SNe Ia) light evolutions have been used as cosmic
clocks to measure TD in distant objects due to their large luminosities and
variability on short timescales. Both light curves (LCs) of distant SNe Ia
Leibundgut et al. (1996); Goldhaber et al. (2001) and their spectral
evolutions Blondin et al. (2008); Riess et al. (1997); Foley et al. (2005)
provide methods to measure the apparent aging rate of faraway objects.
Especially the latter removes the degeneracy between the intrinsic LC width
and TD effect to make it more reliable.

In standard cosmology (SC i.e., $c=$ constant), the ratio of the frequency of
light detected by a local observer to that emitted from a distant source is
given by a factor $1/(1+z)$ Weinberg (1972). However, this ratio is modified
in the minimally extended varying speed of light model (meVSL,
$c=c_{0}a^{b/4}$) as $1/(1+z)^{1-b/4}$ Lee (2021, 2023a, 2023b). Thus, one can
use the TD effect detected in SNe Ia to distinguish these models.

We organize this paper as follows. In section 2, we review the TD formula of
the meVSL model. In section 3, we analyze it with 13 high-redshift ($0.28\leq
z\leq 0.62$) SNe Ia Blondin et al. (2008) data to constrain the meVSL model.
We give our conclusion in section 4.

##  2 Review of cosmological redshift

In the expanding universe, one can obtain it by defining the proper distance
$d_{\textrm{p}}$ as the multiplication of the comoving distance $\chi$ by the
scale factor $d_{\textrm{p}}\equiv a(t)\chi$. Comoving distance is the
distance between two points measured along a path defined at the present
cosmological time. It remains constant in time for objects moving with the
Hubble flow. From this definition, one obtains the relation between proper
distances at $t_{1}$ and $t_{2}$ as

$\displaystyle\chi=\frac{d_{\textrm{p}}(t_{1})}{a(t_{1})}=\frac{d_{\textrm{p}}(t_{2})}{a(t_{2})}\quad\Rightarrow\quad
d_{\textrm{p}}(t_{1})=\frac{a(t_{1})}{a(t_{2})}d_{\textrm{p}}(t_{2})\,.$ (1)

Now one can replace $d_{\textrm{p}}(t_{1})$ as the measured wavelength at
present ($a(t_{0})=1$) and $d_{\textrm{p}}(t_{2})$ as the emitted one in the
past ($a(t_{\textrm{emit}})$) to find their relationship with redshift, $z$ as

$\displaystyle\lambda_{\textrm{obs}}=\frac{1}{a(t_{\textrm{emit}})}\lambda_{\textrm{emit}}\equiv(1+z)\lambda_{\textrm{emit}}\quad\textrm{where}\quad
z\equiv\frac{\lambda_{\textrm{obs}}-\lambda_{\textrm{emit}}}{\lambda_{\textrm{emit}}}\,.$
(2)

One can relate redshift $z$ to the frequency in the RW model given by

$\displaystyle ds^{2}$ $\displaystyle=g_{\mu\nu}dx^{\mu}dx^{\nu}$
$\displaystyle=-\left(dx^{0}\right)^{2}+a(t)^{2}\left(\frac{dr^{2}}{1-kr^{2}}+r^{2}\left(d\theta^{2}+\sin^{2}\theta
d\phi^{2}\right)\right)$
$\displaystyle\equiv-\left(dx^{0}\right)^{2}+a(t)^{2}\left(d\chi^{2}+r^{2}d\Omega^{2}\right)\,,$
(3)

where $x^{0}=ct$, $r$ is the coordinate distance, and $d\Omega$ is the solid
angle. To derive the redshift, one uses the geodesic equation for a light wave
$ds^{2}=0$. Let us assume that an observer at a position $r=0$ and time
$t=t_{0}$ measures the first crest of a light wave emitted at a position $r=R$
and time $t=t_{\textrm{emit}}$. At the time $t_{0}+\delta t_{0}$, the observer
sees the next crest of the same light emitted at a time
$t_{\textrm{emit}}+\delta t_{\textrm{emit}}$. Then, one obtains the following
equation from the RW metric

$\displaystyle\chi$
$\displaystyle=\int_{0}^{R}\frac{dr}{\sqrt{1-kr^{2}}}=\int_{t_{\textrm{emit}}}^{t_{0}}\frac{dx^{0}}{a(t)}=\int_{t_{\textrm{emit}}+\delta
t_{\textrm{emit}}}^{t_{0}+\delta t_{0}}\frac{dx^{0}}{a(t)}$
$\displaystyle\Rightarrow\int_{t_{0}}^{t_{0}+\delta
t_{0}}\frac{dx^{0}}{a(t)}=\int_{t_{\textrm{emit}}}^{t_{\textrm{emit}}+\delta
t_{\textrm{emit}}}\frac{dx^{0}}{a(t)}\,.$ (4)

For tiny variations in time (throughout one cycle of a light wave), scale
factors are essentially constants (i.e., $a(t_{\textrm{emit}}+\delta
t_{\textrm{emit}})\approx a(t_{\textrm{emit}})$ and $a(t_{0}+\delta
t_{0})\approx a(t_{0})$). It yields

$\displaystyle\frac{1}{a(t_{0})}\Bigl{[}x_{0}(t_{0}+\delta
t_{0})-x_{0}(t_{0})\Bigr{]}$
$\displaystyle=\frac{1}{a(t_{\textrm{emit}})}\Bigl{[}x_{0}(t_{\textrm{emit}}+\delta
t_{\textrm{emit}})-x_{0}(t_{\textrm{emit}})\Bigr{]}\,.$ (5)

In SC, $x_{0}(t_{0}+\delta t_{0})-x_{0}(t_{0})=c\delta t_{0}$ and
$x_{0}(t_{\textrm{emit}}+\delta
t_{\textrm{emit}})-x_{0}(t_{\textrm{emit}})=c\delta t_{\textrm{emit}}$ to make
Eq. (5) as

$\displaystyle\frac{c\delta t_{0}}{a(t_{0})}=\frac{c\delta
t_{\textrm{emit}}}{a(t_{\textrm{emit}})}\quad\Rightarrow\quad\lambda(t_{0})=\frac{1}{a(t_{\textrm{emit}})}\lambda(t_{\textrm{emit}})=(1+z)\lambda_{\textrm{emit}}\,,$
(6)

where we define the frequency as $1/\delta
t_{i}=\nu(t_{i})=\nu(t_{0})a(t_{i})^{-1}$ and use the relation
$c=\lambda(t_{i})\nu(t_{i})=$ constant Weinberg (1972). Eq. (6) is identical
to Eq. (2). In the meVSL model (i.e., $c(t_{i})=c_{0}a(t_{i})^{b/4}$) and
$x_{0}(t_{0}+\delta t_{0})-x_{0}(t_{0})\approx c_{0}a(t_{0})^{b/4}\delta
t_{0}$, thus $x_{0}(t_{\textrm{emit}}+\delta
t_{\textrm{emit}})-x_{0}(t_{\textrm{emit}})\approx
c_{0}a(t_{\textrm{emit}})^{b/4}\delta t_{\textrm{emit}}$. From this fact, Eq.
(5) becomes

$\displaystyle\frac{c_{0}\delta t_{0}}{a(t_{0})^{1-b/4}}=\frac{c_{0}\delta
t_{\textrm{emit}}}{a(t_{\textrm{emit}})^{1-b/4}}$
$\displaystyle\Rightarrow\lambda(t_{0})=\frac{1}{a(t_{\textrm{emit}})}\lambda(t_{\textrm{emit}})=(1+z)\lambda_{\textrm{emit}}\,,$
(7)

where we use the frequency as $1/\delta
t_{i}=\nu(t_{i})=\nu(t_{0})a(t_{i})^{-1+b/4}$ and use the relation
$c(t_{i})=\lambda(t_{i})\nu(t_{i})$ Barrow & Magueijo (2000); Balcerzak &
Dąbrowski (2014); Lee (2021, 2023a). Again one can recover the redshift
defined in the expanding universe even in the meVSL model shown in equation
(7). Also, some other VSL models have provided the same equation for the
redshift of the radiation wavelength, $\lambda$, and the scale factor, $a$, as
in SC Gupta (2020); Cuzinatto et al. (2023).

##  3 Analysis

In an expanding universe, the galaxies are all moving away from each other,
and clocks run at different rates to produce the TD in a distant galaxy. In
the RW metric, the space is homogeneous and isotropic, and the expansion is so
regular that the world lines of the galaxies form a nonintersecting and
diverging three-bundle of geodesics. One can define spacelike hypersurfaces
that are orthogonal to all of them, and one can label each hypersurface by a
constant value of the global cosmic time coordinate $t$ Islam (2001); Roos
(2015).

However, the rate of local clocks depends on the hypothesis of the constancy
of the speed of light as shown in the previous section 2. Thus, it is required
to compare models with current TD data to clarify the underlying model.
However, there are few available observational data due to the difficulty of
the direct measurement of TD and its independence of cosmological parameters.
Fortunately, LCs of distant SNe Ia provide data for TD Leibundgut et al.
(1996); Goldhaber et al. (2001), and more reliable data is available from
spectral evolutions of a high redshift SNe Ia Riess et al. (1997); Foley et
al. (2005); Blondin et al. (2008).

###  3.1 Brief review of SNe as clocks

The LC of an SN tracks its brightness over time, showing how its luminosity
changes as it evolves. It begins with a rise in brightness called the "pre-
maximum" phase, reaching peak luminosity, followed by a decline. The shape of
the LC reveals details about the SN, especially for SNe Ia used as standard
candles in cosmology. By studying the LC, astronomers determine peak
brightness, time to reach maximum brightness, and decline rate, aiding the
classification and understanding of energetics, composition, and explosion
mechanisms. Comparing LCs at different distances allows investigation of
cosmic expansion and TD, leading to discoveries like accelerating universe
expansion and dark energy.

Wilson’s method, proposed in 1939, compares the LCs of nearby and distant SNe
Wilson (1939). Those of distant SNe exhibit TD, appearing stretched compared
to nearby ones due to the travel time of light through space. By analyzing the
time dilation effect in LCs, one can determine the time it took for light to
reach the observer from SNe at different distances. This information helps
study the expansion rate of the Universe and test cosmological models.
Practically, this involves gathering data on supernova brightness and
evolution at various redshifts, fitting mathematical models to their LCs, and
comparing observed TD with model predictions.

However, in SNe Ia, there is an inherent correlation between the breadths of
their LCs and their peak luminosities. Brighter SNe Ia tend to have broader
LCs Phillips (1993). It can lead to confusion when distinguishing between the
effects of time dilation and variations in luminosity, especially if the
latter is perceived as a stretching of the LC.

The temporal evolution of SN Ia spectra offers a more reliable approach to
measuring the passage of time for individual SN. One can use the wavelengths
of permanent spectral features to determine the timing of a supernova’s peak
brightness Minkowski (1939). The evolution of spectral characteristics in most
SNe Ia follows a consistent pattern, allowing for an objective assessment of a
supernova’s age based on its spectrum. Comparing two or more spectra of the
same distant supernova to those of nearby SNe measures the rate at which the
distant supernova ages. This method of measuring the aging rate provides a
direct means to verify the presence of time dilation expected in an expanding
Universe Riess et al. (1997).

SN | $z$ | $1/(1+z)$ | aging rate (error)
---|---|---|---
1996bj | 0.574 | 0.635 | 0.527 (0.369)
1997ex | 0.361 | 0.735 | 0.745 (0.076)
2001go | 0.552 | 0.644 | 0.652 (0.062)
2002iz | 0.427 | 0.701 | 0.655 (0.089)
b027 | 0.315 | 0.760 | 0.823 (0.092)
2003js | 0.363 | 0.734 | 0.718 (0.082)
04D2an | 0.621 | 0.617 | 0.567 (0.341)

SN | $z$ | $1/(1+z)$ | aging rate (error)
---|---|---|---
2006mk | 0.475 | 0.678 | 0.753 (0.060)
2006sc | 0.357 | 0.737 | 0.619 (0.121)
2006tk | 0.312 | 0.762 | 0.835 (0.181)
2007tg | 0.502 | 0.666 | 0.687 (0.102)
2007tt | 0.374 | 0.728 | 0.718 (0.108)
2007un | 0.283 | 0.779 | 0.759 (0.135)

Table 1: These are aging rate measurements from SNe Ia observation Blondin et
al. (2008).

###  3.2 The minimum $\chi^{2}$ analysis

We show aging rate measurements of $13$ high-redshift SNe Ia data provided in
the reference Blondin et al. (2008). They are given in Tab. 1. They are based
on a sample of $35$ spectra of $13$ SNe Ia in the redshift range $0.28\leq
z\leq 0.62$ by comparing their differences between the observer-frame ages and
the rest-frame ones.

Model |  $\chi^{2}$/dof | $\bar{b}$ |  $1$-$\sigma$ | GoF (%)
---|---|---|---|---
$(1+z)^{-1}$ | $3.6/13$ | $0$ | $0$ | $99.5$
$(1+z)^{-1+b/4}$ | $3.4/12$ | $0.198$ | $0.415$ | $99.2$
Table 2: A simple $\chi^{2}$ analysis for both the SC and the meVSL model
using TD data Blondin et al. (2008). The best-fit value for the $b$ exponent
is $b=0.198$ with $0.415$ for the $1$-$\sigma$ confidence interval.

We perform a least-square fit to the data for predictions of both the SC and
the meVSL model

$\displaystyle\chi^{2}=\sum_{i=1}^{13}\frac{\left((1+z_{i})_{\textrm{obs}}^{-1}-(1+z_{i})^{-1+b/4}\right)^{2}}{\sigma_{i}^{2}}\,,$
(8)

where we use Eq. (7), $(1+z_{i})_{\textrm{obs}}^{-1}$ denote 13 data points,
and $\sigma_{i}$ represent their errors in Table 1. The SC aging rate
$1/(1+z)$ yields a fit to the data ($\chi^{2}=3.6$ for $13$ degrees of
freedom, d.o.f) with a $99.5$ % goodness-of-fit (GoF). The same $\chi^{2}$
analysis for the meVSL model provides the best-fit value of parameter
$b=0.198$ with the $1$-$\sigma$ error interval being $0.415$. The minimum
$\chi^{2}$ value is $3.4$ with $12$ d.o.f. Its GoF is $99.3$%. These results
are summarized in Table. 2.

Figure 1: Time dilation for two models with data points. The dashed line
depicts the SC prediction of TD. The thick solid line represents the best-fit
value predicted TD of the meVSL model. The shaded region means the
$1$-$\sigma$ confidence interval region.

We also show this in Fig. 1. The vertical lines depict $13$ data points with
their measurement errors. The dashed line represents the SC prediction of the
aging rate $1/(1+z)$. The solid line with shade regions depicts the meVSL
model within the $1$-$\sigma$ confidence interval. This result shows that the
current data cannot distinguish the meVSL from the SC.

Both the SC and the meVSL have almost the same GoF values. It means that both
models provide a similar level of agreement between the observed data and
predictions of models. Thus, both models are almost equally preferred to
current SNe data. Furthermore, the result from the analysis of the meVSL model
cannot be distinguished from that of SC due to its large error. However, there
are about $1500$ SNe LCs available in the redshift range $0\leq z\leq 2.3$
from Pantheon22 data Scolnic et al. (2022); Brout et al. (2022). The
1-$\sigma$ error of $b$ from the Pantheon22 data can be $0.01$ if $\sigma_{i}$
is about $0.1$ and $b\sim 0.1$ by using Fisher matrix analysis (This is the
preliminary result of our ongoing research). Thus, if TD data for Pantheon22
data is available, then we might have a chance to distinguish the meVSL model
from the SC.

##  4 Conclusion

We have driven a formula for the time dilation of a distant object in the
context of the minimally extended varying speed of the light model
$T(z)=T_{0}(1+z)^{1-b/4}$. We obtain the best-fitting values for the b
exponent as $b=0.198\pm 0.415$ at the $1$-$\sigma$ confidence level from $13$
high-redshift SNe Ia time dilation data. This result is worse than one
obtained using the cosmological chronometers Lee (2023b). The best-fitting
values for the $b$ are $b=-0.105\pm 0.178$ and $0.584\pm 0.184$ adopting
($H_{0},\Omega_{m0}$) values from Planck18 Planck Collaboration et al. (2020)
and Pantheon22 Brout et al. (2022), respectively. The result also shows that
the current data are consistent with both the standard model expectation and
that of the meVSL model, and we cannot distinguish between the meVSL model and
the standard one from the time dilation data using SNe. There is also time
dilation data from gamma-ray bursts Singh & Desai (2022). However, the current
unbinned GRB data is too sparse to be directly used for investigating models
with time dilation. Recently, there is a confirmation of the cosmic time
dilation in 190 quasars in the SDSS Stripe 82 region Lewis & Brewer (2023)
through long-term multi-waveband monitoring spanning over two decades using
SDSS, PanSTARRS-1, the Dark Energy Survey, and dedicated follow-up monitoring
with Blanco 4m/DECam Stone et al. (2022). Using Bayesian analysis, they
successfully detect the expected redshift-dependent time dilation in quasar
variability as $n\equiv(1-b/4)=1.28^{+0.28}_{-0.29}$. This result contradicts
previous claims of its absence Hawkins (2010). Previous failure to detect time
dilation in quasar variability may be attributed to limitations such as small
sample sizes and insufficient data sampling and characterization. However, a
significantly larger and diverse sample of quasars across various redshifts
makes it possible to detect time dilation. Thus, we might have a chance to
prefer the specific model in future observations even though the current SNe
Ia data cannot provide the statistically meaningful result to distinguish
models.

## Acknowledgments

SL is supported by Basic Science Research Programs through the National
Research Foundation of Korea (NRF) funded by the Ministry of Science, ICT, and
Future Planning (Grant No. NRF-2017R1A2B4011168 and No.
NRF-2019R1A6A1A10073079).

## Data Availability

The authors confirm that the data supporting the findings of this study are
available within the references and their supplementary materials. No new data
were created or analyzed in this study.

## References

  * Balcerzak & Dąbrowski (2014) Balcerzak A., Dąbrowski M. P., 2014, Physics Letters B, 728, 15
  * Barrow & Magueijo (2000) Barrow J. D., Magueijo J., 2000, ApJ, 532, L87
  * Blondin et al. (2008) Blondin S., et al., 2008, ApJ, 682, 724
  * Brout et al. (2022) Brout D., et al., 2022, ApJ, 938, 110
  * Cuzinatto et al. (2023) Cuzinatto R. R., Holanda R. F. L., Pereira S. H., 2023, MNRAS, 519, 633
  * Foley et al. (2005) Foley R. J., Filippenko A. V., Leonard D. C., Riess A. G., Nugent P., Perlmutter S., 2005, ApJ, 626, L11
  * Goldhaber et al. (2001) Goldhaber G., et al., 2001, ApJ, 558, 359
  * Gupta (2020) Gupta R. P., 2020, MNRAS, 498, 4481
  * Harrison (2000) Harrison E., 2000, Cosmology: The Science of the Universe. Cambridge University Press.
  * Hawkins (2010) Hawkins M. R. S., 2010, MNRAS, 405, 1940
  * Ishak (2019) Ishak M., 2019, Living Reviews in Relativity, 22, 1
  * Islam (2001) Islam N. J., 2001, An Introduction to Mathematical Cosmology. Cambridge University Press.
  * Lee (2021) Lee S., 2021, J. Cosmology Astropart. Phys., 2021, 054
  * Lee (2023a) Lee S., 2023a, Foundations of Physics, 53, 40
  * Lee (2023b) Lee S., 2023b, MNRAS, 522, 3248
  * Leibundgut et al. (1996) Leibundgut B., et al., 1996, ApJ, 466, L21
  * Lewis & Brewer (2023) Lewis G. F., Brewer B. J., 2023, arXiv e-prints, p. arXiv:2306.04053
  * Minkowski (1939) Minkowski R., 1939, ApJ, 89, 156
  * Phillips (1993) Phillips M. M., 1993, ApJ, 413, L105
  * Planck Collaboration et al. (2020) Planck Collaboration et al., 2020, A&A, 641, A6
  * Riess et al. (1997) Riess A. G., et al., 1997, AJ, 114, 722
  * Roos (2015) Roos M., 2015, Introduction to Cosmology. Wiley.
  * Scolnic et al. (2022) Scolnic D., et al., 2022, ApJ, 938, 113
  * Singh & Desai (2022) Singh A., Desai S., 2022, J. Cosmology Astropart. Phys., 2022, 010
  * Stone et al. (2022) Stone Z., et al., 2022, MNRAS, 514, 164
  * Weinberg (1972) Weinberg S., 1972, Gravitation and Cosmology: Principles and Applications of the General Theory of Relativity. John Wiley & Sons, Inc.
  * Wilson (1939) Wilson O. C., 1939, ApJ, 90, 634


# Diffusion as a First Model of Spread of Viral Infection

Paulo H. Acioli<EMAIL_ADDRESS>Department of Physics and Astronomy,
Northeastern Illinois University, Chicago, IL 60625

###### Abstract

The appearance of the coronavirus (COVID-19) in late 2019 has dominated the
news in the last few months as it developed into a pandemic. In many
mathematics and physics classrooms, instructors are using the time series of
the number of cases to show exponential growth of the infection. In this
manuscript we propose a simple diffusion process as the mode of spreading
infections. This model is less sophisticated than other models in the
literature, but it can capture the exponential growth and it can explain it in
terms of mobility (diffusion constant), population density, and probability of
transmission. Students can change the parameters and determine the growth rate
and predict the total number of cases as a function of time. Students are also
given the opportunity to add other factors that are not considered in the
simple diffusion model.

##  I Introduction

The end of 2019 and beginning of 2020 have been dominated by the spread of the
coronavirus (COVID-19). The disease started in the Hubei province in China in
December 2019 and by early January 2020 started to spread. It grew very
rapidly, triggering responses from the Chinese and other governments. On March
11, 2020 the World Health organization declared that COVID-19 was then
characterized as a pandemic.pandemic  The situation triggered different
responses from different governments. In The United States, many colleges took
the initiative to start their own social distance programs, including sending
all students home, extending spring breaks, and ultimately moving all classes
to distance learning.response1 ; response2  Many cities and states followed up
with shelter at home mandates. response3  The situation in many countries
became alarming due to the exponential growth of new cases and deaths. JH ;
WHO

One of the consequences of this pandemic is that instructors at colleges and
universities started to monitor and model the data, using this as a teachable
moment for students and colleagues. The first step in modeling the data is to
plot the number of cases as a function of time and show that it exhibits an
exponential growth. For beginners it is important to introduce them to the
logarithmic scale, where the plot becomes a straight line and students can
extract the exponent and realize that the fit to the US data on March 20, 2020
shows that the number of cases doubled every 2.4 days. Fits to the initial
data for Chicago and New York City show that New York has a higher growth rate
than Chicago FIG. 1, and the total number of cases in Chicago is substantially
smaller. One can speculate that this might be due to difference in population
densities in these cities. After the first 20 days and after mitigation and
social distance measures were imposed, it is clear that the rate of infection
slowed down.

Figure 1: COVID-19 cases in Chicago (CHI,blue)CHI-C19  and New York
CityNYC-C19  (NYC,red). The dashed lines are the best fit to the date after
the initial spike in cases. The $x$-axis is the number of days from the date
that the first case in each city was detected.

In this work we propose to look at some of the factors affecting the spread of
viruses using a simple diffusion model in which each individual in a
population is treated as a Brownian brownian  particle with diffusion constant
$D$. Also added to this model is the incubation period of the virus and a
probability of transmission of the virus if individuals are closer than a
certain distance. This model is to be used as a project in a computational
physics course and verify if it adequately predicts the exponential growth of
the number of cases, as well as if the population density, mobility, and
probability of transmission play roles on the percentage of the population
that will be infected as a function of time. Students, will be asked to modify
a code, analyze the output for different sets of parameters, and write a
critical analysis of its predictive effectiveness.

##  II Computational Project Details

In this section we propose a project to be implemented in computational
physics courses to study the spread of infectious diseases as a simple
diffusion of individuals. The benefits of this project is that its
implementation is simple, but it can lead to a qualitative understanding on
how diseases are spread and it can also allow the student to understand
factors that can affect it. The main ingredients of the model are: individuals
are considered particles that obey a Brownian diffusion process; each
individual will have three possible states, healthy, sick (contagious), and
cured; a healthy individual has a probability of getting infected if its
distance to a sick individual is smaller than a certain threshold; the
incubation and sickness periods are the same; once an individual gets cured,
it cannot be infected or contagious again. VPython (or GlowScript)vpython  is
the language of choice as it allows for a real time visualization of the
infection spread and it will allow for fast simulations of small populations
on a laptop. In the assignment, students will implement the code, analyze the
data generated, critique the initial assumptions, and propose improvements for
more realistic simulations. Below we describe the algorithm starting with the
standard diffusion equation.

The diffusion equation in 2D

$\frac{\partial{f(x,y,t)}}{\partial{t}}=D\left[\frac{\partial^{2}f(x,y,t)}{\partial{x^{2}}}+\frac{\partial^{2}f(x,y,t)}{\partial{y^{2}}}\right]$
(1)

is used to study many phenomenaajp1 ; ajp2 ; ajp3 ; ajp4 ; ajp5 ; ajp6 ; ajp7
; ajp8 ; ajp9 ; ajp10  from diffusion inside the nucleus to population
dynamics to solving the Shrödinger equation. The function $f$ will have
different meanings depending on the application, from temperature to density.
In this work we will use it for the diffusion of individuals, treated as
particles, over closed boundaries subject to contamination of a viral
infection. The normalized solution to Eq. (1) for a single particle is given
by

$f(x,y,t)=f_{0}\frac{1}{\sqrt{4\pi Dt}}e^{-\frac{x^{2}+y^{2}}{4Dt}}.$ (2)

Therefore one can simulate the diffusion of a particle from its previous
position by generating a Gaussian distribution of zero mean and variance
$\sqrt{2Dt}$. For a system on $N$ non-interacting particles with the same
diffusion constant we use Eq. (2) for each particle at each simulation time
step.

The next ingredients in the simulation will be the population density $\rho$,
the number of habitants (particles in the simulation cell) $Npop$, the
diffusion constant $D$, the number of simulation steps $Nstep$, the time step
$dt$, the incubation period ($t_{inc}$), the transmission radius, and the
probability of transmission from an infected to a healthy individual ($prob$).
All of these variables are set at the beginning of the simulation. We preset
that the total simulation time consists of 90 days and that each time step is
0.01 days, therefore each simulation takes 9000 steps. The algorithm is
described below:

  1. 1.

Input $Npop$, $Nstep$, $D$, $dt$, $t_{inc}$, $prob$, $r_{transm}$

  2. 2.

Calculate the size of the square cell as $L=\sqrt{(}Npop/\rho))$

  3. 3.

Initialize the initial population

  4. 4.

Choose a fraction of the initial population to be infected, and set timer for
the sickness ($t_{sick}$)

  5. 5.

Loop over $Nstep$

  6. 6.

Move all individuals according to the Gaussian distribution (Eq. (2)).

  7. 7.

Compute the distance between each healthy and infected individuals

If the distance is less than $r_{transm}$, the healthy individual becomes sick
with probability $prob$

  8. 8.

Subtract the sickness timer by $dt$

  9. 9.

If $t_{sick}<0$ the sick individual gets cured

In order to generate more accurate statistics we suggest the students run
simulations with the same initial parameters multiple times, between 20 and
100, depending on the size of the system and the speed of the student’s
computer. We suggest that the smallest population contains 100 individuals, as
even one sick individual corresponds to an initial infected population of 1%.
Depending on how much time the students have to complete the analysis of the
project, they can use a population of 1000, being aware that each individual
9000 steps simulation can take up to a 30 minutes. They can speed up the
process by using larger time steps, however, they must test if the results of
an individual simulation with the same diffusion constant and different time
steps lead to similar outcomes.

##  III An Example of a Project

In this section we provide an example of a project that students could pursue
and also validate the method by comparing the results of the simulation with
the data presented in FIG. 1. In this project, students will study the
differences of the infection proliferation in New Your City (NYC) and in
Chicago. New York City was chosen because it is experiencing a very rapid
growth on the number of cases and preventive measures such as shelter at home
were taken by the city and state governments at very early stages of the
infection. Chicago was chosen as a local connection to Northeastern Illinois
University (NEIU) students, and is also experiencing an exponential growth in
the number of cases and has also been affected by shelter at home mandates. We
must be mindful that the current model will not be able to study the effects
of the measures that the government are taking, but it might be able to
justify their need.

The population density of New York city is 10,194 people/km2,nycp  while
Chicago’s is 4,665 people/km2.chip  It is natural to offer the students the
hypothesis that if all other variables are the same, the spread of infections
in New York City will be faster than in Chicago and that the number of cases
will be much larger for the same period of time after the first case. With
this hypothesis alone, students should be able to generate a enough data for
the project. They can also discuss if the mobility (diffusion constant) should
be the same in both cases and study the effect of mobility in the spread of
the disease. In the proposed model the period of incubation is the same as the
period of sickness, students can be offered the option to modify this
assumption. The model also assumes that once cured, an individual will attain
immunity, and will not be able to spread the disease. They can discuss
modifications to the model to incorporate relapse. It is clear that with this
very simple model and this very limited two-city project, they can perform a
very thorough study that can give us some qualitative understanding of
infectious disease spread.

In order to validate our model we will try to recreate the same initial
exponential growth as experience by Chicago and New York City. The exponential
fit to the data in FIG. 1 yields exponents of 0.352 and 0.455 for CHI and NYC,
respectively. Although we envision that we can see such growth in a 100
individual population, this growth is unattainable for the same time frame as
the real data, as 100% of the population would be infected in days 13 and 11
if we start with one individual infected in day 0. Deviations from the
exponential behavior will be seen even earlier as there is not enough
susceptible population available to be infected. In order to get similar
exponential growth we start with simulations of populations of 100 individuals
and 1 sick individual chosen at random. The initial diffusion constant was
chosen as $D=100$m2/day with a time step of $dt=0.01$day. Therefore the
variance of the gaussian distribution is $\sqrt{2}$. In this case the number
of simulation steps is 9000 for each individual simulation is 27s in a laptop
with Intel(R) Core(TM) i5-8250U CPU @ 1.60GHz. Therefore 50 simulations will
take about 22.5 minutes of computation. We assume that the radius of
contamination is 2m, and that the probability of contamination is 20% per time
step. The incubation period is taken as 14 days. Below are some of the results
of these simulations.

In FIG. 2a we present the results of the simulation for New York City with the
set of parameters in the previous paragraph. As one can see, under this
assumption 50% of the population is sick on day 11 and 99% of the population
will be infected after 28 days. One can see that after day 18 the number of
sick people starts to decline. However, if the death rate is similar to what
has been observed on the COVID-19 pandemic, about 5.5% of the population of
NYC would perish, and the numbers would be even worse, since no major city in
the world would be able to have hospital beds for 83.8% of its population at
the peak of the infection. In the figure we also observe that the rate of cure
follows the number of infected individuals with a lag time of 14 days, which
is the incubation/sickness period. An exponential fit, shown in FIG. 2b,
yields an exponent of 0.41 which is a bit lower than what was observed in the
real data. Students will need to search for the set of parameters that is able
to fully match the initial growth rate data that they choose to model.

Figure 2: a) Average of 50 90-day simulations of the spread of a virus in a
population of 100 individuals in a square cell with the population density of
NYC ($\rho=0.012$ people/m2). $D=100$m2/day, $prob=0.2$, $dt=0.01$day. b)
Graph of the number of infected individuals (full line) in a logarithmic scale
and the exponential fit for the initial 10 days of the simulation (dashed
line).

In FIG. 3a we report the results for the city of Chicago. Most of the
parameters are the same, with the exception of the density that will be
changed to $\rho=0.0047$ people/m2. For this set of parameters 50 % of the
population will be infected after 26 days reaching a maximum of 87.8% after
day 86. The peak of the number of sick people is 42.3 % and it happens about
32 days after the first 1% of the population is infected. In FIG. 3b we
present a logarithmic scale plot of the total number of infections together
with and exponential fit for the initial stages of the infection. Similarly to
what happened in the simulation for NYC it is clear that the number of
infections deviates from the exponential growth after 13 days. The exponent in
the fit is 0.2 which is also smaller than the value for the data for Chicago.
A closer fit is likely to take place by increasing the probability of
infection.

These two examples do support our initial hypothesis that under the same
conditions, one would expect a slower growth rate in a less dense population.
This leads to the conclusion that keeping all parameters the same, the
percentage of the population infected is correlated to the population density.
In order to determine how mobility affects the rate of infection, one would
change the diffusion constant and repeat the simulations. The expectation is
that the lower the mobility, the lower the total infection rate. In addition,
even with the limitations of a small population simulation we were able to
find growth rates of infections that are similar to what was observed for the
cities of Chicago and New York. Thus validating this approach to simulate the
spread of infectious diseases.

Figure 3: Average of 50 90-day simulations of the spread of a virus in a
population of 100 individuals in a square cell with the population density of
the city of Chicago ($\rho=0.0047$ people/m2). $D=100$m2/day, $prob=0.2$,
$dt=0.01$day.b) Graph of the number of infected individuals (full line) in a
logarithmic scale and the exponential fit for the initial 10 days of the
simulation (dashed line).

In order to show the limitations of small population simulations we run a
simulation for the city of Chicago with a probability of infection of 30% per
time step but increased the number of individuals to 10000. The results are
shown in FIG. 4a. With this set of parameters we obtain the same exponential
growth as observed in FIG. 1 for a period of 15 days. Starting with 0.01 % of
the population infected at the end of 90 days about 40 % of the population
will have been infected. One can see that even a small diffusion constant of
100 m2/day and a probability of transmission of 30 %, the virus will spread to
a large percentage of the population. It is clear that a diffusion model does
a very good job of simulating the spread of the COVID-19 virus and could be
used in a classroom setting, and with the modifications suggested below one
can even use it for more realistic predictions when simulating larger
populations. The main limitation will be access to better computational
resources such as a parallelized code and a computer cluster.

Figure 4: One 90-day simulations of the spread of a virus in a population of
10000 individuals in a square cell with the population density of the city of
Chicago ($\rho=0.0047$ people/m2). $D=100$m2/day, $prob=0.3$, $dt=0.01$day. b)
Graph of the number of infected individuals (full line) in a logarithmic scale
and the exponential fit for the initial 15 days of the simulation (dashed
line).

To finalize, we show snapshots of the simulation cell and its population in
FIG. 5. White spheres represent the susceptible people, red represent sick,
and green are those individuals that recovered. This window is very useful for
a quick analysis of what is happening and can show how the infection is spread
over time. We found it useful to change the radius of the individuals so that
they are visible as the population increases, and the cell size on the screen
remains the same. Because we are displaying the results for 10000 individuals
it is interesting to see that the infection has a clear origin and the
diffusion process spreads the disease outwards. In these simulations we chose
closed boundaries and as a result there will be a limitation on the growth as
there will be less susceptible neighbors to be infected when the disease
reaches a boundary. To avoid this limitation one can use periodic boundary
conditions.

Figure 5: Snapshots of one of the simulation 90-day simulations of the spread
of a virus in a population of 10000 individuals in a square cell with the
population density of CHI ($\rho=0.0047$ people/m2). $D=100$m2/day,
$prob=0.3$, $dt=0.01$day.

The main limitation of this model is that it assumes the same mobility for all
individuals and a constant density in the simulation cell, and no travel
between cities. In addition, the average distance traveled by Brownian
particles is proportional to the square of the simulation time. Thus, there is
a limit on how far they will travel and how many people an individual can
infect. However, depending on the level of the students, instructors can use
the same code and create sub-regions in the simulation cell that would have
different population densities and subsets of the populations can be confined
to these regions. The students could allow individuals to go from one region
to another with a given probability. This would allow the possibility of an
infected individual to move to a region that otherwise would have no
infections. This modification would allow for the spread of the infection
across borders. One could also include quarantine effects, by creating small
cells where a single infected individual will be confined for a period of time
and no other individuals are allowed in. Social distancing can be implemented
by a small repulsive potential attributed to some individuals. These
modifications will bring this simple model closer to more realistic
simulations as those of references ndlib ; infdiffusion . To help the
implementation of these projects, the original code can be downloaded from the
author’s website code .

##  IV Conclusion

This project on itself is simple enough that it can be implemented and
analyzed in a first course on computational physics. However, it is rich
enough that can lead to a large amount of data that can be used to
qualitatively and quantitatively analyze the spread of viral infection. In
addition, it can help students understand the phenomenon of particle diffusion
and Brownian motion. In the example presented, the simulation was limited to
the same 20% chance of infection, the same time step, the size of the
population was fixed, and the length of the disease was the same as the
incubation period. Under these conditions students will have a larger number
of parameters to change and verify if they can reproduce data readily
available in the news. In addition, they will have the opportunity to
determine that the initial growth rate in the number of cases follows an
exponential trend and notice that once it reaches a certain threshold if it
will follow the expected logistic behavior.

In conclusion, we hope to have convinced the reader that a simple diffusion
model can be used to qualitatively explain the spread of disease and even in
some cases quantify it. It is our judgment that it will allow students to work
on a problem that is directly affecting them and that these simulations will
help them offer valuable insight on the issue.

###### Acknowledgements.

I would like to thank Dr. Greg Anderson for sparking the discussion of using
the COVID-19 pandemic as a teachable moment and Dr. Orin Harris for inspiring
the idea of comparing the data for different populations.

## References

  * (1) World Health Organization Director-General’s March 11, 2020 COVID-19 Media Briefing. (https://www.who.int/dg/speeches/detail/who-director-general-s-opening-remarks-at-the-media-briefing-on-covid-19—11-march-2020)
  * (2) https://www.wsj.com/articles/coronavirus-prompts-colleges-to-send-students-home-11583862936.
  * (3) https://www.usatoday.com/videos/news/2020/03/14/coronavirus-colleges-spring-break/5046501002/
  * (4) https://www.usatoday.com/story/news/nation/2020/03/21/coronavirus-lockdown-orders-shelter-place-stay-home-state-list/2891193001/
  * (5) https://coronavirus.jhu.edu/map.html
  * (6) https://www.who.int/emergencies/diseases/novel-coronavirus-2019/situation-reports
  * (7) https://www.chicago.gov/city/en/sites/covid-19/home.html
  * (8) https://www.ny1.com/nyc/all-boroughs/news/health-and-medicine/coronavirus-blog
  * (9) Feynman, R. (1964). ”The Brownian Movement”. The Feynman Lectures of Physics, Volume I. pp. 41–1.
  * (10) https://vpython.org/
  * (11) F. A. Williams, Amer. J. Phys. 26, 467 (1958).
  * (12) D. C. Kelly, Amer. J. Phys. 36, 585 (1968).
  * (13) Karmeshu, and L. S. Kothari, Amer. J. Phys. 40, 1264 (1972).
  * (14) C. Domb and E. L. Offenbacher, Amer. J. Phys. 46, 49 (1978).
  * (15) J. D’Arruda, and E. W. Larsen, Amer. J. Phys. 46, 392 (1978).
  * (16) A. M. Albano, N. B. Abraham, D. E. Chyba, and M. Martelli, Amer. J. Phys. 52, 161 (1984).
  * (17) S. Redner and P. L. Krapivsky, Amer. J. Phys. 67, 1277 (1999).
  * (18) K. Ghosh, K. A. Dill, M. M. Inamdar, E. Seitaridou, and R. Phillips, Amer. J. Phys. 74, 123 (2006).
  * (19) D. R. Spiegel and S. Tuli, Amer. J. Phys. 79, 747 (2011).
  * (20) T. Pang, Amer. J. Phys. 82, 980 (2014).
  * (21) ”NYC Population: Current and Projected Populations”. www1.nyc.gov. Retrieved May 2, 2019.
  * (22) https://www.opendatanetwork.com/entity/1600000US1714000/Chicago_IL/geographic.population.density?year=2018
  * (23) G. Rossetti, L. Milli, S. Rinzivillo, A. Sîrbu, D. Pedreschi, and F. Giannotti Int. J. Data Sci. and Anal. 5, 61(2018).
  * (24) https://towardsdatascience.com/modelling-the-coronavirus-epidemic-spreading-in-a-city-with-python-babd14d82fa2
  * (25) The code can be downloaded from http://physics.neiu.edu/research/viraldiffusion/Viral_Diffusion.py. Instructors and students using the code should cite the present manuscript in any presentation or publication that uses this code as downloaded or modified.