# Epidemics on Hypergraphs: Spectral Thresholds for Extinction

Desmond John Higham School of Mathematics, University of Edinburgh, Edinburgh,
EH9 3FD, UK (d.j.higham@ed.ac.uk). Both authors were supported by Engineering
and Physical Sciences Research Council grant EP/P020720/1.  Henry-Louis de
Kergorlay School of Mathematics, University of Edinburgh, Edinburgh, EH9 3FD,
UK<EMAIL_ADDRESS>
###### Abstract

Epidemic spreading is well understood when a disease propagates around a
contact graph. In a stochastic susceptible-infected-susceptible setting,
spectral conditions characterise whether the disease vanishes. However,
modelling human interactions using a graph is a simplification which only
considers pairwise relationships. This does not fully represent the more
realistic case where people meet in groups. Hyperedges can be used to record
such group interactions, yielding more faithful and flexible models, allowing
for the rate of infection of a node to vary as a nonlinear function of the
number of infectious neighbors. We discuss different types of contagion models
in this hypergraph setting, and derive spectral conditions that characterize
whether the disease vanishes. We study both the exact individual-level
stochastic model and a deterministic mean field ODE approximation. Numerical
simulations are provided to illustrate the analysis. We also interpret our
results and show how the hypergraph model allows us to distinguish between
contributions to infectiousness that (a) are inherent in the nature of the
pathogen and (b) arise from behavioural choices (such as social distancing,
increased hygiene and use of masks). This raises the possibility of more
accurately quantifying the effect of interventions that are designed to
contain the spread of a virus.

##  1 Introduction

Compartmental models for disease propagation have a long and illustrious
history [2, 16], and they remain a fundamental predictive tool [9, 13]. For a
stochastic, individual-level model, it has been suggested recently that
hyperedge information should be incorporated [8, 15, 17]. Hyperedges allow us
to account directly for group interactions of any size, rather than, as in,
for example, [11, 18, 21], treating them as a collection of essentially
independent pairwise encounters.

In this work we contribute to the modeling and analysis of disease spreading
on a hypergraph. We include the case where the number of infected nodes in a
hyperedge contributes nonlinearly to the overall infection rate; this covers
the so-called _collective contagion model_ setting and a new alternative that
we call a _collective suppression model_. The main contributions of our work
are

  * •

a mean field approximation (6)–(7) with a spectral condition for local
asymptotic stability of the zero-infection state (Theorem 6.1) and an
extension to global asymptotic stability (Theorem 6.4) when the nonlinear
infection function is concave,

  * •

for the exact, individual-level model, a spectral condition for exponential
decay of the non-extinction probability in the concave case (Theorem 8.1) and
a spectral bound on the expected time to extinction (Corollary 8.2),

  * •

extensions of these results to more general partitioned hypergraph models,
where distinct infection rates apply to different categories of hyperedge (4),

  * •

results for the non-concave collective contagion model (Theorem 9.1 and
Theorem 9.2),

  * •

a complementary condition that rules out extinction of the disease (Theorem
8.5),

  * •

interpretations of these mathematical results: the spectral thresholds for
disease extinction naturally distinguish between the inherent biological
infectiousness of the disease and behavioural choices of the individuals in
the population, allowing us to account for intervention strategies (Section
10).

The manuscript is organized as follows. In section 2 we introduce the
traditional graph-based susceptible-infected-susceptible (SIS) model and quote
a spectral condition that characterizes control of the disease. We then
discuss the generalization to hyperedges, and motivate the use of infection
rates that do not scale linearly with respect to the number of infected
neighbors. Section 3 formalizes the hypergraph model and shows how it may be
simulated. In section 4 we derive a mean field approximation to characterize
the behaviour of the model, and in section 5 we give computational results to
illustrate its relevance. Section 6 analyses the deterministic mean field
setting and gives a spectral condition for long-term decay of the disease. The
result is local for general infection rates and global (independent of the
initial condition) for the concave case. This spectral condition generalizes a
well-known result concerning disease propagation on a graph. Section 7
provides further computational simulations to illustrate the spectral
threshold. The full stochastic model is then studied in section 8, where we
extend the analysis in [11] to our hypergraph setting. Here we study
extinction of the disease in the case where the nonlinearity in the infection
rate is concave. We also derive conditions for non-extinction. In section 9 we
extend our analysis to the so-called collective contagion model proposed in
[8, 15, 17]. In section 10 we summarize and interpret our results, and discuss
follow-on work.

For a review of recent studies of spreading processes on hypergraphs,
including the dissemination of rumours, opinions and knowledge, we recommend
[3, subsection 7.1.2]. The model that we study fits into the framework of [6].
This work introduced the idea of a nonlinear “infection pressure” from each
hyperedge, and derived a mean field approximation that was compared with
microscale-level simulation results. In [15], the authors studied this type of
model on simplicial complexes of degree up to two (a subclass of the more
general hypergraph setting) and also studied a mean field approximation. These
authors examined the mean field system from a dynamical systems perspective
and analysed issues such as bistability, hysteresis and discontinuous
transitions. Similarly, in [8, 17], a hypergraph version was considered. Our
work differs from these studies in (a) focusing on the derivation of spectral
thresholds for extinction of a disease in both the exact and mean field
settings and (b) seeking to interpret the results from a mathematical
modelling perspective. We mention that it would also be of interest to develop
corresponding thresholds for the mean field models in [6, 8, 15].

##  2 Stochastic SIS models

###  2.1 Stochastic SIS model on a graph

Classical ODE compartmental models are based on the assumption that any pair
of individuals is equally likely to interact—this is the homogeneous mixing
case [2]. If, instead, we have knowledge of all possible pairwise interactions
between individuals, then this information may be incorporated via a contact
graph and used in a stochastic model. Here, each node represents an
individual, and an edge between nodes $\displaystyle i$ and $\displaystyle j$
indicates that individuals $\displaystyle i$ and $\displaystyle j$ interact.
For a population with $\displaystyle n$ individuals, we may let $\displaystyle
A\in\mathbb{R}^{n\times n}$ denote the corresponding symmetric adjacency
matrix, so nodes $\displaystyle i$ and $\displaystyle j$ interact if and only
$\displaystyle A_{ij}=1$. In this setting, a stochastic SIS model uses the
two-state random variable $\displaystyle X_{i}(t)$ to represent the § status
of node $\displaystyle i$ at time $\displaystyle t$, with $\displaystyle
X_{i}(t)=0$ for a susceptible node and $\displaystyle X_{i}(t)=1$ for an
infected node. Each $\displaystyle X_{i}(t)$ then follows a continuous time
Markov process where the infection rate is given by

$\beta\sum_{j=1}^{n}A_{ij}X_{j}(t)$ (1)

and the recovery rate is $\displaystyle\delta$. Here, $\displaystyle\beta>0$
and $\displaystyle\delta>0$ are parameters governing the strength of the two
effects. In this model, we see from (1) that the current chance of infection
increases linearly in proportion to the current number of infected neighbors.

This model was studied in [21, Theorem 1], where it was argued that the
condition

$\lambda(A)\frac{\beta}{\delta}<1$ (2)

guarantees the disease will die out. Here, $\displaystyle\lambda(A)$ denotes
the largest eigenvalue of the symmetric matrix $\displaystyle A$. Further
justification for this result may be found, for example, in [11, 18]. We note
that (2) gives an elegant generalization of the homogeneous mixing case (where
$\displaystyle A$ corresponds to the complete graph).

###  2.2 Why use a hypergraph?

It has been argued [1, 3, 4, 5, 10] that in many network science applications
we lose information by recording only pairwise interactions. For example,
emails can be sent to groups of recipients, scholarly articles may have
multiple coauthors, and many proteins may interact to form a complex. In such
cases, recording the relevant lists of interacting nodes gives a more
informative picture than reducing these down to a collection of edges.

In the setting of an SIS model, we may argue that individuals typically come
together in well-defined groups, for example, in a household, a workplace or a
social setting. Such groups may be handled by the use of hyperedges, leading
to a hypergraph; these concepts are formalized in the next section.

With a classic graph model, as described in section 2.1, the rate of infection
of a node is linearly proportional to the number of infectious neighbors. With
a hypergraph we may consider more intricate contagion mechanisms. For example,
using the terminology of [17], the collective contagion model is used in [8,
15, 17]). Here, infection only starts spreading within a hyperedge after a
certain threshold number of infectious neighbors has been reached. This type
of behaviour is relevant, for example, in an office environment. A small
number of workers may be able to socially distance in way that effectively
eliminates the risk of infection. However, if the number of individuals (size
of the hyperedge) is too large, then the disease may spread.

We mention here that an alternative type of mechanism may also operate, which
we call collective suppression. Imagine that a disease may be contracted
through contact with a surface that was previously touched by an infected
individual. Now suppose that a group of individuals is likely to use the same
physical object, such as a door handle, hand rail, cash machine, or water
cooler. If an infected individual contaminates the object, then further
contamination by other individuals is less relevant. In this case, doubling
the number of common users will increase the risk of infection by a factor
less than two; generally risk grows sublinearly as a function of the size of
the hyperedge.

These arguments motivate us to study the case where the rate of infection of a
node within a hyperdege is dependent on a generic function $\displaystyle f$
of the number of infectious neighbors in a hyperedge; this approach was also
taken in [6]. We will be particularly concerned with the case where
$\displaystyle f$ is concave, since this is tractable for analysis and allows
us to draw conclusions about the collective contagion model.

We note that if $\displaystyle f$ is the identity, then we recover linear
dependence on the number of infectious neighbors and the hypergraph model is
equivalent to a virus spreading on the clique graph of the hypergraph.

##  3 SIS on a Hypergraph

###  3.1 Background

We continue with some standard definitions [7].

######  Definition 1.

A _hypergraph_ is a tuple
$\displaystyle\mathcal{H}\mathrel{\mathop{\ordinarycolon}}=(V,E)$ of _nodes_
$\displaystyle V$ and _hyperedges_ $\displaystyle E$ such that $\displaystyle
E\subset{\mathcal{P}}(V)$. Here, $\displaystyle{\mathcal{P}}(V)$ denotes the
power set of $\displaystyle V$.

We will let $\displaystyle n$ and $\displaystyle m$ denote the number of nodes
and hyperedges, respectively; that is, $\displaystyle|V|=n$ and
$\displaystyle|E|=m$. Loosely, a hypergraph generalizes the concept of a graph
by allowing an “edge” to be a list of more than two nodes.

######  Definition 2.

Consider a hypergraph
$\displaystyle\mathcal{H}\mathrel{\mathop{\ordinarycolon}}=(V,E)$. The
_incidence matrix_ , $\displaystyle{\mathcal{I}}$, is the $\displaystyle
n\times m$ matrix such that $\displaystyle\mathcal{I}_{ih}=1$ if node
$\displaystyle i$ belongs to hyperedge $\displaystyle h$ and
$\displaystyle\mathcal{I}_{ih}=0$ otherwise.

It is also useful to introduce $\displaystyle
W\mathrel{\mathop{\ordinarycolon}}={\mathcal{I}}{\mathcal{I}}^{T}$. This
$\displaystyle n\times n$ matrix has the property that $\displaystyle W_{ij}$
records the number of hyperedges containing both nodes $\displaystyle i$ and
$\displaystyle j$. In particular, if $\displaystyle\mathcal{H}$ is a graph
then $\displaystyle W$ is the affinity matrix of the graph.

###  3.2 General infection model

In our context, the nodes represent individuals and a hyperedge records a
collection of individuals who are known to interact as a group. As in the
graph case introduced in subsection 2.1, we use a state vector $\displaystyle
X(t)$ which follows a continuous time Markov process, where, for each
$\displaystyle 1\leq i\leq n$, $\displaystyle X_{i}(t)=1$ if node
$\displaystyle i$ is infected at time $\displaystyle t$ and $\displaystyle
X_{i}(t)=0$ otherwise. We continue to assume that an infectious node becomes
susceptible with constant recovery rate $\displaystyle\delta>0$. However,
generalizing (1), we now assume that a susceptible node $\displaystyle i$
becomes infectious with rate

$\beta\sum_{h\in
E}\mathcal{I}_{ih}f\big{(}\sum_{j=1}^{n}X_{j}(t)\mathcal{I}_{jh}\big{)},$ (3)

where $\displaystyle\beta>0$ is a constant. In (3), $\displaystyle
f\mathrel{\mathop{\ordinarycolon}}\mathbb{R}_{+}\to\mathbb{R}_{+}$ specifies
the manner in which the contribution to the overall level of infectiousness
from each hyperedge involving node $\displaystyle i$ is assumed to increase in
proportion to the number of infected nodes in that hyperedge. Throughout our
analysis we will always assume that $\displaystyle f(0)=0$ and $\displaystyle
f$ is $\displaystyle C^{1}$ in a neighborhood of $\displaystyle 0$.

If $\displaystyle f$ is the identity, the rate of infection reduces to
$\displaystyle\beta\sum_{j=1}^{n}W_{ij}X_{j}(t)$. This gives a weighted
version of the infection rate of an SIS model on a graph. As discussed in
subsection 2.2, it may be appropriate to choose nonlinear $\displaystyle f$ in
certain circumstances. We note that in [6] the authors have in mind functions
which behave like the identity near the origin and have a horizontal
asymptote. Instances of such functions are $\displaystyle
x\mapsto\mathrm{arctan}(x)$ and $\displaystyle x\mapsto\min\\{x,c\\}$ for some
$\displaystyle c>0$. Relaxing these conditions, we may ask more generally in
such a setting that the function be concave. On the other hand, the authors in
[8, 15] consider a _collective contagion model_ , where infection spreads
within a hyperdge only if a certain threshold of infectious vertices is
reached in that hyperedge. A collective contagion model may be represented via
the function $\displaystyle x\mapsto c_{2}\mathbbm{1}(x\geq c_{1})$ for some
$\displaystyle c_{1},c_{2}>0$, or $\displaystyle x\mapsto\max\\{0,x-c\\}$ for
some $\displaystyle c>0$.

###  3.3 Partitioned hypergraph model

We also introduce a more general case where we partition the hyperedges into
$\displaystyle K$ disjoint categories with each category $\displaystyle 1\leq
k\leq K$ having its own distinct rate of infection in response to the number
infected nodes in a hyperedge, represented by a function $\displaystyle
f_{k}$. For example, the categories may correspond to different types of
housing, workplaces, hospitality venues or sports facilities. We may then
represent the infection rate of node $\displaystyle i$ as

$\beta\sum_{k=1}^{K}\sum_{h\in
E}{\mathcal{I}}^{(k)}_{ih}f_{k}\left(\sum_{j=1}^{n}X_{j}(t){\mathcal{I}}^{(k)}_{jh}\right),$
(4)

where we let $\displaystyle{\mathcal{I}}^{(k)}_{ih}=1$ if $\displaystyle i$
belongs to hyperedge $\displaystyle h$ in the category $\displaystyle k$ and
$\displaystyle{\mathcal{I}}^{(k)}_{ih}=0$ otherwise; so
$\displaystyle{\mathcal{I}}^{(k)}$ is the incidence matrix of the
subhypergraph consisting of only the hyperedges from category $\displaystyle
k$. We will refer to this as a _partitioned hypergraph model_.

In this generalized case, a collective contagion model could be defined by
first organizing the hyperedges into categories depending on their size, so
that category $\displaystyle k$ is the set of hyperedges of size
$\displaystyle k+1$. A collective contagion model may then represented, for
example, via the functions $\displaystyle
f_{1}\mathrel{\mathop{\ordinarycolon}}x\mapsto x$, and $\displaystyle
f_{k}\mathrel{\mathop{\ordinarycolon}}x\mapsto c_{2,k}\mathbbm{1}(x\geq
c_{1,k})$, $\displaystyle k\in\\{2,\dots,K\\}$.

##  4 Mean Field Approximation

A classic approach to studying processes such as (3), where infection rates
are random, is to develop a mean field approximation for the expected process

$\displaystyle\big{(}\mathbb{E}[X_{i}(t)]\big{)}_{t\geq
0}=\big{(}\mathbb{P}(X_{i}(t)=1)\big{)}_{t\geq
0}=\mathrel{\mathop{\ordinarycolon}}(p_{i}(t))_{t\geq 0},$

with deterministic rates. In our case, the rate of recovery
$\displaystyle\delta$ is constant, so can remain unchanged. Let us express the
rate of infection (3) of a node solely in terms of the expected processes
$\displaystyle\\{p_{i}(t)\\}_{i=1}^{n}$. To do this we can substitute the
$\displaystyle X_{j}(t)$ appearing in (3) by their expected values
$\displaystyle p_{j}(t)$. The approximate rate of infection for node
$\displaystyle i$ then becomes

$\beta\sum_{h\in
E}\mathcal{I}_{ih}f\big{(}\sum_{j=1}^{n}p_{j}(t)\mathcal{I}_{jh}\big{)}.$ (5)

We arrive at the deterministic mean field ODE

$\frac{dP(t)}{dt}=g(P(t)),$ (6)

where $\displaystyle
g_{i}\mathrel{\mathop{\ordinarycolon}}\mathbb{R}^{n}\to\mathbb{R}$ is defined
by

$g_{i}(P(t))\mathrel{\mathop{\ordinarycolon}}=\beta\sum_{h\in
E}\mathcal{I}_{ih}f\left(\sum_{j=1}^{n}p_{j}(t)\mathcal{I}_{jh}\right)(1-p_{i}(t))-\delta
p_{i}(t).$ (7)

##  5 Simulations and Comparison between Exact and Mean Field Models

Let us emphasize that the approximate infection rates in (5) differ in general
from the expectation of the random rates in (3). When the function
$\displaystyle f$ is concave, however, Jensen’s reverse inequality indicates
that the rates in (5) are greater than the expectation of the rates in (3).
Hence, in this case the expected quantities $\displaystyle p_{i}(t)$ are
overestimated by (6)–(7). This is fine since we are looking for conditions for
the disease to vanish. If $\displaystyle f$ is not concave (e.g., for a
collective contagion model), these expected quantities are underestimated and
it is not clear a priori whether the exact model is well approximated by the
mean field ODE.

In this section we therefore present results of computational simulations in
order to gain insight into the accuracy of our mean field approximation.

###  5.1 Simulation algorithm

Before presenting numerical results, we summarize our approach for simulating
the individual-level stochastic model, which is based on a standard time
discretization; see, for example, [6]. Using a small fixed time step
$\displaystyle\Delta t$, we advance from time $\displaystyle t$ to
$\displaystyle t+\Delta t$ as follows. First, let $\displaystyle
r\in[0,1]^{n}$ be a random vector of i.i.d. values uniformly sampled from
$\displaystyle[0,1]$. For every node $\displaystyle 1\leq i\leq n$,

  * •

when $\displaystyle X_{i}(t)=0$, set $\displaystyle X_{i}(t+\Delta t)=1$ if

$\displaystyle
r_{i}<1-\exp\left(-\beta\sum_{h}\mathcal{I}_{ih}f(\sum_{j}X_{j}(t)\mathcal{I}_{jh})\Delta
t\right),$

and set $\displaystyle X_{i}(t+\Delta t)=0$ otherwise;

  * •

when $\displaystyle X_{i}(t)=1$, set $\displaystyle X_{i}(t+\Delta t)=0$ if

$\displaystyle r_{i}<1-\exp\left(-\delta\Delta t\right),$

and set $\displaystyle X_{i}(t+\Delta t)=1$ otherwise.

###  5.2 Computational results

In the simulations we chose $\displaystyle n=400$ nodes with fixed recovery
rate $\displaystyle\delta=1$. We look at results for different choices of
infection strength $\displaystyle\beta$ and $\displaystyle i_{0}$, the latter
denoting the (independent) initial probability for each node to be infectious.
We simulated the mean field ODE using Euler’s method with time step
$\displaystyle\Delta t=0.05$. The largest size of a hyperedge was
$\displaystyle 5$ and we distributed the number of hyperedges for the
hypergraph randomly as follows: $\displaystyle 300$ edges, $\displaystyle 200$
hyperedges of size $\displaystyle 3$, $\displaystyle 100$ hyperedges of size
$\displaystyle 4$ and $\displaystyle 50$ hyperedges of size $\displaystyle 5$.
To give a feel for the level of fluctuation, the individual-level paths are
averaged over $\displaystyle 10$ runs, each with the same hypergraph
connectivity and initial state.

Figures 1, 2 and 3 show results for three concave choices of $\displaystyle
f$; respectively,

  * •

$\displaystyle f(x)=\min(x,3)$,

  * •

$\displaystyle f(x)=\log(1+x)$,

  * •

$\displaystyle f(x)=\mathrm{arctan}(x)$.

For Figure 4 we used a collective contagion model on a partitioned hypergraph.
Assigning each hyperedge to a category in $\displaystyle\\{1,2,3,4\\}$, where
category $\displaystyle k$ contains the hyperedges of size $\displaystyle
k+1$, we chose the following associated functions to determine the infection
rates: $\displaystyle f_{1}(x)\mathrel{\mathop{\ordinarycolon}}=x$, and for
$\displaystyle k\in\\{2,3,4\\}$, $\displaystyle
f_{k}(x)\mathrel{\mathop{\ordinarycolon}}=(k-1)\mathbbm{1}(x\geq k-1)$.

The four figures show the proportion of infectious individuals as a function
of time. In these simulations, and others not reported here, we observe that
the initial value $\displaystyle i_{0}$ does not affect the asymptotic
behavior of the process: the process vanishes or converges to a non-zero
equilibrium depending on the value of $\displaystyle\beta$ but regardless of
the value of $\displaystyle i_{0}$. In Figures 1, 2 and 3, where
$\displaystyle f$ is concave, we know that the mean field model gives an upper
bound on the expected proportion of infected individuals in the microscale
model. We also see that the mean field model provides a reasonably sharp
approximation. Moreover, we see a similar level of sharpness in Figure 4 for
the collective contagion model, where $\displaystyle f$ is not concave.

Figure 1: Here, $\displaystyle f(x)=\min(x,3)$. Red dashed line: mean field
approximation from (6)–(7). Blue solid line: proportion of infected
individuals, $\displaystyle\sum_{i}X_{i}(t)/n$, from the individual-level
stochastic model (3), averaged over $\displaystyle 10$ runs.  Figure 2: Here,
$\displaystyle f(x)=\log(1+x)$. Red dashed line: mean field approximation from
(6)–(7). Blue solid line: proportion of infected individuals,
$\displaystyle\sum_{i}X_{i}(t)/n$, from the individual-level stochastic model
(3), averaged over $\displaystyle 10$ runs.  Figure 3: Here, $\displaystyle
f(x)=\mathrm{arctan}(x)$. Red dashed line: mean field approximation from
(6)–(7). Blue solid line: proportion of infected individuals,
$\displaystyle\sum_{i}X_{i}(t)/n$, from the individual-level stochastic model
(3), averaged over $\displaystyle 10$ runs.  Figure 4: Collective contagion
model on a partitioned hypergraph. Red dashed line: mean field approximation
from (6) with (9). Blue solid line: proportion of infected individuals,
$\displaystyle\sum_{i}X_{i}(t)/n$, from the individual-level stochastic model
(3), averaged over $\displaystyle 10$ runs.

A key advantage of the mean field approximation is that it gives rise to a
deterministic autonomous dynamical system for which there exists a rich theory
to study the asymptotic stability of equilibrium points. This motivates the
analysis in the next section.

##  6 Stability Analysis

We provide below spectral conditions which imply that the infection-free
solution $\displaystyle 0\in\mathbb{R}^{n}$ is a locally or globally
asymptotically stable equilibrium of (6)–(7). We will find that local
asymptotic stability can be shown with no structural assumptions on
$\displaystyle f$. We will also find that global asymptotic stability follows
_under the same conditions_ when $\displaystyle f$ is concave. Our conclusions
fit into a framework that generalizes the graph case (2): the spectral
threshold takes the form

$\displaystyle\lambda(W)\frac{c_{f}\beta}{\delta}<1$

for some constant $\displaystyle c_{f}>0$ depending only on the choice of
$\displaystyle f$.

Throughout this work, to be concrete we let $\displaystyle\|\cdot\|$ denote
the Euclidean norm.

###  6.1 Local asymptotic stability

######  Theorem 6.1.

If

$\lambda(W)\frac{f^{\prime}(0)\beta}{\delta}<1$ (8)

then $\displaystyle 0\in\mathbb{R}^{n}$ is a locally asymptotically stable
equilibrium of (6)–(7); that is, there exists a positive $\displaystyle\gamma$
such that
$\displaystyle||P(0)||<\gamma\Rightarrow\lim_{t\to\infty}||P(t)||=0$.

###### Proof.

We see that $\displaystyle g(0)=0$, so $\displaystyle 0\in\mathbb{R}^{n}$ is
an equilibrium for (6). It remains to show that this solution is locally
asymptotically stable. Appealing to a standard linearization result [20], it
suffices to show that every eigenvalue of the Jacobian matrix
$\displaystyle\nabla g(0)$ has a negative real part. We compute

$\displaystyle\displaystyle\frac{\partial g_{i}}{\partial p_{j_{0}}}$
$\displaystyle\displaystyle=\begin{cases}\beta\sum_{h}\mathcal{I}_{ih}\mathcal{I}_{j_{0}h}f^{\prime}(\sum_{j}p_{j}\mathcal{I}_{jh})(1-p_{i}),&\text{
}j_{0}\neq i,\\\
\beta\sum_{h}\mathcal{I}_{ih}\mathcal{I}_{j_{0}h}f^{\prime}(\sum_{j}p_{j}\mathcal{I}_{jh})(1-p_{i})-\beta\sum_{h}\mathcal{I}_{ih}f(\sum_{j}p_{j}\mathcal{I}_{jh})-\delta,&\text{
}j_{0}=i.\end{cases}$

We see that $\displaystyle\nabla g(0)=\beta f^{\prime}(0)W-\delta I$. This
matrix is symmetric and therefore has real eigenvalues. Hence, it suffices
that the largest eigenvalue of $\displaystyle\beta f^{\prime}(0)W$ does not
exceed $\displaystyle\delta$, and the result follows. ∎

Theorem 6.1 extends to the partitioned model in (4). In this case
$\displaystyle g_{i}(P(t))$ in the mean field ODE (6) is defined as

$g_{i}(P(t))\mathrel{\mathop{\ordinarycolon}}=\beta\sum_{k=1}^{K}\sum_{h\in
E}\mathcal{I}_{ih}^{(k)}f_{k}\big{(}\sum_{j=1}^{n}p_{j}(t)\mathcal{I}_{jh}^{(k)}\big{)}(1-p_{i}(t))-\delta
p_{i}(t),$ (9)

and we let $\displaystyle
W^{(k)}\mathrel{\mathop{\ordinarycolon}}={\mathcal{I}}^{(k)}{{\mathcal{I}}^{(k)}}^{T}$.

######  Theorem 6.2.

If

$\lambda\left(\sum_{k=1}^{K}f_{k}^{\prime}(0)W^{(k)}\right)\frac{\beta}{\delta}<1$
(10)

then $\displaystyle 0\in\mathbb{R}^{n}$ is a locally asymptotically stable
equilibrium of (6), with $\displaystyle g_{i}$ defined in (9).

###### Proof.

The proof of Theorem 6.1 extends straightforwardly. We compute

$\displaystyle\displaystyle\frac{\partial g_{i}}{\partial p_{j_{0}}}$
$\displaystyle\displaystyle=\begin{cases}\beta\sum_{k}\sum_{h}\mathcal{I}^{(k)}_{ih}\mathcal{I}^{(k)}_{j_{0}h}f_{k}^{\prime}(\sum_{j}p_{j}\mathcal{I}^{(k)}_{jh})(1-p_{i}),&\text{
}j_{0}\neq i,\\\
\beta\sum_{k}\sum_{h}\mathcal{I}^{(k)}_{ih}\mathcal{I}^{(k)}_{j_{0}h}f_{k}^{\prime}(\sum_{j}p_{j}\mathcal{I}^{(k)}_{jh})(1-p_{i})-\beta\sum_{k}\sum_{h}\mathcal{I}^{(k)}_{ih}f_{k}(\sum_{j}p_{j}\mathcal{I}^{(k)}_{jh})-\delta,&\text{
}j_{0}=i,\end{cases}$ $\displaystyle\displaystyle\Rightarrow\ $
$\displaystyle\displaystyle\frac{\partial g_{i}}{\partial
p_{j_{0}}}|_{P=0}=\begin{cases}\beta\sum_{k}W^{(k)}_{ij}f_{k}^{\prime}(0),&\text{
}j_{0}\neq i,\\\ \beta\sum_{k}W^{(k)}_{ij}f_{k}^{\prime}(0)-\delta,&\text{
}j_{0}=i,\end{cases}$

and note that

$\displaystyle\lambda(\beta\sum_{k}f^{\prime}_{k}(0)W^{(k)}-\delta
I)<0\quad\Leftrightarrow\quad\lambda(\sum_{k}f^{\prime}_{k}(0)W^{(k)})<\frac{\delta}{\beta}.$

∎

###  6.2 Global asymptotic stability for the concave infection model

We now show that when $\displaystyle f$ is concave the condition in Theorem
6.1 ensures global stability of the zero equilibrium, and hence guarantees
that the disease dies out according to the mean field approximation.

######  Definition 3.

Given a matrix $\displaystyle A$, define its symmetric version to be

$\displaystyle A^{(S)}\mathrel{\mathop{\ordinarycolon}}=(A+A^{T})/2.$

######  Lemma 6.3.

Suppose that $\displaystyle A$ and $\displaystyle B$ are $\displaystyle
n\times n$ real matrices, and suppose that there exists a diagonal matrix
$\displaystyle\Lambda$ such that for all $\displaystyle
i\in\\{1,2,\dots,n\\},\ \Lambda_{ii}\geq 0$, and

$\displaystyle A=B-\Lambda.$

Then the largest eigenvalues of $\displaystyle A$ and $\displaystyle B$
satisfy $\displaystyle\lambda(A)\leq\lambda(B)$, and the largest eigenvalues
of $\displaystyle A^{(S)}$ and $\displaystyle B^{(S)}$ also satisfy
$\displaystyle\lambda(A^{(S)})\leq\lambda(B^{(S)})$.

###### Proof.

Let $\displaystyle x$ be a unit eigenvector associated with
$\displaystyle\lambda(A^{(S)})$. We have

$\displaystyle\displaystyle 2\lambda(A^{(S)})$
$\displaystyle\displaystyle=x^{T}Ax+x^{T}A^{T}x$
$\displaystyle\displaystyle=x^{T}Bx+x^{T}B^{T}x-2x^{T}\Lambda x$
$\displaystyle\displaystyle=\sum_{i,j=1}^{n}(b_{ij}+b_{ji})x_{i}x_{j}-2\sum_{i=1}^{n}\Lambda_{ii}x_{i}^{2}$
$\displaystyle\displaystyle\leq x^{T}(B+B^{T})x$
$\displaystyle\displaystyle\leq\max\\{x^{T}(B+B^{T})x\ |\
x^{T}x=1\\}=2\lambda(B^{(S)}).$

The inequality $\displaystyle\lambda(A)\leq\lambda(B)$ may be shown similarly.
∎

######  Theorem 6.4.

Suppose $\displaystyle f$ is concave. If

$\displaystyle\lambda(W)\frac{f^{\prime}(0)\beta}{\delta}<1,$

then $\displaystyle 0\in\mathbb{R}^{n}$ is a globally asymptotically stable
equilibrium of (6)–(7); so $\displaystyle\lim_{t\to\infty}||P(t)||=0$ for any
valid initial condition (that is, with $\displaystyle 0\leq p(0)_{i}\leq 1$).

###### Proof.

From the global asymptotic stability result in [14, Lemma $\displaystyle
1^{\prime}$ ] it is sufficient to show that all eigenvalues of the symmetric
matrix $\displaystyle(\nabla g(P))^{(S)}$ are strictly less than
$\displaystyle 0$, for all $\displaystyle P\neq 0$. We have

$\displaystyle\displaystyle\frac{\partial g_{i}}{\partial p_{j_{0}}}$
$\displaystyle\displaystyle=\begin{cases}\beta\sum_{h}\mathcal{I}_{ih}\mathcal{I}_{j_{0}h}f^{\prime}(\sum_{j}p_{j}\mathcal{I}_{jh})(1-p_{i}),&\
j_{0}\neq i,\\\
\beta\sum_{h}\mathcal{I}_{ih}\mathcal{I}_{j_{0}h}f^{\prime}(\sum_{j}p_{j}\mathcal{I}_{jh})(1-p_{i})-\beta\sum_{h}\mathcal{I}_{ih}f(\sum_{j}p_{j}\mathcal{I}_{jh})-\delta,&\
j_{0}=i.\\\ \end{cases}$

Letting $\displaystyle B$ denote the $\displaystyle n\times n$ matrix given by

$\displaystyle\displaystyle
B_{ij_{0}}=\begin{cases}\beta\sum_{h}\mathcal{I}_{ih}\mathcal{I}_{j_{0}h}f^{\prime}(\sum_{j}p_{j}\mathcal{I}_{jh})(1-p_{i}),&\
j_{0}\neq i,\\\
\beta\sum_{h}\mathcal{I}_{ih}\mathcal{I}_{j_{0}h}f^{\prime}(\sum_{j}p_{j}\mathcal{I}_{jh})(1-p_{i})-\delta,&\
j_{0}=i,\\\ \end{cases}$

we have $\displaystyle\nabla g(P)=B-\Lambda$, where $\displaystyle\Lambda$ is
the $\displaystyle n\times n$ diagonal matrix where for all $\displaystyle
i\in\\{1,2,\dots\\}$,
$\displaystyle\Lambda_{ii}\mathrel{\mathop{\ordinarycolon}}=\beta\sum_{h}\mathcal{I}_{ih}f(\sum_{j}p_{j}\mathcal{I}_{jh})\geq
0$. On the one hand Lemma 6.3 now yields $\displaystyle\lambda((\nabla
g(P))^{(S)})\leq\lambda(B^{(S)})$; on the other hand, note that $\displaystyle
B+\delta I\leq\nabla g(0)+\delta I$, where we interpret the inequality in a
componentwise sense, and where we use $\displaystyle
f^{\prime}(\sum_{j}p_{j}\mathcal{I}_{jh})\leq f^{\prime}(0)$, since
$\displaystyle f$ is concave. Hence $\displaystyle B^{(S)}+\delta I\leq\nabla
g(0)+\delta I$, and since $\displaystyle B^{(S)}+\delta I$ has only non-
negative entries, appealing to the Perron-Frobenius theorem, we have
$\displaystyle\lambda(B^{(S)})\leq\lambda(\nabla g(0))$.

Combining these inequalities and using the spectral condition in the statement
of the theorem, we deduce that

$\displaystyle\lambda\left((\nabla g(P))^{(S)}\right)\leq\lambda(\nabla
g(0))=\lambda\left(\beta f^{\prime}(0)W-\delta I\right)<0,$

as required. ∎

A straightforward adaptation of the proof of Theorem 6.4 yields the following
global asymptotic stability result for the more general partitioned model.

######  Theorem 6.5.

Suppose all $\displaystyle f_{k}$ are concave. If

$\displaystyle\lambda(\sum_{k=1}^{K}f_{k}^{\prime}(0)W^{(k)})\frac{\beta}{\delta}<1,$

then $\displaystyle 0\in\mathbb{R}^{n}$ is a globally asymptotically stable
equilibrium of (6), with $\displaystyle g_{i}$ defined in (9).

##  7 Simulations to Test the Spectral Condition

We now show the results of experiments that test the sharpness of our spectral
vanishing condition. Here, we used the concave functions $\displaystyle
f(x)=2\log(1+x)$ (on the left of Figure 5 and in Figure 6) and $\displaystyle
f(x)=\mathrm{arctan}(x)$ (on the left in Figure 5 and in Figure 7) to
construct partitioned models with $\displaystyle f_{1}(x)=x$ and
$\displaystyle f_{k}(x)=f(x)$ for all $\displaystyle k\geq 2$. We fixed a
hypergraph with $\displaystyle n=400$ nodes, $\displaystyle 400$ edges,
$\displaystyle 200$ hyperedges of size $\displaystyle 3$, $\displaystyle 100$
hyperedges of size $\displaystyle 4$ and $\displaystyle 50$ hyperedges of size
$\displaystyle 5$. At time zero, each node was infected with independent
probability $\displaystyle i_{0}=0.5$ and we used a recovery rate of
$\displaystyle\delta=1$. In addition to the mean field ODE, we also simulated
the microscale model, averaged over $\displaystyle 5$ runs, using the
discretization scheme described in Section 5.2, with $\displaystyle\Delta
t=0.1$.

In Figure 5 the asterisks (red) show the corresponding proportion of infected
individuals according to the mean field model,
$\displaystyle\sum_{i}p_{i}(t)/n$, at time $\displaystyle t=150$, for a range
of different $\displaystyle\beta$ between $\displaystyle 0$ and $\displaystyle
0.2$: $\displaystyle\beta\in\\{(0.01)k\ |\ k\in\\{0,\dots,20\\}\\}$. The
crosses (blue) show the corresponding proportion of infected individuals from
the microscale model, $\displaystyle\sum_{i}X_{i}(t)/n$, at time
$\displaystyle t=150$. The vertical green line represents the critical value
$\displaystyle\beta_{c}\mathrel{\mathop{\ordinarycolon}}=\delta/\lambda(\sum_{k=1}^{K}f_{k}^{\prime}(0)W^{(k)})$
(we have $\displaystyle\beta_{c}\approx 0.0265$ on the left and
$\displaystyle\beta_{c}\approx 0.0490$ on the right).

For the mean field model, we know from Theorem 6.5 that
$\displaystyle\beta<\beta_{c}$ guarantees global stability of the zero-
infection state. We see that $\displaystyle\beta_{c}$ also lies close to the
threshold beyond which extinction of the disease is lost in the mean field
model. For the individual-level stochastic model, Theorem 8.3 below shows that
$\displaystyle\beta<\beta_{c}$ is also sufficient for eventual extinction of
the disease. This is consistent with the results in Figure 5.

Figure 5: Left: infection function based on $\displaystyle 2\log(1+x)$. Right:
infection function based on $\displaystyle\mathrm{arctan}(x)$. For different
choices of infection strength $\displaystyle\beta$ (horizontal axis), we show
the proportion of infection individuals at time $\displaystyle t=150$
(vertical axis) for the mean field approximation (6) with (9) in red asterisks
and for the individual-level stochastic model (3) in blue crosses. The
spectral bound arising from our analysis is show as a green vertical line.

In the left of Figures 6 and 7 we show individual trajectories of the
proportion of infected individuals, $\displaystyle\sum_{i}p_{i}(t)/n$,
according to the mean field model, for a range of $\displaystyle\beta$ values.
For the same range of $\displaystyle\beta$ values, the plots on the right of
these figures show the corresponding proportion of infected individuals from
the microscale model, $\displaystyle\sum_{i}X_{i}(t)/n$. The curves are
colored in red if the spectral vanishing condition
$\displaystyle\beta<\beta_{c}$ is satisfied. We see qualitative agreement
between the mean field and individual-level models, and extinction for the
$\displaystyle\beta$ values below the spectral threshold.

Figure 6: Results with the arctan infection function. Left: proportion of
infected individuals using the mean field model. Right: proportion of infected
individuals using the individual-level model. From bottom to top, the
$\displaystyle\beta$ values used were From bottom to top, the
$\displaystyle\beta$ values used were $\displaystyle\beta\in\\{(0.02)k\ |\
k\in\\{0,\dots,10\\}\\}$. Cases where $\displaystyle\beta$ is below the
spectral bound are colored in red.  Figure 7: Results with the arctan
infection function. Left: proportion of infected individuals using the mean
field model. Right: proportion of infected individuals using the individual-
level model. From bottom to top, the $\displaystyle\beta$ values used were
$\displaystyle\beta\in\\{(0.02)k\ |\ k\in\\{0,\dots,10\\}\\}$. Cases where
$\displaystyle\beta$ is below the spectral bound are colored in red.

Having derived and tested spectral conditions that concern extinction of the
disease at the mean field approximation level, in the next section we study
the microscale model directly.

##  8 Exact Model

To proceed, we recall our assumption that at time zero each node has the same,
independent, probability, $\displaystyle i_{0}$, of being infectious; so
$\displaystyle\mathbb{P}(X_{j}(0)=1)=i_{0}$ for all $\displaystyle 1\leq j\leq
n$. This implies that $\displaystyle n\,i_{0}$ is the expected number of
infectious individuals at time zero.

We are interested in the stochastic process $\displaystyle\sum_{i}X_{i}(t)$,
which records the number of infected individuals. Our analysis generalizes
arguments in [11], which considered a stochastic SIS model on a graph with
$\displaystyle f$ as the identity map.

###  8.1 Extinction

Our first result shows that the spectral condition arising from the mean field
analysis in Theorems 6.1 and 6.4 is also relevant to the probability of
extinction in the individual-level model.

######  Theorem 8.1.

Suppose $\displaystyle f$ is concave in the hypergraph infection model (3).
Then

$\displaystyle\mathbb{P}\left(\sum_{i}X_{i}(t)>0\right)\leq
n\,i_{0}\,\exp\left((\beta f^{\prime}(0)\lambda(W)-\delta)t\right).$

Hence, if $\displaystyle\lambda(W)f^{\prime}(0)\beta/\delta<1$ then the
disease vanishes at an exponential rate.

###### Proof.

Consider the continuous time Markov process $\displaystyle\\{(Y_{i}(t))_{t\geq
0}\\}_{i=1}^{n}$ taking values in $\displaystyle\mathbb{N}^{n}$, with
transition of states defined for every $\displaystyle 1\leq i\leq n$ and
$\displaystyle\ t\geq 0$ by

$\displaystyle\begin{cases}k&\rightarrow\ k+1\text{, with rate }\beta
f^{\prime}(0)\sum_{j}W_{ij}Y_{j}(t),\\\ k&\rightarrow\ k-1\text{, with rate
}\delta.\end{cases}$

This new process is introduced here purely for the purpose of analysis.
However, it may be interpreted as a disease model where the state of each
individual is represented by a non-negative integer that indicates severity of
infection. Here, exposure to highly infected individuals raises the chance of
an increase in infection severity.

Suppose also that $\displaystyle X_{i}(0)=Y_{i}(0)$ for all $\displaystyle
1\leq i\leq n$. Since $\displaystyle f$ is concave,

$\displaystyle\beta\sum_{h}\mathcal{I}_{ih}f(\sum_{j}X_{j}(t)\mathcal{I}_{jh})\leq\beta\sum_{h}\mathcal{I}_{ih}f^{\prime}(0)\sum_{j}X_{j}(t)\mathcal{I}_{jh}=\beta
f^{\prime}(0)\sum_{ij}W_{ij}X_{j}(t),$

from which we see that $\displaystyle Y_{i}$ stochastically dominates
$\displaystyle X_{i}$. Hence

$\displaystyle\mathbb{P}\left(\sum_{i}X_{i}(t)>0\right)\leq\mathbb{P}\left(\sum_{i}Y_{i}(t)>0\right)\leq\sum_{i}q_{i}(t),$

where $\displaystyle
q_{i}(t)\mathrel{\mathop{\ordinarycolon}}=\mathbb{E}[Y_{i}(t)]$. In terms of
the Chapman–Kolmogorov Equation, or Chemical Master Equation, [12], we have

$\displaystyle\frac{dq_{i}(t)}{dt}=\beta
f^{\prime}(0)\sum_{j}W_{ij}q_{j}(t)-\delta q_{i}(t).$

Letting $\displaystyle Q(t)=[q_{1}(t),q_{2}(t),\ldots,q_{n}(t)]^{T}$, this
linear ODE system solves to give

$\displaystyle Q(t)=\exp\left(t(\beta f^{\prime}(0)W-\delta I)\right)Q(0).$

The matrix $\displaystyle\exp\left(t(\beta f^{\prime}(0)W-\delta I)\right)$ is
symmetric and has spectral radius $\displaystyle\exp\left((\beta
f^{\prime}(0)\lambda(W)-\delta I)t\right)$. Hence, in Euclidean norm,

$\displaystyle||Q(t)||\leq\exp((\beta f^{\prime}(0)\lambda(W)-I)t)||Q(0)||.$

Since $\displaystyle||Q(0)||=\sqrt{n}\,i_{0}$ and, by Cauchy–Schwarz,

$\displaystyle\sum_{i}q_{i}(t)\leq\sqrt{n}||Q(t)||,$

the proof is complete. ∎

We deduce, analogously to [11], the following corollary.

######  Corollary 8.2.

Suppose $\displaystyle f$ is concave in the hypergraph infection model (3).
Let $\displaystyle\tau$ denote the time of extinction of the disease and
suppose $\displaystyle\lambda(W)f^{\prime}(0)\beta<\delta$, then

$\displaystyle\mathbb{E}[\tau]\leq\frac{\log
n+1}{\delta-f^{\prime}(0)\beta\lambda(W)}.$

###### Proof.

Using Theorem 8.1,

$\displaystyle\displaystyle\mathbb{E}[\tau]$
$\displaystyle\displaystyle=\int_{0}^{\infty}\mathbb{P}(\tau>t)dt$
$\displaystyle\displaystyle=\int_{0}^{\infty}\mathbb{P}(\sum_{i}X_{i}(t)>0)dt$
$\displaystyle\displaystyle\leq\frac{\log
n}{\delta-f^{\prime}(0)\beta\lambda(W)}+\int_{(\log
n)/(\delta-f^{\prime}(0)\beta\lambda(W))}^{\infty}n\exp((\beta
f^{\prime}(0)\lambda(W)-\delta)t)dt$ $\displaystyle\displaystyle\leq\frac{\log
n+1}{\delta-f^{\prime}(0)\beta\lambda(W)}.$

∎

Likewise the partitioned case yields the following result.

######  Theorem 8.3.

Suppose every $\displaystyle f_{k}$ is concave in the partitioned hypergraph
model with infection rate (4). Then

$\displaystyle\mathbb{P}\left(\sum_{i}X_{i}(t)>0\right)\leq
n\,i_{0}\,\exp\left(\beta\lambda(\sum_{k=1}^{K}f^{\prime}_{k}(0)W^{(k)})-\delta\right).$

Hence, if
$\displaystyle\lambda(\sum_{k=1}^{K}f^{\prime}_{k}(0)W^{(k)})\beta/\delta<1$
then the disease vanishes at an exponential rate.

We also have the following analogue of Corollary 7.2 on the expected time to
extinction for the partitioned case.

######  Corollary 8.4.

Suppose every $\displaystyle f_{k}$ is concave in the partitioned hypergraph
model with infection rate (4). Let $\displaystyle\tau$ denote the time of
extinction of the disease and suppose
$\displaystyle\lambda(\sum_{k=1}^{K}f^{\prime}_{k}(0)W^{(k)})\beta/\delta<1$,
then

$\displaystyle\mathbb{E}[\tau]\leq\frac{\log
n+1}{\delta-\beta\lambda(\sum_{k=1}^{K}f^{\prime}_{k}(0)W^{(k)})}.$

###  8.2 Conditions that preclude extinction

So far, we have focused on deriving thresholds that imply extinction. In this
subsection, following ideas from [11], we derive a condition under which the
disease will persist.

Note that our analysis does not require the graph associated with
$\displaystyle W$ to be connected. The disconnected setting is relevant, for
example, when interventions have been imposed in order to limit interactions.
We let $\displaystyle\Delta\mathrel{\mathop{\ordinarycolon}}=D-W$ denote the
Laplacian, and let $\displaystyle\lambda_{c}(\Delta)>0$ denote the smallest
non-zero eigenvalue of $\displaystyle\Delta$. We also let $\displaystyle
e_{\max}$ denote the size of the largest hyperedge.

######  Definition 4.

Given a hypergraph $\displaystyle\mathcal{H}$, a function $\displaystyle
f\mathrel{\mathop{\ordinarycolon}}\mathbb{R}_{+}\to\mathbb{R}_{+}$ and a
subset of the nodes $\displaystyle S\subset V$, let

$\displaystyle E(S,f)\mathrel{\mathop{\ordinarycolon}}=\sum_{i\in S}\sum_{h\in
E}\mathcal{I}_{ih}f(\sum_{j\in S^{c}}\mathcal{I}_{jh}),$

where $\displaystyle S^{c}\mathrel{\mathop{\ordinarycolon}}=V\backslash S$ is
the complement of $\displaystyle S$. Also define for integer $\displaystyle
1\leq m\leq\lfloor n/2\rfloor$

$\displaystyle\eta(\mathcal{H},m,f)\mathrel{\mathop{\ordinarycolon}}=\inf\left\\{\frac{E(S,f)}{|S|}\
|\ \ 1\leq|S|\leq m\right\\},$

and let
$\displaystyle\eta(\mathcal{H},m)\mathrel{\mathop{\ordinarycolon}}=\eta(\mathcal{H},m,Id)$.

Notice that when $\displaystyle S$ consists of those nodes for which
$\displaystyle X_{i}(t)=0$, we can write the infection transition rate of
$\displaystyle\sum_{i=1}^{n}X_{i}(t)$ as $\displaystyle\beta E(S,f)$. More
generally, $\displaystyle\beta E(S,f)$ may be regarded as the rate at which
nodes in the set $\displaystyle S$ may be infected by nodes in the remainder
of the network. When $\displaystyle f=Id$ and $\displaystyle m=\lfloor
n/2\rfloor$, $\displaystyle\eta(\mathcal{H},m)$ is the Cheeger constant, or
isoperimetric number, associated with the weighted graph induced by
$\displaystyle W=\mathcal{I}\mathcal{I}^{T}$. We may also regard
$\displaystyle\eta(\mathcal{H},m,f)$ as the smallest average infection rate
over all subsets consisting of no more than half of the network.

The next theorem gives a probabilistic lower bound on the time to extinction.

######  Theorem 8.5.

Recall that $\displaystyle\tau$ denotes the hitting time of the state
$\displaystyle 0$ for the process $\displaystyle(\sum_{j}X_{j}(t))_{t\geq 0}$
in the hypergraph model (3). If $\displaystyle f$ is concave and
$\displaystyle\lambda_{c}(\Delta)>\left(2\frac{e_{\max}-1}{f(e_{\max}-1)}\right)\frac{\delta}{\beta}$,
then

$\displaystyle\mathbb{P}\left(\tau>\frac{\lfloor
r^{-m+1}\rfloor}{2m}\right)\geq\frac{1-r}{e}(1+O(r^{m})),$

where $\displaystyle
r\mathrel{\mathop{\ordinarycolon}}=\frac{(e_{\max}-1)\delta}{f(e_{\max}-1)\,\beta\,\eta(\mathcal{H},m)}<1$
and $\displaystyle
m\mathrel{\mathop{\ordinarycolon}}=\lfloor\frac{n}{2}\rfloor$.

From standard Cheeger inequalities [19], we know that the Cheeger constant of
the graph induced by $\displaystyle W$ satisfies $\displaystyle
2\eta(\mathcal{H},m)\geq\lambda_{c}(\Delta)$, hence using the assumptions on
$\displaystyle\lambda_{c}(\Delta)$ in the theorem, we see that $\displaystyle
r<1$ indeed.

In order to prove this result, we introduce the following lemma.

######  Lemma 8.6.

If $\displaystyle f$ is concave and non-decreasing, then

$\displaystyle\frac{f(e_{\max}-1)}{e_{\max}-1}\eta(\mathcal{H},m)\leq\eta(\mathcal{H},m,f).$

###### Proof.

The proof is immediate once we see that by concavity of $\displaystyle f$, for
all $\displaystyle x\in\\{0,\dots,e_{\max}-1\\}$

$\displaystyle\frac{f(e_{\max}-1)}{e_{\max}-1}x\leq f(x).$

∎

Now, to prove Theorem 8.5 consider the Markov process
$\displaystyle(Z(t))_{t\geq 0}$ valued in $\displaystyle\\{0,\dots,m\\}$, with
transition of states given by

$\displaystyle\begin{cases}k\ \rightarrow\ k+1\text{ with transition rate
}k\beta\frac{f(e_{\max}-1)}{e_{\max}-1}\eta(\mathcal{H},m),\\\ k\ \rightarrow\
k-1\text{ with transition rate }k\delta.\end{cases}$

This Markov process is stochastically dominated by
$\displaystyle((\sum_{i=1}^{n}X_{i}(t))_{t\geq 0})$, which has the same
downward transition rate, and an upward transition rate that is at least as
large:

$\displaystyle\displaystyle
k\beta\frac{f(e_{\max}-1)}{e_{\max}-1}\eta(\mathcal{H},m)$
$\displaystyle\displaystyle\leq k\beta\eta(\mathcal{H},m,f)$
$\displaystyle\displaystyle\leq\beta E(S,f),$

where $\displaystyle
S\mathrel{\mathop{\ordinarycolon}}=\\{i\in\\{1,2,\ldots,n\\}\ |\
X_{i}(t)=0\\}$. Thus, to show Theorem 8.5 it suffices to find a suitable lower
bound for $\displaystyle\mathbb{P}\left(\widehat{\tau}>\frac{\lfloor
r^{-m+1}\rfloor}{2m}\right)$, where $\displaystyle\widehat{\tau}$ is the
hitting time of $\displaystyle 0$ for the process $\displaystyle(Z(t))_{t\geq
0}$. This follows by applying Theorem $\displaystyle 4.1$ of [11] to the
process $\displaystyle(Z(t))_{t\geq 0}$.

##  9 Collective Contagion Models

We now consider the collective contagion models from [8, 15, 17], where
infection only starts spreading within a hyperedge once a threshold number of
infectious nodes in that hyperedge has been reached. As discussed in
subsection 3.2, collective contagion models can be represented by nonlinear
functions of the form $\displaystyle
f(x)\mathrel{\mathop{\ordinarycolon}}=\max\\{0,x-c\\}$ for some $\displaystyle
c>0$, or $\displaystyle
f(x)\mathrel{\mathop{\ordinarycolon}}=c_{2}\mathbbm{1}(x\geq c_{1})$ for some
$\displaystyle c_{1},c_{2}>0$. In these cases it is obvious that the zero-
infection state for the mean field approximation is locally asymptotically
stable (and, indeed, Theorem 6.1 applies). However, because the functions are
not concave, the theory found in Section 8 for the exact model does not
directly apply. Nonetheless, we can still derive similar spectral conditions
for the vanishing of the disease by finding concave functions which serve as
upper bounds for $\displaystyle f$. For instance using $\displaystyle
c_{2}\mathbbm{1}(x\geq c_{1})\leq\frac{c_{2}}{c_{1}}x\mathbbm{1}(x\leq
c_{1})+c_{2}\mathbbm{1}(x\geq c_{1})$, the bounds in Theorem 8.1 and Corollary
8.2 lead to the the following result.

######  Theorem 9.1.

Suppose that $\displaystyle
f(x)\mathrel{\mathop{\ordinarycolon}}=c_{2}\mathbbm{1}(x\geq c_{1})$ for some
$\displaystyle c_{1},c_{2}>0$. Then

$\displaystyle\mathbb{P}\left(\sum_{i}X_{i}(t)>0\right)\leq
n\,i_{0}\,\exp\left(\beta\frac{c_{2}}{c_{1}}\lambda(W)-\delta\right).$

In particular if
$\displaystyle\lambda(W)<\frac{c_{1}}{c_{2}}\frac{\delta}{\beta}$, then the
disease asymptotically vanishes with exponential decay and the extinction time
$\displaystyle\tau$ satisfies

$\displaystyle\mathbb{E}[\tau]\leq\frac{\log
n+1}{\delta-\beta\lambda(W)c_{2}/c_{1}}.$

Likewise note that
$\displaystyle\max\\{0,x-c\\}\leq\frac{e_{\max}-1-c}{e_{\max}-1}x,$ where we
recall that $\displaystyle e_{\max}$ is the largest size of a hyperedge of
$\displaystyle\mathcal{H}$. Hence we deduce the following result.

######  Theorem 9.2.

Suppose that $\displaystyle
f(x)\mathrel{\mathop{\ordinarycolon}}=\max\\{0,x-c\\}$ for some $\displaystyle
c>0$. Then

$\displaystyle\mathbb{P}\left(\sum_{i}X_{i}(t)>0\right)\leq
n\,i_{0}\,\exp\left(\beta\frac{e_{\max}-1-c}{e_{\max}-1}\lambda(W)-\delta\right).$

In particular if

$\displaystyle\lambda(W)<\left(\frac{e_{\max}-1}{e_{\max}-1-c}\right)\frac{\delta}{\beta},$

then the disease asymptotically vanishes with exponential decay and the
extinction time $\displaystyle\tau$ satisfies

$\displaystyle\mathbb{E}[\tau]\leq\frac{\log
n+1}{\delta-\beta\frac{e_{\max}-1-c}{e_{\max}-1}\lambda(W)}.$

##  10 Discussion

In this work we derived several spectral conditions that control the spread of
disease in an SIS model on a hypergraph. The conditions have the general form

$\beta\,\lambda(W)\,c_{f}/\delta<1,$ (11)

where $\displaystyle c_{f}>0$ is a constant depending on the function
$\displaystyle f$ that determines the nonlinear infection rate within a
hyperedge.

We note that in the special case where (i) the hypergraph is an undirected
graph and hence $\displaystyle W$ becomes the binary adjacency matrix, and
(ii) we have linear dependence on the number of infectious neighbors for the
infection rate of a node, so $\displaystyle f$ is the identity function, the
condition (11) reduces to the well-known vanishing spectral condition studied
in, for example, ([11, 18, 21]).

There are two important points to be made about the general form of (11).
First, the hypergraph structure appears only via the presence of the symmetric
matrix $\displaystyle W\in\mathbb{R}^{n\times n}$. Recall that $\displaystyle
W_{ij}$ records the number of times that $\displaystyle i$ and $\displaystyle
j$ both appear in the same hyperedge. Such weighted but _pairwise_ information
is all that feeds into this spectral threshold. On a positive note, this
implies that useful predictions can be made about disease spread on a
hypergraph without full knowledge of the types of hyperedge present and the
distribution of nodes within them. (For example, when collecting human
interaction data it is more reasonable to ask an individual to list each
neighbour and state how many different ways they interact with that neighbour
than to ask an individual to list all hyperedges they take part in.) However,
this observation also raises the possibility that more refined analysis might
lead to sharper bounds, perhaps at the expense of simplicity and
interpretability.

Our second point is that the new vanishing condition (11) neatly separates
three aspects:

(a)

The biologically-motivated infection parameter, $\displaystyle\beta$.

(b)

The interaction structure, captured in $\displaystyle\lambda(W)$.

(c)

The coefficient $\displaystyle c_{f}$ that arises from modelling the nonlinear
infection process. For instance, Theorem 6.1 and Theorem 6.4 have
$\displaystyle c_{f}=f^{\prime}(0)$. In the collective contagion model case
$\displaystyle f(x)=c_{2}\mathbbm{1}(x\geq c_{1})$, Theorem 9.1 indicates that
we can take $\displaystyle c_{f}=c_{2}/c_{1}$.

We may view $\displaystyle\beta$ as an invariant biological constant that
reflects the underlying virulence of the disease and is not affected by human
behaviour. The factor $\displaystyle\lambda(W)$, which arises from the
interaction structure, will be determined by regional and cultural issues,
including population density, age demographics, typical household sizes, and
the nature of prevalent commercial and manufacturing activities.
Interventions, including full or partial lockdowns, could be modeled through a
change in $\displaystyle\lambda(W)$. The third factor, $\displaystyle c_{f}$,
is strongly dependent upon human behaviour and may be adjusted to reflect
individual-based containment strategies such as social distancing, mask
wearing or more frequent hand washing.

This work has focused on modelling, analysis and interpretation at the
abstract level, concentrating on the fundamental question of disease
extinction. Having developed this theory, it would, of course, now be of great
interest to perform practical experiments using realistic interaction and
infection data, with the aim of

  * •

calibrating model parameters,

  * •

testing hypotheses about the appropriate functional form of the infection
rate,

  * •

testing the predictive power of the modeling framework, especially in
comparison with simpler homogeneous mixing and pairwise interaction versions,

  * •

quantifying the effect of different interventions.

## References

  * [1] U. Alvarez-Rodriguez, F. Battiston, G. F. de Arruda, Y. Moreno, M. Perc, and V. Latora, Evolutionary dynamics of higher-order interactions in social networks, Nat. Hum. Behav., (2021).
  * [2] R. M. Anderson and R. M. May, Infectious Diseases of Humans: Dynamics and Control, Oxford University Press, Oxford, 1992.
  * [3] F. Battiston, G. Cencetti, I. Iacopini, V. Latora, M. Lucas, A. Patania, J.-G. Young, and G. Petri, Networks beyond pairwise interactions: Structure and dynamics, Physics Reports, 874 (2020), pp. 1–92.
  * [4] A. R. Benson, R. Abebe, M. T. Schaub, A. Jadbabaie, and J. Kleinberg, Simplicial closure and higher-order link prediction, Proceedings of the National Academy of Sciences, 115 (2018), pp. E11221–E11230.
  * [5] A. R. Benson, D. F. Gleich, and J. Leskovec, Higher-order organization of complex networks, Science, 353 (2016), pp. 163–166.
  * [6] A. Bodó, G. Katona, and P. Simon, SIS epidemic propagation on hypergraphs, Bulletin of Mathematical Biology, 78 (2016), pp. 713–735.
  * [7] A. Bretto, Hypergraph Theory: An introduction, Springer, Berlin, 2013\.
  * [8] G. F. de Arruda, G. Petri, and Y. Moreno, Social contagion models on hypergraphs, Phys. Rev. Res., 2 (2020).
  * [9] E. Estrada, COVID-19 and SARS-CoV-2. Modeling the present, looking at the future, Physics Reports, 869 (2020), pp. 1–51.
  * [10] E. Estrada and J. A. Rodríguez-Velázquez, Subgraph centrality and clustering in complex hyper-networks, Physica A: Statistical Mechanics and its Applications, 364 (2006), pp. 581–594.
  * [11] A. Ganesh, L. Massoulié, and D. Towsley, The effect of network topology on the spread of epidemics, Proceedings - IEEE INFOCOM, 2 (2005), pp. 1455–1466.
  * [12] D. T. Gillespie, Markov Processes: An Introduction for Physical Scientists, Academic Press, Camnbridge, MA, 1991.
  * [13] G. Giordano, F. Blanchini, R. Bruno, P. Colaneri, A. Di Filippo, A.and Di Matteo, and M. Colaneri, Modelling the COVID-19 epidemic and implementation of population-wide interventions in Italy, Nature Medicine, 26 (2020), pp. 855–860.
  * [14] P. Hartman, On the stability in the large for systems of ordinary differential equations, Canadian Journal of Mathematics, 13 (1961), pp. 480–492.
  * [15] I. Iacoponi, G. Petri, A. Barrat, and V. Latora, Simplicial models of social contagion, Nature Communications, 10 (2019).
  * [16] W. O. Kermack and A. G. McKendrick, A contribution to the mathematical theory of epidemics, Proceedings of the Royal Society A, 115 (1927), pp. 700–721.
  * [17] N. W. Landry and J. G. Restrepo, The effect of heterogeneity on hypergraph contagion models, Chaos, 30 (2020).
  * [18] P. V. Mieghem, J. Omic, and R. Kooij, Virus spread in networks, IEEE Transactions on Networking, 17 (2009).
  * [19] R. Montenegro and P. Tetali, Mathematical aspects of mixing times in Markov chains, Foundations and Trends in Theoretical Computer Science, 1 (2006), pp. 237 – 354.
  * [20] F. Verhulst, Nonlinear Differential Equations and Dynamical Systems, Springer-Verlag, Berlin, 1990.
  * [21] Y. Wang, D. Chakrabarti, C. Wang, and C. Faloutsos, Epidemic spreading in real networks: an eigenvalue point of view, Proceedings 22nd International Symposium on Reliable Distributed Systems, (2003).


# A Malmquist-like bias in the inferred areas of diamond caustics and the
resulting bias in inferred time delays for gravitationally lensed quasars

Derek Baldwin MIT Department of Physics
Cambridge, MA 02139, USA  Paul L. Schechter MIT Department of Physics
Cambridge, MA 02139, USA  MIT Kavli Institute for Astrophysics and Space
Research
Cambridge, MA 02139, USA

###### Abstract

Quadruply lensed quasars are visible only when the source quasar lies within
the diamond caustic of the lensing galaxy. This condition creates a Malmquist-
like selection effect in the population of observed quadruply lensed quasars,
increasing the true caustic area. The inferred areas are therefore biased low
by an amount proportional to the square of the uncertainty in the natural
logarithm of the area, $\sigma_{\ln{A}}$. The inferred time delays, $\tau$,
vary with the caustic area of the lensing galaxy as $\sqrt{A}$. As a
consequence, model time delays are biased low, leading to underestimates of
the Hubble constant. Classical Malmquist bias also plays a role, with inferred
magnifications biased low. As magnification is strongly anti-correlated with
the area of the caustic, this classical Malmquist bias has an effect opposite
to, but smaller than, the effect of the area bias. The overall effect is
estimated to be small (median bias of less than 1%), but has the potential to
be larger. Current posterior plots do not include caustic area, which limits
the capability to accurately estimate the bias. With anticipated increases in
sample sizes (in particular with the completion of LSST), correction for these
biases becomes increasingly important.

galaxies: quasars — gravitational lensing: strong, Malmquist bias, time delay
cosmography

##  1 Introduction

The “era of precision cosmology” (e.g. Kaplinghat & Turner 2001) was widely
heralded by cosmologists (but not all; see Bridle et al. 2003) and the
occasional astronomer. Much of the precision in “precision cosmology” can be
attributed to the high degree of homogeneity and small amplitude of
perturbations in the early universe. But in today’s universe the populations
of astronomical objects used for cosmological inference (e.g Cepheid variable
stars, type Ia sypernovae, gravitationally lensed quasars) exhibit wide ranges
of properties that must be taken into full account. These properties introduce
uncertainties that become all the more important with the ever-increasing
precision of our cosmological measurements. The present paper discusses
systematic biases in the measurements of several quantities associated with
gravitational lenses (specifically caustic area, time delay, and
magnification) that depend on these uncertainties.
We consider time delay measurements of quadruply lensed quasars, touted as a
high precision method of measurement for the Hubble constant, $H_{0}$ (Refsdal
1964; Treu & Marshall 2016; Wong et al. 2020). It is well-known that these
systems have a “diamond” or “astroidal111We note that these caustics form
perfect astroids only in special cases. The caustics considered here are
almost, but not quite, stretched astroids.” caustic (Ohanian, 1983), which
gives rise to a “quad” of images only when the source quasar is located inside
this caustic (Finch et al., 2002). We demonstrate that this condition leads to
an analogue of the well-known Malmquist effect (Malmquist, 1922; Binney &
Merrifield, 1998), in which the true luminosities of stars in flux-limited
samples are higher than inferred due to a selection bias in the sample. In §2
we show how a similar selection effect occurs in the study of quadruply-lensed
quasars, leading to a systematic underestimate of the diamond caustic area
with respect to the true mean. In §3, we estimate the magnitude of this bias
from uncertainties in the measurements of the external convergences of seven
systems.
While such an underestimate in the inferred caustic areas of lensing galaxies
might seem benign, we argue in §4 that the area and associated relative time
delay in the system are strongly correlated, and so such an underestimate of
caustic area introduces an underestimate of the time delay for each system. We
claim that this in turn creates a finite underestimate of the Hubble constant
that must be accounted for when making precise measurements. In §5 we
calculate an estimate for the bias in the Hubble constant. Finally, in §6, we
explore the classical Malmquist bias introduced by uncertainties in the
magnification.
The present paper serves to establish the existence of this systematic error
rather than produce a definitive calculation of it (e.g. accounting for
multiple correlated model parameters).

##  2 The Malmquist-like effects for lens caustic areas

In the same way that uncertainties in the derived luminosities of stars cause
the true luminosities of the stars in a magnitude-limited sample to be
brighter than the derived luminosities, the uncertainties in the derived areas
of lens caustics cause the true caustic area to be larger than the derived
caustic areas. This is due to the fact that a quadruply-lensed quasar requires
the source to be within the diamond caustic (Finch et al., 2002). The systems
we see as “quads” are likely to have true astroidal caustic areas larger than
inferred because lensing galaxies with true astroidal caustics that are 1
standard deviation larger than inferred are over-represented and lensing
galaxies with true astroidal caustics that are 1 standard deviation smaller
than inferred are under-represented. Since we select only systems with their
source quasar located inside the diamond caustic, the probability that a
system is included in a sample of “quads” is proportional to the true area of
its caustic. This is illustrated in Figure 1.

Figure 1: Plots of diamond caustics with areas
$\ln{A_{inf}}\pm\sigma_{\ln{A}}$, where $\sigma_{\ln{A}}=\ln{2}$, and randomly
generated source quasar locations. Quads with higher caustic areas are over-
represented in a sample because more source locations fall within the larger
areas. The average caustic area in a sample of quads is $A_{peak}\approx
1.5A_{inf}$.

The derivation here of the size of the effect is similar to the derivation of
Malmquist’s effect in Binney & Merrifield (1998). We introduce a Gaussian
distribution of systems with true caustic areas $A$ about an inferred mean
value $A_{inf}$. Then the contribution to the total number of observed
quadruply-lensed systems, $N$, with respect to area is given by:

$\,dN=\frac{1}{\sqrt{2\pi\sigma_{\ln{A}}^{2}}}n\left(\frac{A}{A_{inf}}\right)e^{\textstyle\textstyle\frac{-(\ln{A}-\ln{A_{inf}})^{2}}{2\sigma_{\ln{A}}^{2}}}\,d(\ln{A})$
(1)

where $n$ is a number density (numbers per steradian) and $\sigma_{\ln{A}}$ is
the standard deviation in $\ln{A}$. The ratio $A/A_{inf}$ gives the
probability of observing a lens with true caustic area $A$ relative to one
with $A_{inf}$. Substituting $\exp{\left(\ln{A}-\ln{A_{inf}}\right)}$ for
$A/A_{inf}$ we have

$\displaystyle\,dN$
$\displaystyle=\frac{1}{\sqrt{2\pi\sigma_{\ln{A}}^{2}}}ne^{\textstyle\frac{-(\ln{A}-\ln{A_{inf}})^{2}}{2\sigma_{\ln{A}}^{2}}+(\ln{A}-\ln{A_{inf}})}\,d(\ln{A})$
(2)
$\displaystyle=\frac{1}{\sqrt{2\pi\sigma_{\ln{A}}^{2}}}ne^{\textstyle\frac{-(\ln{A}-\ln{A_{inf}}-\sigma_{\ln{A}}^{2})^{2}+\sigma_{\ln{A}}^{4}}{2\sigma_{\ln{A}}^{2}}}\,d(\ln{A})$
(3)
$\displaystyle=\frac{1}{\sqrt{2\pi\sigma_{\ln{A}}^{2}}}ne^{\textstyle\frac{\sigma_{\ln{A}}^{2}}{2}}e^{\textstyle\frac{-(\ln{A}-\ln{A_{inf}}-\sigma_{\ln{A}}^{2})^{2}}{2\sigma_{\ln{A}}^{2}}}\,d(\ln{A})$
(4)

which is a new Gaussian distribution with a shifted peak at

$\ln{A_{peak}}=\ln{A_{inf}}+\sigma_{\ln{A}}^{2}\,.$ (5)

We can then define the “logarithmic area bias” that reflects this systematic
error:

$\Delta_{\ln{A}}\equiv\ln{A_{inf}}-\ln{A_{peak}}=-\sigma_{\ln{A}}^{2}\,.$ (6)

The logarithmic mean of true caustic areas for lenses for which the inferred
caustic area is $A_{inf}$ is larger than inferred by $\sigma_{\ln{A}}^{2}$, or
the inferred caustic areas are biased low by $\sigma_{\ln{A}}^{2}$. Similarly,
Malmquist’s effect causes absolute magnitudes to be biased faint,

$\Delta_{{M}}=M_{0}-M_{peak}=1.38\sigma_{M}^{2}\,,$ (7)

where the factor of 1.38 is $\frac{3}{2}\ln{2.5}$. In both cases, there is an
additive bias in the logarithm of the relevant quantity – luminosity in the
Malmquist case and caustic area for the quads – that is proportional to the
square in the uncertainty in the logarithm of that quantity.
In the next section we explore how to estimate $\sigma_{\ln{A}}\,.$

##  3 Estimating the uncertainty in caustic area

###  3.1 Background

The “effective”, 2D lensing potential, $\psi$, is given by

$\psi(\vec{r})=\frac{D_{ds}}{D_{d}D_{s}}\frac{2}{c^{2}}\int_{0}^{\infty}\Phi(D_{d}\vec{r},z)\,dz$
(8)

where $D_{ds},D_{d}$ and $D_{s}$ are the angular diameter distances between
the lens and the source, the distance between the lens and the observe, and
the distance between the observer and the source, respectively. The Newtonian
gravitational potential of the lens is $\Phi$.
Typical of the potentials we consider here is the power-law potential subject
to an external shear field,

$\psi(r,\theta)=br^{\alpha}+\frac{\gamma r^{2}}{2}\cos{2\theta}$ (9)

where _b_ gives the mass scale of the lens, _$\alpha$_ is the power-law slope
of the lens, and _$\gamma$_ is the strength of the external tidal shear.
The deflection of the image is given by:

$\vec{r}-\vec{r}_{s}=\nabla\psi(\vec{r})\,,$ (10)

where $\vec{r}$ is the angular image position and $\vec{r}_{s}$ is the angular
source position (Bourassa & Kantowski, 1975).
Each image has an associated time delay due to its varying path length. We
define the time delay difference $\tau$ to be the largest difference in time
delays between the 4 images in a quadruply-lensed quasar. Treu & Marshall
(2016) and Refsdal (1964) tell us that $\tau\propto 1/H_{0}$.
We later examine external convergences, and so we also consider the surface
mass density of a lensing galaxy, defined by the projection of the lensing
galaxy onto the plane perpendicular to the line-of-sight to the quasar. That
is,

$\Sigma(\vec{r})=\int{\rho(\vec{r},z)}\,dz$ (11)

where $\Sigma$ is the surface mass density, $\vec{r}$ is the radial coordinate
of the projected galaxy, $z$ is the line-of-sight coordinate, and $\rho$ is
the conventional mass density profile of the lensing galaxy.
The dimensionless convergence, $\kappa$, is defined as a scaled version of the
surface mass density,

$\kappa(\vec{r})=\frac{\Sigma}{\Sigma_{cr}}\,,$ (12)

where $\Sigma_{cr}$ is the critical surface mass density needed to produce
multiple images,

$\Sigma_{cr}=\frac{c^{2}}{4\pi G}\frac{D_{s}}{D_{d}D_{ds}}\,.$ (13)

Here $c$ and $G$ are the conventional speed of light and Newtonian
gravitational constant, respectively, and $\kappa$ satisfies a dimensionless
version of Poisson’s equation.

###  3.2 Dependence of caustic area on external convergence

The mass sheet degeneracy is well-known in the literature (see Falco et al.
1985). Here we explore some of its consequences with respect to caustic areas
and time delays.
We follow the development of Falco et al. (1985) and Schneider & Sluse (2013).
We consider a convergence of the form

$\kappa^{\prime}(\vec{r})=\lambda\kappa(\vec{r})+(1-\lambda)\,,$ (14)

where $\kappa$, as usual, is the convergence characterizing the mass
distribution of the lensing galaxy and $\lambda$ is a dimensionless quantity.
We now define the quantity

$\kappa_{ext}=1-\lambda$ (15)

which is called the external convergence in the literature. This makes the
convergence

$\kappa^{\prime}(\vec{r})=(1-\kappa_{ext})\kappa(\vec{r})+\kappa_{ext}\,.$
(16)

According to Schneider & Sluse (2013), if we scale the source plane
coordinates by $1-\kappa_{ext}\,,$ the observables in the image plane remained
unchanged. This means that if we model a galaxy without including the external
convergence and then include the external convergence after, the corrected
coordinates $\vec{x}_{corr}$ of the caustic curve are calculated from the
inferred coordinates $\vec{x}_{inf}$ by the equation:

$\vec{x}_{corr}=\frac{\vec{x}_{inf}}{1-\kappa_{ext}}\,.$ (17)

We are left with the corrected caustic area

$A_{corr}=\frac{A_{inf}}{(1-\kappa_{ext})^{2}}\,.$ (18)

###  3.3 Markov chain Monte-Carlo simulations

Logarithmic caustic area has not, until now, seemed sufficiently important to
warrant its inclusion in lens modelers’ posterior plots (e.g. Figure 2 of
Birrer et al. 2020). In its absence, we must infer $\sigma_{\ln{A}}$ from the
uncertainties in other derived parameters. We urge investigators to include it
in their MCMC simulations and posterior plots, obviating the need for
calculations like those that follow.

###  3.4 Caustic area as a function of model parameters

We start with generic model parameters, $\alpha_{i}.$ The standard deviation
in the natural logarithm of the area is given by

$\sigma_{\ln{A}}^{2}=\sum_{i}\sum_{j}\frac{\partial\ln{A}}{\partial\alpha_{i}}\frac{\partial\ln{A}}{\partial\alpha_{j}}\sigma_{ij}^{2}\,,$
(19)

where $\sigma_{ij}^{2}$ is the covariance between $i$ and $j$.
The covariances $\sigma_{ij}^{2}=\rho(i,j)\sigma_{i}\sigma_{j}$, where
$\rho(i,j)$ is the correlation coefficient, can be estimated from published
posterior plots, but $A(\alpha_{1},...,\alpha_{i})$ is needed to proceed
further. In what follows we examine this functional dependence for a few
simple potentials, giving a sense of how it might vary for more realistic
models.

####  3.4.1 Singular isothermal potential in an external shear field

Finch et al. (2002) give an expression for the caustic area of a singular
isothermal potential in an external shear field,

$\psi(r,\theta)=br+\frac{\gamma r^{2}}{2}\cos{2\theta}\,,$ (20)

and find

$A=\frac{3\pi}{2}b^{2}\frac{\gamma^{2}}{1-\gamma^{2}}\,.$ (21)

This is one of the few lensing potentials for which an exact analytic solution
is possible and, by virtue of this fact, it avoids much of the complication
that follows. It is a limiting case of the next potential we consider.

####  3.4.2 Singular isothermal elliptical potential in an external shear
field

Luhtaru et al. (2021) consider the SIEP + XS model,

$\psi(r,\theta)=b\sqrt{q_{pot}x^{2}+\frac{y^{2}}{q_{pot}}}-\frac{\gamma}{2}(x^{2}-y^{2})\,.$
(22)

They arrive at the approximate expression

$A\approx\frac{3\pi b^{2}}{2}(\gamma+2\eta)^{2}\,,$ (23)

where $\eta$ is the semi-ellipticity, $q_{pot}=(1-\eta)/(1+\eta)$. We can
rewrite this in terms of the “effective quadrupole,”

$\Gamma_{eff}=\frac{\gamma+\eta}{1+\gamma\eta}\,,$ (24)

which is tightly constrained by the elongation of the image configuration. The
parameter $b$, roughly the Einstein radius, $\theta_{E}$, is also well
constrained. This gives:

$A\approx\frac{3\pi\theta_{E}^{2}}{2}\left(2\Gamma_{eff}-\gamma\right)^{2}\,.$
(25)

The only poorly constrained parameter is then $\gamma$. Calculation of the
uncertainty in $\ln{A}$ is straightforward because we need not consider
covariances.

####  3.4.3 Power-law + external shear potential

We are interested in the power-law + external shear potential, given by
equation (9). We calculate the approximate area of the diamond caustic for
this case in Appendix A. For this potential the Einstein radius is given by
$\theta_{E}=(b\alpha)^{1/(2-\alpha)}$, which gives

$A\approx\frac{3\pi}{2}\theta_{E}^{2}\gamma^{2}\left(\frac{1}{1-\gamma^{2}}\right)^{\frac{1}{2-\alpha}}\,.$
(26)

As noted previously, the Einstein radius is very well determined for any given
system and so is relatively constant even as the slope is allowed to vary.
The covariance between the slope and external shear will enter the calculation
of the uncertainty in the area. Some systems, such as PG1115+080 and
HE0435-1223, have rings that can, in the absence of external convergence,
determine the power-law slope. Thus one would expect a weak correlation
between the slope and the external shear, which is supported by the posterior
plots in Figures 8 and 13 of Chen et al. (2019). Other systems, however, do
not have these rings to constrain the slope (e.g. WFI2033-4723). The external
shear and power-law slope are thus more strongly correlated for these systems,
which is reflected in the corner plot in Rusu et al. (2020) Figure 7.
The complications arising from correlated parameters leave us with no easy way
to proceed with an estimate in area uncertainties, especially in the absence
of published correlation coefficients and varying circumstances for each
individual system. We thus forego such a calculation while once more urging
modelers to include area in their posterior plots to permit straightforward
calculations.

###  3.5 A special case: caustic area as a function of external convergence

Most investigators estimate the external convergence from counts of galaxies
near the line of sight to lensed quasars. Were this strictly true, it would be
possible to determine that part of the the uncertainty in the caustic area due
to uncertainty in the uncertainty in $\kappa_{ext}$ from equation (18). But in
practice modelers (e.g. Rusu et al. 2020; Chen et al. 2019) use the external
shear derived from the lens model as an additional constraint on the external
convergence.
While this makes sense – there is mass surface density associated with the
same groups and clusters that produce the external shear – it introduces a
covariance, not explicitly presented, between the two parameters. We
nonetheless proceed to calculate caustic area uncertainties on the assumption
that caustic areas do not depend on other model parameters in order to obtain
an estimate of the effect.
We examine the seven systems collected in Birrer et al. (2020) Figure 6. Using
equation (18) we find the uncertainties in the logarithms of the areas due to
the external shear as

$\displaystyle\sigma_{\ln{A}}$
$\displaystyle=\frac{\,d(\ln{A})}{\,d\kappa_{ext}}\sigma_{\kappa_{ext}}$ (27)
$\displaystyle=\frac{2}{1-\kappa_{ext}}\sigma_{\kappa_{ext}}\,.$ (28)

The confidence intervals for $\kappa_{ext}$ are all lopsided, which is not
consistent with our simplifying assumption of a Gaussian distribution. We
approximate a Gaussian distribution by averaging the two bounds of the
confidence interval to get an ad-hoc standard deviation. The approximate
logarithmic area biases for these seven systems are listed in Table 1.

Table 1: External convergences and area biases for the seven systems in Birrer et al. (2020) Figure 6 System | $\kappa_{ext}$ | $\Delta_{\ln{A}}$
---|---|---
B1608+656 | $0.103^{+0.084}_{-0.045}$ | $-2.1\cdot 10^{-2}$
RXJ1131-1231 | $0.069^{+0.043}_{-0.026}$ | $-5.5\cdot 10^{-3}$
HE0435-1223 | $0.004^{+0.033}_{-0.021}$ | $-2.8\cdot 10^{-3}$
SDSS1206+4332 | $-0.004^{+0.036}_{-0.022}$ | $-3.3\cdot 10^{-3}$
WFI2033-4723 | $0.059^{+0.077}_{-0.044}$ | $-1.7\cdot 10^{-2}$
PG1115+080 | $-0.006^{+0.032}_{-0.021}$ | $-6.4\cdot 10^{-3}$
DES0408-5354 | $-0.040^{+0.038}_{-0.024}$ | $-3.6\cdot 10^{-3}$

###  3.6 Combining the effects of uncertainties in two uncorrelated parameters

Suppose that two uncorrelated model parameters, call them $f$ and $g$,
contribute to the uncertainty in caustic area. We take the propagating
uncertainty in the logarithm of the area due to the uncertainty in $f$ to be
$\sigma_{\ln{A},f}$. That is,

$\sigma_{\ln{A},f}=\frac{\partial\ln{A}}{\partial f}\sigma_{f}\,.$ (29)

We similarly define the uncertainty in the log of the area due to the
uncertainties in $g$ to be $\sigma_{\ln{A},g}$. We start with equation (1),

$\,dN\propto\Big{(}\frac{A}{A_{inf}}\Big{)}e^{\textstyle\frac{-(\ln{A}-\ln{A_{inf}})^{2}}{2\sigma_{\ln{A}}^{2}}}\,d(\ln{A})\,.$
(30)

According to equation (19), we can express $\sigma_{\ln{A}}^{2}$ as

$\sigma_{\ln{A}}^{2}=\sigma_{\ln{A},f}^{2}+\sigma_{\ln{A,g}}^{2}$ (31)

for uncorrelated parameters $(\alpha_{1},\alpha_{2})=(f,g)$, because the
covariance of two uncorrelated parameters is zero. Then

$\displaystyle\,dN\propto\Big{(}\frac{A}{A_{inf}}\Big{)}e^{\textstyle\frac{-(\ln{A}-\ln{A_{inf}})^{2}}{2(\sigma_{\ln{A},f}^{2}+\sigma_{\ln{A},g}^{2})}}\,d(\ln{A})\,.$
(32)

The algebra in this scenario is the exact same as in §2 upon the substitution
$\sigma_{\ln{A}}^{2}=\sigma_{\ln{A},f}^{2}+\sigma_{\ln{A},g}^{2}$, so we omit
it. This gives an area bias of

$\Delta_{\ln{A}}^{both}=A_{inf}-A_{peak}=-(\sigma_{\ln{A},f}^{2}+\sigma_{\ln{A,g}}^{2})$
(33)

for uncorrelated $f$ and $g$.

###  3.7 Combining the effects of uncertainties in two perfectly correlated
parameters

We now suppose that $f$ and $g$ are perfectly correlated parameters, following
the exact same procedure. We again expand the $\sigma_{\ln{A}}^{2}$ term
according to equation (19), with $(\alpha_{1},\alpha_{2})=(f,g)$. We then get

$\displaystyle\sigma_{\ln{A}}^{2}$ $\displaystyle=\left(\frac{\partial
A}{\partial f}\right)^{2}\sigma_{f}^{2}+2\frac{\partial A}{\partial
f}\frac{\partial A}{\partial g}\sigma_{fg}^{2}+\left(\frac{\partial
A}{\partial g}\right)^{2}\sigma_{g}^{2}$ (34)
$\displaystyle=\left(\frac{\partial A}{\partial
f}\right)^{2}\sigma_{f}^{2}+2\frac{\partial A}{\partial f}\frac{\partial
A}{\partial g}\sigma_{f}\sigma_{g}+\left(\frac{\partial A}{\partial
g}\right)^{2}\sigma_{g}^{2}$ (35)

in the case of perfect correlation. Putting this in terms of
$\sigma_{\ln{A},f}$, defined in equation (29), and $\sigma_{\ln{A},g}$, we are
left with

$\,dN\propto\left(\frac{A}{A_{inf}}\right)e^{\textstyle\frac{-(\ln{A}-\ln{A_{inf}})}{2(\sigma_{\ln{A},f}+\sigma_{\ln{A},g})^{2}}}\,.$
(36)

The algebra in this case is again the same as in §2 with the substitution
$\sigma_{\ln{A}}=\sigma_{\ln{A},f}+\sigma_{\ln{A},g}$, so we don’t repeat it.
We arrive at

$\ln{A_{peak}}=\ln{A_{inf}}+(\sigma_{\ln{A},f}+\sigma_{\ln{A},g})^{2}$ (37)

and

$\Delta_{\ln{A}}^{both}=-(\sigma_{\ln{A},f}+\sigma_{\ln{A},g})^{2}\,.$ (38)

##  4 The effect of caustic areas on time delays

###  4.1 The effect of variable power-law slope on caustic areas and time
delays

We examine 3 different systems that figure prominently in the TDCOSMO (Birrer
et al., 2020) results: HE0435-1223, WFI2033-4723, and PG1115+080, and then
analyze how varying diamond caustic areas influence the time delay. We adopt
the power-law + external shear lensing potential of equation (9), as
implemented in Keeton’s lensmodel program (Keeton, 2001).

####  4.1.1 Methodology

Lensmodel’s times are non-dimensionalized, and the program assigns a
standardized redshift to the lens and source in the absence of a specification
(which is the case here).

  1. 1.

For each of the 3 systems, we find the lens model that best fits the
observations (image/galaxy position, image fluxes, etc.), all with a power-law
slope ($\alpha$) value of 1.

  2. 2.

We then create fictitious “perfect twin systems” with image positions, fluxes,
etc. altered to exactly match the simulated values from lensmodel. This
removes uncertainties due to image/galaxy positions so as to isolate the
effects of changing $\alpha$ on area.

  3. 3.

We then create new models for each system with set $\alpha$ values of
0.5,0.9,1.0,1.1, and 1.5. Each new model is generated with the “perfect twin”
data values for the positions, fluxes, etc. This creates different-sized
diamond caustics with areas given by equation (26).

  4. 4.

We then find the time delays for each of these modeled systems. The caustic
areas are calculated numerically from the curves generated by lensmodel.

####  4.1.2 Results

The results for each of the 3 systems are contained below in Table 2 and
Figure 2.

Table 2: Areas and Time Delays for HE0435-1223, WFI2033-4723,
and PG1115+080 | HE0435-1223 | WFI2033-4723 | PG1115+080
---|---|---|---
$\alpha$ | $A$ | $\tau$ | $A$ | $\tau$ | $A$ | $\tau$
(Number) | $(\prime\prime)^{2}$ | (Time) | $(\prime\prime)^{2}$ | (Time) | $(\prime\prime)^{2}$ | (Time)
0.5 | 0.0913 | 0.3544 | 0.305 | 0.9740 | 0.181 | 0.8048
0.9 | 0.0491 | 0.2605 | 0.170 | 0.7278 | 0.100 | 0.5975
1.0 | 0.0405 | 0.2370 | 0.140 | 0.6640 | 0.083 | 0.5444
1.1 | 0.0328 | 0.2134 | 0.113 | 0.5994 | 0.067 | 0.4909
1.5 | 0.0100 | 0.1187 | 0.036 | 0.3354 | 0.021 | 0.2739
Figure 2: Log-Log Plots for Caustic Area ($A$) vs. Time Delay ($\tau$) with
lines of best fit. Each line of best fit has slope $\frac{1}{2}\,.$

We see in Figure 2 that all three cases have a straight line log-log plot with
slope $0.5\pm.01$. Without further proof, we shall assume that the time delay
follows the relationship $\tau\propto\sqrt{A}$.

###  4.2 The effects of external convergences on time delays

For the specific case of external convergence, it is possible to show
analytically that the time varies as the square root of the area. Birrer et
al. (2016) tell us that the true time delay scales inversely with
$1-\kappa_{ext}$. That is,

$\tau_{corr}=\frac{\tau_{inf}}{1-\kappa_{ext}}\,,$ (39)

where $\tau_{inf}$ is the inferred time delay, calculated ignoring line-of-
sight mass effects.
Comparing with equation (18), we see that

$\tau\propto\sqrt{A}$ (40)

when varying $\kappa_{ext}$, the same relationship found in §4.1.2.
Together, these results allows us to use the biases in the area calculated in
§3 to calculate the “logarithmic delay bias”, given by

$\displaystyle\Delta_{\ln{\tau}}\equiv\frac{1}{2}\Delta_{\ln{A}}$ (41)

or half of the logarithmic area bias. The delay biases for the three systems
are given in Table 3.

##  5 The effect of area bias on the Hubble constant

The negative bias in caustic area produces a negative bias in time delay. The
Hubble constant is proportional to the ratio of the modeled time delay to the
observed time delay. That is,

$\frac{H_{0}^{inferred}}{H_{0}^{assumed}}=\frac{\tau_{mod}}{\tau_{obs}}\,,$
(42)

where $\tau_{obs}$ is the observed time delay and $\tau_{mod}$ is the modelled
time delay. If the modeled time delays are indeed biased low, the inferred
Hubble constant will also be biased low.
Defining the bias in the logarithm of the Hubble constant to be:

$\Delta_{\ln{H_{0}}}=\ln{H_{0}}-\ln{H_{0}^{inf}}\,,$ (43)

we have

$\Delta_{\ln{H_{0}}}=\Delta_{\ln{\tau}}\,.$ (44)

The logarithm of the Hubble constant is biased by the same amount as the
modelled time delay.
We can use this result to roughly estimate percentage errors in the Hubble
constant using the area biases in Table 1. These estimates are included in
Table 3 below.

Table 3: Percent errors in the Hubble constant due to uncertainties in $\kappa_{ext}$ System | $\Delta_{\ln{A}}$ | $\Delta_{\ln{\tau}}$ | % Error in $H_{0}$
---|---|---|---
B1608+656 | $-2.1\cdot 10^{-2}$ | $-1.1\cdot 10^{-2}$ | 1.1%
RXJ1131-1231 | $-5.5\cdot 10^{-3}$ | $-2.8\cdot 10^{-3}$ | 0.3%
HE0435-1223 | $-2.8\cdot 10^{-3}$ | $-1.4\cdot 10^{-3}$ | 0.1%
SDSS1206+4332 | $-3.3\cdot 10^{-3}$ | $-1.7\cdot 10^{-3}$ | 0.2%
WFI2033-4723 | $-1.7\cdot 10^{-2}$ | $-8.5\cdot 10^{-3}$ | 0.9%
PG1115+080 | $-6.4\cdot 10^{-3}$ | $-3.2\cdot 10^{-3}$ | 0.3%
DES0408-5354 | $-3.6\cdot 10^{-3}$ | $-1.8\cdot 10^{-3}$ | 0.2%

These biases are in all cases quite small compared to the random uncertainties
in the reported measurements of $H_{0}$. But they should be treated as lower
limits, as we have not assessed the effects of other parameters. In the
absence of covariances, their inclusion would increase the errors. Taking
proper account of the covariances, however, might conceivably decrease the
errors. We once again urge authors to include caustic areas in their
posteriors to accurately capture the true magnitude of the bias, as the
calculations become increasingly muddy when including more parameters and
their correlations. We also caution authors that under- and/or over-estimates
in parameter uncertainties can greatly affect the magnitudes of the time and
area biases, as they grow as the _square_ of said uncertainties.

##  6 The effect of magnification on selection: Malmquist redux

###  6.1 Malmquist bias in mean magnification

In this section we examine the actual Malmquist bias by considering selection
effects caused by uncertainties in the inferred magnification. We start with
the Gaussian distribution of systems

$\,dN\propto\Big{(}10^{\textstyle\beta\ln{\mu}}\Big{)}e^{\textstyle\frac{-(\ln{\mu}-\ln{\mu_{0}})^{2}}{2\sigma_{\ln{\mu}}^{2}}}\,d(\ln{\mu})\,,$
(45)

where $\mu$ is the magnification and $\sigma_{\ln{\mu}}$ is the uncertainty in
the logarithm of magnification. The $10^{\beta\ln{\mu}}$ is is meant to
represent the increase in the surface density of quasars with decreasing flux,
with a rate controlled by the generic parameter $\beta$. A reasonable value
might be $\beta=0.3$, roughly half the Euclidean expectation of
$10^{0.6(2.5\log{\mu})}$, but we leave it to investigators to choose a value
appropriate to the case at hand. It has the same effect as the $A/A_{inf}$
term in §2, acting as a selection function. This expression is equivalent to

$\,dN\propto\left(\frac{\mu}{\mu_{0}}\right)^{2.3\beta}e^{\textstyle\frac{-(\ln{\mu}-\ln{\mu_{0}})^{2}}{2\sigma_{\ln{\mu}}^{2}}}\,d(\ln{\mu})\,,$
(46)

if $\ln{10}$ is approximated as 2.3. We can then complete the square in the
exponent and get the expression

$\,dN\propto
e^{\textstyle\frac{-(\ln{\mu}-\ln{\mu_{0}}-2.3\beta\sigma_{ln{\mu}}^{2})^{2}}{2\sigma_{\ln{\mu}}^{2}}}\,d(\ln{\mu})\,.$
(47)

This is another Gaussian distribution with a shifted average of

$\ln{\mu_{peak}}=\ln{\mu_{0}}+2.3\beta\sigma_{\ln{\mu}}^{2}\,.$ (48)

We can define yet another bias, the “magnification uncertainty
bias222Magnification uncertainty bias is not to be confused with
“magnification bias”, the term introduced by Turner (1980) to explain why the
ratio of lensed quasars to unlensed quasars is higher at brighter apparent
magnitudes”, as

$\Delta_{\ln{\mu}}=\ln{\mu_{0}}-\ln{\mu_{peak}}=-2.3\beta\sigma_{\ln{\mu}}^{2}\,.$
(49)

This has the same form as the biases calculated above, and so carries the same
unfortunate computational complications. We therefore urge that modelers
include magnification in addition to caustic area in their posterior plots so
as to accurately compute the uncertainties in the magnification.

###  6.2 Magnification and caustic area

Magnification is well known to be related to caustic area (e.g. Finch et al.
2002). The exact relation, however, changes depending on the model used. For
lenses with superposed external convergence,

$A\propto\mu^{-1}\,.$ (50)

This is the case because the areas of both the caustic and source are
multiplied by the same factor. We also claim without proof that this is a good
approximation for lens profiles similar to equation (9), at least for slopes
close to unity, as shown in Figure 3. This plot was created using the same
method in §4.1, but this time with average magnification of the 4 images,
$\mu$, instead of $\tau$.

Figure 3: Log-Log Plots for Caustic Area ($A$) vs. Magnification ($\mu$) with
lines of best fit. Each line of best fit has slope $-1\,.$

In constrast, for the singular isothermal sphere with external shear, Finch et
al. (2002) find

$A\propto\mu^{-2}\,.$ (51)

In both cases, the magnification is strongly correlated with the caustic area,
and so we cannot compute the bias in the time delay as if they were
uncorrelated.
For simplicity, we assume that the magnification and area are perfectly
correlated. We proceed exactly as in §2, but adding the $10^{\beta\ln{\mu}}$
term in front of the Gaussian. With this term we have

$\displaystyle\,dN$
$\displaystyle\propto\Big{(}10^{\beta\ln{\mu}}\Big{)}\left(\frac{A}{A_{inf}}\right)e^{\textstyle\frac{-(\ln{A}-\ln{A_{inf}})^{2}}{2\sigma_{\ln{A}}^{2}}}\,d(\ln{A})$
(52)
$\displaystyle\propto\left(\frac{\mu}{\mu_{0}}\right)^{2.3\beta}\left(\frac{A}{A_{inf}}\right)e^{\textstyle\frac{-(\ln{A}-\ln{A_{inf}})^{2}}{2\sigma_{\ln{A}}^{2}}}\,d(\ln{A})\,.$
(53)

There are two cases to consider, $A\propto\mu^{-1}$ and $A\propto\mu^{-2}$. We
first examine the former. Substituting in for the magnification in equation
(53), we get

$\,dN\propto\left(\frac{A}{A_{inf}}\right)^{1-2.3\beta}e^{\textstyle\frac{-(\ln{A}-\ln{A_{inf}})^{2}}{2\sigma_{\ln{A}}^{2}}}\,d(\ln{A})\,,$
(54)

from which we arrive at a logarithmic area bias of

$\Delta_{\ln{A}}^{comb}=-(1-2.3\beta)\sigma_{\ln{A}}^{2}\,,$ (55)

or, for our proposed value of $\beta=0.3$,

$\Delta_{\ln{A}}^{comb}=-0.31\sigma_{\ln{A}}^{2}\,.$ (56)

For the other case, $A\propto\mu^{-2}$, the calculation is the exact same
following the substitution. We then get

$\Delta_{\ln{A}}^{comb}=-(1-1.15\beta)\sigma_{\ln{A}}^{2}\,,$ (57)

and

$\Delta_{\ln{A}}^{comb}=-0.655\sigma_{\ln{A}}^{2}\,.$ (58)

for $\beta=0.3$. As before, these biases have an associated logarithmic time
bias with magnitude one-half the area bias.
We note that the inclusion of the magnification works _against_ the area bias
and reduces it. The exact magnitude of the reduction is case-dependent, as the
relationship between the area and magnification varies for each case. While we
believe that $A\propto\mu^{-1}$ is a good approximation for lens profiles of
the form of equation (9), we suggest that authors still investigate the area-
magnification relation for their systems of interest to accurately measure
this effect.

##  7 Conclusion

Due to the nature of the selection of quadruply lensed quasars in our analysis
of lensing systems, it is inevitable that, on average, known quads have
associated lensing galaxies with inferred caustic areas biased low with
respect to the true caustic areas, because the cross section for quadruplicity
is proportional to caustic area. A Malmquist bias due to the limiting
magnitude of the survey produces a smaller bias in the caustic area that has
the opposite sign. While it is in principle possible to compute an approximate
correction for this selection effect from the uncertainties reported for the
model parameters, in practice it requires too many simplifying assumptions. We
therefore suggest that investigators report uncertainties in caustic areas in
their analysis via MCMC simulations.
We claim that there is a straightforward relationship between the caustic area
of a given system and its associated time delay, namely $\tau\propto\sqrt{A}$.
A systematic bias in the caustic areas therefore entails a bias in the time
delays measured in a given system. The Hubble constant is inversely
proportional to the derived time delays in these systems, and so a bias low in
area introduces a systematic underestimate in the value of the Hubble
constant. We therefore claim that current measurements of the Hubble constant
via the study of quadruply lensed quasars are biased low due to the over-
representation of larger caustic areas in a sample.
Using reported uncertainties in estimated external convergences, we compute
biases of less than 1%. Uncertainties in other model parameters may increase
this. Of greatest concern, caustic area bias puts a premium on neither under-
nor over-estimating uncertainties, as bias scales with the square of
uncertainty.

ACKNOWLEDGEMENTS

We thank Drs. Philip Marshall, Sherry Suyu and Thomas Collett for comments on
the manuscript. DB thanks the MIT UROP program for its support.

##  Appendix A Approximate caustic area for a power law + external shear
potential

Here we follow closely the development in Finch et al. (2002) to derive an
approximate analytic expression for astroidal caustic areas. We begin with the
lensing potential:

$\psi(r,\theta)=br^{\alpha}+\frac{\gamma r^{2}}{2}\cos{2\theta}\,.$ (A1)

It will be easier to work with this in Cartesian form:

$\displaystyle\psi(x,y)$
$\displaystyle=b\left(\sqrt{x^{2}+y^{2}}\right)^{\alpha}+\frac{\gamma\left(x^{2}+y^{2}\right)}{2}\cos{2\arctan{\frac{y}{x}}}$
(A2)
$\displaystyle=b\left(\sqrt{x^{2}+y^{2}}\right)^{\alpha}+\frac{\gamma\left(x^{2}+y^{2}\right)}{2}\left(\frac{x^{2}}{x^{2}+y^{2}}-\frac{y^{2}}{x^{2}+y^{2}}\right)$
(A3)
$\displaystyle=b\left(\sqrt{x^{2}+y^{2}}\right)^{\alpha}+\frac{\gamma}{2}\left(x^{2}-y^{2}\right)\,.$
(A4)

The astroidal caustic is where the inverse magnification matrix,

$\mu^{-1}=\left(1-\frac{\partial^{2}\psi}{\partial
x^{2}}\right)\left(1-\frac{\partial^{2}\psi}{\partial
y^{2}}\right)-\left(\frac{\partial^{2}\psi}{\partial x\partial y}\right)$ (A5)

disappears.
As for the shape of the caustic, according to An & Evans (2006), the “series
[for the area] are also Taylor-series expansions with respect to $\gamma$ at
$\gamma=0$. If we truncate the expansion of the caustic after [the first]
term, the expression reduces to the equation of the tetra-cuspi-hypocycloid,
or the astroid.” We make the simplifying assumption that $\gamma<<1$, so as to
simplify the expression to that of approximately an astroid. The area of an
astroid is given by

$A=\frac{3\pi}{8}xy\,,$ (A6)

where $x$ is the $x$-intercept of the caustic and $y$ is the $y$-intercept
(Finch et al. 2002). We then only consider the intercepts of the caustic in
this calculation. On the axes, the cross term goes to zero so we can ignore
it. We are then left with:

$\displaystyle\mu^{-1}$
$\displaystyle=\left(1-b\left(\alpha\sqrt{x^{2}+y^{2}}^{\alpha-1}\frac{y^{2}}{(x^{2}+y^{2})^{\frac{3}{2}}}+\alpha(\alpha-1)\sqrt{x^{2}+y^{2}}^{\alpha-2}\frac{x}{\sqrt{x^{2}+y^{2}}}\right)-\gamma\right)$
(A7)
$\displaystyle\times\left(1-b\left(\alpha\sqrt{x^{2}+y^{2}}^{\alpha-1}\frac{x^{2}}{(x^{2}+y^{2})^{\frac{3}{2}}}+\alpha(\alpha-1)\sqrt{x^{2}+y^{2}}^{\alpha-2}\frac{y}{\sqrt{x^{2}+y^{2}}}\right)+\gamma\right)\,.$
(A8)

On the x-axis, it is the second term that goes to zero. Thus the x-intercept
solves the equation:

$\displaystyle 0=1-b\alpha x^{\alpha-2}+\gamma\,,$ (A9)

giving us

$x_{c}=\pm
b^{\frac{1}{2-\alpha}}\left(\frac{\alpha}{1+\gamma}\right)^{\frac{1}{2-\alpha}}\,.$
(A10)

A parallel calculation for the y-intercept gives us

$y_{c}=\pm
b^{\frac{1}{2-\alpha}}\left(\frac{\alpha}{1-\gamma}\right)^{\frac{1}{2-\alpha}}\,.$
(A11)

Remapping back to the source plane via the lens equation,

$\vec{r}-\vec{r_{s}}=\nabla\psi(\vec{r})\,,$ (A12)

we get

$\displaystyle x_{a}$ $\displaystyle=x_{c}\left(1-b\alpha
x_{c}^{\alpha-2}-\gamma\right)$ (A13) $\displaystyle=-2\gamma x_{c}$ (A14)
$\displaystyle=\pm 2\gamma
b^{\frac{1}{2-\alpha}}\left(\frac{\alpha}{1+\gamma}\right)^{\frac{1}{2-\alpha}}\,.$
(A15)

The same process for the y-intercept yields

$\displaystyle y_{c}=\pm 2\gamma
b^{\frac{1}{2-\alpha}}\left(\frac{\alpha}{1-\gamma}\right)^{\frac{1}{2-\alpha}}\,.$
(A16)

By the assumption $\gamma<<1$, we know the area of the caustic is proportional
to the product of these two intercepts. The constant of proportionality for an
astroid shape is $3\pi/8$. Thus we arrive at the final expression

$A=\frac{3\pi}{2}\gamma^{2}b^{\frac{2}{2-\alpha}}\left(\frac{\alpha^{2}}{1-\gamma^{2}}\right)^{\frac{1}{2-\alpha}}\,.$
(A17)

## References

  * An & Evans (2006) An, J. H., & Evans, N. W. 2006, MNRAS, 369, 317, doi: 10.1111/j.1365-2966.2006.10303.x
  * Binney & Merrifield (1998) Binney, J., & Merrifield, M. 1998, Galactic Astronomy (Princeton University Press)
  * Birrer et al. (2016) Birrer, S., Amara, A., & Refregier, A. 2016, J. Cosmology Astropart. Phys, 2016, 020, doi: 10.1088/1475-7516/2016/08/020
  * Birrer et al. (2020) Birrer, S., Shajib, A. J., Galan, A., et al. 2020, A&A, 643, A165, doi: 10.1051/0004-6361/202038861
  * Bourassa & Kantowski (1975) Bourassa, R. R., & Kantowski, R. 1975, ApJ, 195, 13, doi: 10.1086/153300
  * Bridle et al. (2003) Bridle, S. L., Lahav, O., Ostriker, J. P., & Steinhardt, P. J. 2003, Science, 299, 1532, doi: 10.1126/science.1082158
  * Chen et al. (2019) Chen, G. C. F., Fassnacht, C. D., Suyu, S. H., et al. 2019, MNRAS, 490, 1743, doi: 10.1093/mnras/stz2547
  * Falco et al. (1985) Falco, E. E., Gorenstein, M. V., & Shapiro, I. I. 1985, ApJ, 289, L1, doi: 10.1086/184422
  * Finch et al. (2002) Finch, T. K., Carlivati, L. P., Winn, J. N., & Schechter, P. L. 2002, ApJ, 577, 51, doi: 10.1086/342163
  * Kaplinghat & Turner (2001) Kaplinghat, M., & Turner, M. S. 2001, Phys. Rev. Lett., 86, 385, doi: 10.1103/PhysRevLett.86.385
  * Keeton (2001) Keeton, C. R. 2001, arXiv e-prints, astro.  https://arxiv.org/abs/astro-ph/0102341
  * Luhtaru et al. (2021) Luhtaru, R., Schechter, P. L., & de Soto, K. M. 2021, ApJ, 915, 4, doi: 10.3847/1538-4357/abfda1
  * Malmquist (1922) Malmquist, K. G. 1922, Meddelanden fran Lunds Astronomiska Observatorium Serie I, 100, 1
  * Ohanian (1983) Ohanian, H. C. 1983, ApJ, 271, 551, doi: 10.1086/161221
  * Refsdal (1964) Refsdal, S. 1964, MNRAS, 128, 307, doi: 10.1093/mnras/128.4.307
  * Rusu et al. (2020) Rusu, C. E., Wong, K. C., Bonvin, V., et al. 2020, MNRAS, 498, 1440, doi: 10.1093/mnras/stz3451
  * Schneider & Sluse (2013) Schneider, P., & Sluse, D. 2013, A&A, 559, A37, doi: 10.1051/0004-6361/201321882
  * Treu & Marshall (2016) Treu, T., & Marshall, P. J. 2016, A&A Rev., 24, 11, doi: 10.1007/s00159-016-0096-8
  * Turner (1980) Turner, E. L. 1980, ApJ, 242, L135, doi: 10.1086/183418
  * Wong et al. (2020) Wong, K. C., Suyu, S. H., Chen, G. C. F., et al. 2020, MNRAS, 498, 1420, doi: 10.1093/mnras/stz3094