# Looking forwards and backwards: dynamics and genealogies of locally
regulated populations

Alison M. Etheridge Thomas G. Kurtz Department of Statistics Departments of
Mathematics and Statistics Oxford University University of Wisconsin - Madison
24-29 St Giles 480 Lincoln Drive Oxford OX1 3LB Madison, WI 53706-1388 UK USA
<EMAIL_ADDRESS><EMAIL_ADDRESS>http://www.stats.ox.ac.uk/~etheridg/ http://www.math.wisc.edu/~kurtz/ Ian
Letter Peter L. Ralph Department of Statistics Departments of Mathematics and
Biology Oxford University University of Oregon 24-29 St Giles Fenton Hall
Oxford OX1 3LB Eugene, OR 97403-1222 UK USA<EMAIL_ADDRESS><EMAIL_ADDRESS>https://www.stats.ox.ac.uk/~restucci/
https://math.uoregon.edu/profile/plr Terence Tsui Ho Lung Department of
Statistics Oxford University 24-29 St Giles Oxford OX1 3LB UK
<EMAIL_ADDRESS>https://www.maths.ox.ac.uk/people/terence.tsui

###### Abstract

We introduce a broad class of mechanistic spatial models to describe how
spatially heterogeneous populations live, die, and reproduce. Individuals are
represented by points of a point measure, whose birth and death rates can
depend both on spatial position and local population density, defined at a
location to be the convolution of the point measure with a suitable non-
negative integrable kernel centred on that location. We pass to three
different scaling limits: an interacting superprocess, a nonlocal partial
differential equation (PDE), and a classical PDE. The classical PDE is
obtained both by a two-step convergence argument, in which we first scale time
and population size and pass to the nonlocal PDE, and then scale the kernel
that determines local population density; and in the important special case in
which the limit is a reaction-diffusion equation, directly by simultaneously
scaling the kernel width, timescale and population size in our individual
based model.

A novelty of our model is that we explicitly model a juvenile phase. The
number of juveniles produced by an individual depends on local population
density at the location of the parent; these juvenile offspring are thrown off
in a (possibly heterogeneous, anisotropic) Gaussian distribution around the
location of the parent; they then reach (instant) maturity with a probability
that can depend on the local population density at the location at which they
land. Although we only record mature individuals, a trace of this two-step
description remains in our population models, resulting in novel limits in
which the spatial dynamics are governed by a nonlinear diffusion.

Using a lookdown representation, we are able to retain information about
genealogies relating individuals in our population and, in the case of
deterministic limiting models, we use this to deduce the backwards in time
motion of the ancestral lineage of an individual sampled from the population.
We observe that knowing the history of the population density is not enough to
determine the motion of ancestral lineages in our model. We also investigate
(and contrast) the behaviour of lineages for three different deterministic
models of a population expanding its range as a travelling wave: the Fisher-
KPP equation, the Allen-Cahn equation, and a porous medium equation with
logistic growth.

Key words: population model, interacting superprocess, lookdown construction,
porous medium equation, reaction-diffusion equation, travelling waves,
genealogies, Fisher-KPP equation

MSC 2010 Subject Classification: Primary:
Secondary:

###### Contents

  1. 1 Introduction
  2. 2 Model and main results
    1. 2.1 Scaling limits of the population process
    2. 2.2 Ancestral lineages in the scaling limit
  3. 3 Examples and applications
    1. 3.1 Beyond linear diffusion
    2. 3.2 Ancestry in different types of travelling waves
    3. 3.3 Clumping from nonlocal interactions
    4. 3.4 Lineage motion is not uniquely determined by population density
  4. 4 Heuristics
    1. 4.1 The population density
    2. 4.2 Motion of ancestral lineages
  5. 5 The lookdown process
    1. 5.1 Lookdown representation of the model of Definition 2.4
    2. 5.2 Explicit construction of lines of descent
    3. 5.3 Limiting processes for lines of descent
  6. 6 Proofs of convergence for nonlocal models
    1. 6.1 Preliminaries
    2. 6.2 Proof of Theorem 2.10: convergence for the nonlocal process
    3. 6.3 Convergence of some nonlocal equations to classical PDEs
      1. 6.3.1 Reaction–diffusion equation limits
      2. 6.3.2 Porous Medium Equation
  7. 7 Simultaneous scaling with interaction distance
    1. 7.1 Moment bounds for $\rho_{\epsilon}*\eta$
    2. 7.2 Continuity estimates for $\rho_{\epsilon}*\eta$
    3. 7.3 Identification of the limit
  8. 8 Proofs of results for the lookdown process and ancestral lineages
    1. 8.1 Tightness of the Lookdown Process
    2. 8.2 Motion of ancestral lineages
  9. A Markov Mapping Theorem
    1. A.1 Lookdown Generators
  10. B Technical Lemmas
    1. B.1 Constraints on kernel widths
    2. B.2 Tightness of processes

##  1 Introduction

As one takes a journey, long or short, the landscape changes: forests thicken
or thin or change their composition; even in flat plains springtime grasslands
host intergrading mosaics of different types of flowers. The aim of this paper
is to introduce and study a broad class of mechanistic spatial models that
might describe how spatially heterogeneous populations live, die, and
reproduce. Questions that we (start to) address include: How does population
density change across space and time? How might we learn about the underlying
dynamics from genealogical or genetic data? And, how does genetic ancestry
spread across geography when looking back through time in these populations?

Reproduction of individuals naturally leads to spatial branching process
models, including branching random walk, branching Brownian motion, and the
Dawson-Watanabe superprocesses. However, as a result of the branching
assumption (once born, individuals behave independently of one another), a
population evolving according to any of these models will either die out or
grow without bound and, in so doing, can develop clumps of arbitrarily large
density and extent. Our starting point here is an individual-based model of a
single species in continuous space in which birth, death, and establishment
may all depend on local population density as well as on spatial location,
allowing for stable populations through density-dependent feedback.

Although it is often mathematically convenient to assume that individuals
follow Brownian motion during their lifetime, in our model, offspring are
thrown off according to some spatial distribution centred on the location of
the parent and do not subsequently move. This is particularly appropriate for
modelling plant populations, in which this dispersal of offspring around the
parent is the only source of spatial motion.

Often models do not distinguish between juveniles and adults, so, for example,
the number of adults produced by a single parent is determined only by the
degree of crowding at the location of the parent. Although we shall similarly
only follow the adult population, in formulating the dynamics of the models we
shall distinguish between production of juveniles, which will depend upon the
location of the adult, and their successful establishment, which will depend
on the location in which a juvenile lands. The result is that not only the
absolute number, but also the spatial distribution around their parent, of
those offspring that survive to adulthood will depend upon the local
population density.

We shall consider three different classes of scaling limits for our model. The
first yields a class of (generalised) superprocesses in which coefficients
governing both the spatial motion and the branching components of the process
can depend on local population density; the second is a corresponding class of
deterministic non-local differential equations; and the third are classical
PDEs. We measure local population density around a point by convolving with a
smooth kernel $\rho(\cdot)$, which may differ for the two stages of
reproduction. When the limiting population process is deterministic, it is a
(weak) solution of an equation of the form

$\partial_{t}\varphi_{t}(x)=r\left(x,\varphi_{t}\right)\mathcal{B}^{*}\left[\varphi_{t}(\cdot)\gamma\big{(}\cdot,\varphi_{t}\big{)}\right](x)+\varphi_{t}(x)F\left(x,\varphi_{t}\right),$
(1.1)

where $\varphi_{t}(x)$ can be thought of as the population density at $x$
(although the limit may be a measure without a density), and $\mathcal{B}^{*}$
is (the adjoint of) a strictly uniformly elliptic second order differential
operator, typically the Laplacian. The dependence of each of the terms $r$,
$\gamma$, and $F$ on $\varphi$ is only through the local density at $x$, e.g.,
$F(x,\varphi)=F(x,\rho\\!*\\!\varphi(x))$. We shall be more specific about the
parameters below.

By replacing $\rho$ by
$\rho^{\epsilon}(\cdot)=\rho(\cdot/\epsilon)/\epsilon^{d}$, we can also scale
the “width” of the region over which we measure local population density. When
the population follows (1.1), we expect that if we take a second limit of
$\epsilon\to 0$, thus scaling the kernels appearing in $r$, $\gamma$, and $F$
and making interactions pointwise, we should recover a nonlinear PDE. We
verify that this is indeed the case in two important examples: a special case
of the porous medium equation with a logistic growth term, in which the
limiting equation takes the form

$\partial_{t}\varphi=\Delta(\varphi^{2})+\varphi(1-\varphi);$ (1.2)

and a wide class of semi-linear PDEs of the form

$\partial_{t}\varphi=\mathcal{B}^{*}\varphi+\varphi F(\varphi),$ (1.3)

which includes the Fisher-KPP equation and the Allen-Cahn equation. Equations
of this form have been studied extensively in the context of spatial ecology
(see for instance Lam and Lou (2023) and Cantrell and Cosner (2004)) and in
many other fields; for instance, Ghosh and Good (2022) derive a stochastic
version of (1.3) to describe abundances of mutant bacteria strains along the
human gut, while Li et al. (2022) study the effects of nonlinear diffusion on
long-term survival of a lattice-based interacting particle system. However, we
do not study the effect of movement of adults, which can additionally affect
the limiting equations: see for instance Holmes et al. (1994) or Potts and
Börger (2023).

It is of interest to understand under what conditions we can replace the two-
step limiting process described above by one in which we simultaneously scale
the kernels and the other parameters in our population model to arrive at the
PDE limit. This is mathematically much more challenging, but we establish such
one-step convergence in cases for which the limit is a classical reaction-
diffusion equation of the form (1.3) with $\mathcal{B}=\Delta$, and $\rho$ is
a Gaussian density. We allow a wide class of reaction terms, $F$, so that the
Fisher-KPP equation (that is equation (1.3) with $\mathcal{B}=\Delta$ and
$F(\varphi)=1-\varphi$) emerges as a special case.

Such results on (one-step) convergence to reaction-diffusion equation limits
have been achieved for a variety of interacting particle systems. Following
the now classical contributions of De Masi et al. (1986); DeMasi and Presutti
(2006); Oelschläger (1985), much of this work has focused on lattice based
models with one particle per site, or on systems with a fixed number, $N$, of
interacting diffusions as $N\to\infty$. For systems of proliferating
particles, as considered for example by Oelschläger (1989); Flandoli et al.
(2019); Flandoli and Huang (2021), an additional challenge (also apparent in
our models), is the control of concentration of particles. We follow
Oelschläger (1989); Flandoli et al. (2019) in considering ‘moderate
interactions’, meaning that the number of individuals in the neighbourhood
over which we measure local population density tends to infinity, whereas
Flandoli and Huang (2021) also consider the situation in which that number
remains finite. We refer to Flandoli and Huang (2021) for a more thorough
literature review, but note that both our model and scaling differ from those
considered in the body of work discussed there: whereas in those settings, the
only scalings are the number of particles in the system and the size of the
neighbourhood over which individuals interact with one another, in keeping
with the vast literature on continuous state branching models, we also scale
time and so must ensure that births are adequately compensated by deaths to
prevent the population from exploding.

The history of a natural population is often only accessible indirectly,
through patterns of genetic diversity that have been laid down; from genetic
data, one can try to infer the genealogical trees that relate individuals in a
sample from the population, and these have been shaped by its history (see
e.g., Neigel and Avise, 1993; Kelleher et al., 2019). It is therefore of
interest to establish information about the distribution of genealogical trees
under our population model, which we do with a lookdown construction. Lookdown
constructions were first introduced in Donnelly and Kurtz (1996) to provide a
mechanism for retaining information about genealogical relationships between
individuals sampled from a population evolving according to the Moran model
when passing to the infinite population limit. Since then, they have been
extended to a wide range of models. Of particular relevance to our work here
are the papers of Kurtz and Rodrigues (2011) and Etheridge and Kurtz (2019),
in which lookdown constructions are provided for a wide variety of population
models, including spatially structured branching processes.

In general, even armed with a lookdown construction, calculation of relevant
statistics of the genealogy remains a difficult question. However, in special
circumstances, some progress can be made. As an illustration, we shall
consider a scenario that has received a great deal of attention in recent
years, in which a population is expanding into new territory as a travelling
wave. In Section 3.2 we shall describe the motion of a single ancestral
lineage relative to three different (deterministic) wavefronts across
$\mathbb{R}^{1}$.

Most work on the topic of “waves” of expanding populations has focused on
models that caricature the classical Fisher-KPP equation with a stochastic
term, i.e.

$dw=\big{(}\Delta w+sw(1-w)\big{)}dt+\sqrt{\frac{\alpha(w)}{N}}W(dt,dx),$

where $W$ is space-time white noise, and $N$ is a measure of the local
population density. The coefficient $\alpha(w)$ is generally taken to be
either $w$, corresponding to a superprocess limit, or $w(1-w)$ giving a
spatial analogue of a Wright-Fisher diffusion. Starting with the pioneering
work of Brunet et al. (2006), a considerable body of evidence has been amassed
to underpin the conjecture that for this, and a wide class of related models,
genealogies converge on suitable timescales in the infinite density limit to a
Bolthausen-Sznitman coalescent. This reflects the fact that, for this
equation, ancestral lineages become trapped in the wavefront, where the growth
rate of the population is highest. Once there, they will experience rapid
periods of coalescence corresponding to significant proportions of individuals
in the front being descended from particularly reproductively successful
ancestors.

If one replaces the logistic growth term of the classical Fisher-KPP equation
with a nonlinearity that reflects cooperative behaviour in the population,
such as

$wF(w)=w(1-w)(Cw-1),$ (1.4)

then, for sufficiently large $C$ (strong cooperation), the nature of the
deterministic wave changes from “pulled” to “pushed”, (Birzu et al., 2018,
2021), and so the genealogies will be quite different from the Fisher-KPP
case. For example, Etheridge and Penington (2022) show that for a discrete
space model corresponding to this nonlinearity with $C>2$, after suitable
scaling, the genealogy of a sample converges not to a Bolthausen-Sznitman
coalescent, but to a Kingman coalescent. The reason, roughly, is that
ancestral lineages settle to a stationary distribution relative to the
position of the wavefront which puts very little weight close to the ‘tip’ of
the wave, so that when ancestral lineages meet it is typically at a location
in which population density is high, where no single ancestor produces a
disproportionately large number of descendants in a short space of time.

The shape of the wave is not determined solely by the reaction term. For
example, as a result of the nonlinear diffusion, for suitable initial
conditions, the solution to the one-dimensional porous medium equation with
logistic growth (1.2) converges to a travelling wave with a sharp cut-off;
i.e., in contrast to the classical Fisher KPP equation, the solution at time
$t$ vanishes beyond $x=x_{0}+ct$ for some constant wavespeed $c>0$ (Kamin and
Rosenau, 2004). As a first step towards understanding what we should expect in
models with nonlinear diffusion, one can ask about the position of an
ancestral lineage relative to the wavefront in the deterministic models. In
Section 3.2 we shall see that in our framework, even with logistic growth, the
nonlinear diffusion corresponding to the porous medium equation results in a
stationary distribution for the ancestral lineage that is concentrated behind
the wavefront, leading us to conjecture that in the stochastic equation the
cooperative behaviour captured by the nonlinear diffusion will also result in
a qualitatively different pattern of coalescence to that seen under the
stochastic Fisher-KPP equation. Indeed, we believe that it should be feasible
to show that in an appropriate limit one recovers a Kingman coalescent.

##### Structure of the paper

In this paper we study scaling limits of spatial population models, obtaining
convergence of both the population process (i.e., the population density as a
function of time, although strictly speaking it is a measure that may not have
a density) and of lineages traced back through such a population. We retain
information about lineages as we pass to the scaling limit by means of a
lookdown construction.

In what follows we first study various scaling limits of the spatial
population process, and then turn our attention to lineages traced back
through these populations. First, in Section 2, we describe the model and the
main results, Theorems 2.10, 2.20, and 2.23. Next, in Section 3, we discuss a
few striking consequences of these results regarding the behavior of
genealogies in traveling waves, the appearance of periodic “clumps” in
seemingly homogeneous population models, and identifiability of the underlying
dynamics from a stationary population profile. In Section 4, we provide
heuristic explanations of why the theorems ought to be true, and some key
ideas behind them, and in Section 5 we define and discuss the lookdown
construction. Proofs of the results begin in Section 6, which proves results
for population models with nonlocal interactions, while Section 7 gives the
more difficult proof for the case when interaction distances also go to zero
in the limit. Finally, Section 8 gives proofs for convergence of the lookdown
process and the associated results for the motion of lineages. The Appendix
contains a few more technical and less central lemmas. The results are
illustrated in a few places with individual-based simulations, made using SLiM
(Haller and Messer, 2019), but these are provided for visualization and we do
not embark on numerical study.

##### Acknowledgements:

Thanks go to Gilia Patterson for identifying the “clumping” phenomenon, and to
Marcin Bownick and David Levin for useful discussions. AME thanks everyone in
MAPS at Université Paris Cité for their hospitality during the period in which
much of this research took place. AME and PLR also thank the Kavli Institute
for Theoretical Physics for their hospitality and birdwatching opportunities.
PLR was supported by the NIH NHGRI (grant #HG011395), IL by the ANID/Doctorado
en el extranjero doctoral scholarship, grant #2018-72190055, and TTHL by the
EPSRC Centre for Doctoral Training in Mathematics of Random Systems: Analysis,
Modelling and Simulation (EP/S023925/1) the Deutsche Forschungsgemeinschaft
under Germany’s Excellence Strategy, EXC-2047/1-390685813, the Rhodes Trust
and St. John’s College, Oxford.

##  2 Model and main results

Our model is one of individuals distributed across a continuous space which we
shall take to be $\mathbb{R}^{d}$. For applications, $d=1$ or $d=2$ (or even
$d=3$ for cells within the body), but our main results apply more generally.
At time zero, the population is distributed over a bounded region, with
$\mathcal{O}(N)$ individuals per unit area in that region, so the total number
of individuals will also be $\mathcal{O}(N)$. The population changes in
continuous time, and we encode the state of the population at time $t$ by a
counting measure $X(t)$, which assigns one unit of mass to the location of
each individual.

Population dynamics are controlled by three quantities, birth ($\gamma$),
establishment ($r$), and death ($\mu$), each of which can depend on spatial
location and local population density in a way specified below. Each
individual gives birth at rate $\gamma$ to a single (juvenile) offspring,
which is dispersed according to a kernel $q(x,\cdot)$ away from the location
$x$ of the parent. We assume that $q$ is the density of a multivariate
Gaussian, allowing a nonzero mean and anisotropic variance. Both the mean and
covariance of $q$ can change across space, but do not depend on population
density. The offspring does not necessarily survive to be counted in the
population: it “establishes” with probability $r$, or else it dies
immediately. Independently, each individual dies with rate $\mu$.

We aim to capture universal behaviour by passing to a scaling limit.
Specifically, we shall take the “density”, $N$, to infinity, and also scale
time by a factor of $\theta=\theta(N)$, in such a way that defining
$\eta^{N}(t)=X(\theta t)/N$, the process $\\{\eta^{N}(t)\\}_{t\geq 0}$ will
converge to a suitable measure-valued process as $N$ and $\theta$ tend to
infinity, with the nature of the limit depending on how they tend to infinity
together. Evidently, we also need to scale the dispersal kernel if we are to
obtain a nontrivial limit, for which we use $q_{\theta}(x,\cdot)$, the density
of the multivariate Gaussian obtained by multiplying the mean and variance
components of $q(x,\cdot)$ by $1/\theta$.

Birth, establishment, and death can depend on the location of the individual
and the local population density. Since we would like the population density
to scale with $N$, these are functions of $X/N$, i.e., the counting measure
with mass $1/N$ placed at the location of each individual. First consider
birth rates, defined by a nonnegative function
$\gamma(x,m):\mathbb{R}^{d}\times\mathbb{R}_{\geq 0}\to\mathbb{R}_{\geq 0}$ of
location $x$ and local population density $m$. Local population density is
defined as the convolution of $X/N$ with a smooth (non-negative, integrable)
kernel $\rho_{\gamma}(\cdot)$. We write this convolution as
$\rho_{\gamma}\\!*\\!X/N$. Then, when the state of the population is $X$, an
individual at location $x$ gives birth to a single juvenile offspring at rate
$\gamma(x,\rho_{\gamma}\\!*\\!X(x)/N)$. Similarly, the establishment
probability of an offspring at location $y$ is is
$r(y,\rho_{r}\\!*\\!X(y)/N)$, where
$r(y,m):\mathbb{R}^{d}\times\mathbb{R}_{\geq 0}\to[0,1]$ and again
$\rho_{r}\\!*\\!X$ is the convolution of $X/N$ with the smooth kernel
$\rho_{r}$.

We shall write $\mu_{\theta}(x,X/N)$ for the per-capita death rate of mature
individuals in the population. In order for the population density to change
over timescales of order $\theta$, we should like the net per capita
reproductive rate to scale as $1/\theta$. In classical models, in which $r$,
$\gamma$, and $\mu$ are constant, this quantity is simply $r\gamma-\mu$. Here,
because production of juveniles and their establishment are mediated by
population density measured relative to different points, the net reproductive
rate will take a more complicated form. In particular, the total rate of
production of mature offspring by an individual at $x$ will be

$\gamma\big{(}x,\rho_{\gamma}\\!*\\!X(x)/N\big{)}\int
r\big{(}y,\rho_{r}\\!*\\!X(y)/N\big{)}q_{\theta}(x,dy).$ (2.1)

Nonetheless, it will be convenient to define the death rate $\mu_{\theta}$ in
terms of its deviation from $r\gamma$. To this end, we define the death rate
of an individual at $x$ to be

$\displaystyle\mu_{\theta}(x)=\max\left\\{0,r(x,\rho_{r}\\!*\\!X(x)/N)\gamma(x,\rho_{\gamma}\\!*\\!X(x)/N)-\frac{1}{\theta}F(x,\rho_{F}\\!*\\!X(x)/N)\right\\},$
(2.2)

where this equation defines $F(x,m):\mathbb{R}^{d}\times\mathbb{R}_{\geq
0}\to\mathbb{R}$, and $\rho_{F}$ is again a smooth kernel. The function $F$ is
nearly the net per capita reproductive rate, scaled by $\theta$, and would be
equal to it in a nonspatial model; but, as can be seen from (2.1), differs
because an offspring’s establishment probability is measured at their new
location rather than that of their parent. For the most part, we work with $F$
instead of $\mu_{\theta}$.

So, each of the three demographic parameters $r$, $\gamma$, and $F$, depends
on local density, measured by convolution with a smooth kernel, each of which
can be different. As a result, death rate depends (in principle) on population
densities measured in three different ways, so that we could write
$\mu_{\theta}(x)=\mu_{\theta}(x,\rho_{\gamma}\\!*\\!X(x)/N,\rho_{r}\\!*\\!X(x)/N,\rho_{F}\\!*\\!X(x)/N)$.
This may seem unnecessarily complex. However, not only is it natural from a
biological perspective, it also turns out to be convenient for capturing
nontrivial examples in the scaling limit.

######  Remark 2.1.

Although this model allows fairly general birth and death mechanisms, there
are a number of limitations. Perhaps most obviously, to simplify the notation
individuals give birth to only one offspring at a time, although this
restriction could be easily lifted (as in Section 3.4 of Etheridge and Kurtz,
2019). Furthermore, individuals do not move during their lifetime, and the age
of an individual does not affect its fecundity or death rate. Finally, there
is no notion of mating (although limitations on reproduction due to
availability of mates can be incorporated into the birth rate, $\gamma$), so
the lineages we follow will be uniparental. For these reasons, the model is
most obviously applicable to bacterial populations or selfing plants, although
we do not anticipate that incorporation of these complications will change the
general picture.

For each $N$ and $\theta$, we study primarily the process with mass scaled by
$N$ and time scaled by $\theta$,

$\big{(}\eta^{N}_{t}\big{)}_{t\geq 0}:=\big{(}X(\theta t)/N\big{)}_{t\geq 0},$

which takes values in the space of càdlàg paths in
$\mathcal{M}_{F}(\mathbb{R}^{d})$ (the space of finite measures on
$\mathbb{R}^{d}$ endowed with the weak topology). In fact $\eta_{t}^{N}$ will
be a purely atomic measure comprised of atoms of mass $1/N$.

######  Notation 2.2.

Expressions like $\gamma(x,\rho_{\gamma}\\!*\\!\eta(x))$ will appear
repeatedly in what follows. To make formulae more readable, we overload
notation to define

$\gamma(x,\eta):=\gamma(x,\rho_{\gamma}\\!*\\!\eta(x)),$

and similarly write $r(x,\eta)$ for $r(x,\rho_{r}\\!*\\!\eta(x))$, $F(x,\eta)$
for $F(x,\rho_{F}\\!*\\!\eta(x))$, and $\mu_{\theta}(x,\eta)$ for the
expression of equation (2.2). When convenient, we may also suppress the
arguments completely, writing simply $\gamma$, $r$, $F$, and $\mu_{\theta}$
for these quantities.

######  Terminology 2.3.

In our prelimiting model, the population is represented by a point measure in
which each individual is assigned a mass $1/N$. We use the term “population
density” for this process, as it is supposed to measure population size
relative to a nominal occupancy of $N$ individuals per unit area. There is no
implication that the measure representing the population is absolutely
continuous with respect to Lebesgue measure; indeed in the prelimit it is
certainly not.

In summary, at each time $t$, $\eta^{N}_{t}$ is purely atomic, consisting of
atoms of mass $1/N$ (which are the individuals). At instantaneous rate
$\theta\gamma(x,\eta^{N}_{t})N\eta^{N}_{t}(dx)$ an offspring of mass $1/N$ is
produced at location $x$; it disperses to a location $y$ offset from $x$ by an
independent Gaussian random variable with mean $\vec{b}(x)/\theta$ and
covariance matrix $\mathbf{C}(x)/\theta$, and once there establishes
instantaneously with probability $r(y,\eta^{N}_{t})$, or else dies. At
instantaneous rate $\theta\mu_{\theta}(x,\eta^{N}_{t})N\eta^{N}_{t}(dx)$ an
individual at location $x$ dies. Note that the process
$\left(\eta^{N}_{t}\right)_{t\geq 0}$, which records numbers and locations of
adult individuals, is just a scaled spatial birth and death process. If, for
example, we insist that $\gamma(x,m)$ is bounded, then existence (and in
particular non-explosion) is guaranteed by comparison with a pure birth
process. We do not dwell on this, as we shall require more stringent
conditions if we are to pass to the limit as $\theta$ and $N$ tend to
infinity.

It is convenient to characterise the process as a solution to a martingale
problem. We write $C_{b}^{\infty}(\mathbb{R}^{d})$ for the space of bounded
smooth functions on $\mathbb{R}^{d}$, and, where convenient, we write $\langle
f,\eta\rangle=\int_{\mathbb{R}^{d}}f(x)\eta(dx)$.

######  Definition 2.4 (Martingale Problem Characterisation).

For each value of $N$ and $\theta$, and each purely atomic
$\eta_{0}^{N}\in\mathcal{M}_{F}(\mathbb{R}^{d})$ with atoms of mass $1/N$,
$(\eta^{N}_{t})_{t\geq 0}$ is the (scaled) empirical measure of a birth-death
process with càdlàg paths in $\mathcal{M}_{F}(\mathbb{R}^{d})$ for which, for
all $f\in C^{\infty}_{b}(\mathbb{R}^{d})$, writing $q_{\theta}(x,dy)$ for the
Gaussian kernel with mean $x+\vec{b}(x)/\theta$ and covariance
$\mathbf{C}(x)/\theta$,

$\displaystyle M^{N}_{t}(f)$ $\displaystyle:=\langle
f,\eta^{N}_{t}\rangle-\langle f,\eta^{N}_{0}\rangle$ (2.3)
$\displaystyle\qquad{}-\int_{0}^{t}\bigg{\\{}\bigg{\langle}\left(\int\theta\left(f(z)r(z,\eta^{N}_{s})-f(x)r(x,\eta^{N}_{s})\right)q_{\theta}(x,dz)\right)\gamma(x,\eta^{N}_{s}),\eta^{N}_{s}(dx)\bigg{\rangle}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad{}+\bigg{\langle}f(x)F(x,\eta^{N}_{s}),\eta^{N}_{s}(dx)\bigg{\rangle}\bigg{\\}}ds$

is a martingale (with respect to the natural filtration), with angle bracket
process

$\displaystyle\left\langle M^{N}(f)\right\rangle_{t}=$
$\displaystyle\frac{\theta}{N}\int_{0}^{t}\bigg{\\{}\Big{\langle}\gamma(x,\eta^{N}_{s})\int
f^{2}(z)r(z,\eta^{N}_{s})q_{\theta}(x,dz),\eta^{N}_{s}(dx)\Big{\rangle}$ (2.4)
$\displaystyle\qquad\qquad{}+\Big{\langle}\mu_{\theta}(x,\eta^{N}_{s})f^{2}(x),\eta^{N}_{s}(dx)\Big{\rangle}\bigg{\\}}ds.$

The angle bracket process (or, “conditional quadratic variation”) is the
unique previsible process making $(M^{N}(f)_{t})^{2}-\langle
M^{N}(f)\rangle_{t}$ a martingale with respect to the natural filtration. It
differs from the usual quadratic variation (usually denoted $[M^{N}(f)]_{t}$)
because the process has jumps; for the (continuous) limit the two notions will
coincide. The use of angle brackets for both integrals and this process is
unfortunately standard but should not cause confusion, since the angle bracket
process always carries a subscript for time.

The form of (2.3) and (2.4) is explained in Section 4. Note that since
(juvenile) individuals are produced at rate $N\gamma\eta$, but each has mass
$1/N$, these factors of $N$ cancel in (2.3). Under our scaling, $N$ and
$\theta=\theta(N)$ will tend to infinity in such a way that
$\alpha:=\lim_{N\to\infty}\theta(N)/N$ exists and is finite. From the
expression (2.4) it is easy to guess that whether the limiting processes will
be deterministic or stochastic is determined by whether $\alpha$ is zero or
nonzero.

It is convenient to record some notation for the generator of the diffusion
limit of a random walk with jump distribution determined by
$q_{\theta}(x,dy)$.

######  Definition 2.5 (Dispersal generator).

As above, we define the dispersal kernel, $q_{\theta}(x,dy)$, to be the
density of a multivariate Gaussian with mean $\vec{b}(x)/\theta$ and
covariance matrix $\mathbf{C}(x)/\theta$ (although often we omit the
dependence of $\vec{b}$ and $\mathbf{C}$ on $x$). Furthermore, we define for
$f\in C_{b}^{\infty}(\mathbb{R}^{d})$,

$\displaystyle\mathcal{B}f(x)=\sum_{ij}\mathbf{C}(x)_{ij}\partial_{x_{i}}\partial_{x_{j}}f(x)+\sum_{i}\vec{b}(x)_{i}\partial_{x_{i}}f(x)$
(2.5)

and denote the adjoint of $\mathcal{B}$ by

$\displaystyle\mathcal{B}^{*}f(x)$
$\displaystyle=\sum_{ij}\partial_{x_{i}}\partial_{x_{j}}(\mathbf{C}(x)_{ij}f(x))-\sum_{i}\partial_{x_{i}}(f(x)\vec{b}(x)_{i})$
$\displaystyle=\sum_{ij}C_{ij}(x)\partial_{x_{i}}\partial_{x_{j}}f(x)+\sum_{i}\left(\sum_{j}\partial_{x_{j}}C_{ij}(x)-\vec{b}_{i}(x)\right)\partial_{x_{i}}f(x)$
$\displaystyle\qquad{}+\left(\sum_{ij}\partial_{x_{i}}\partial_{x_{j}}C_{ij}(x)-\sum_{i}\partial_{x_{i}}\vec{b}_{i}(x)\right)f(x).$

######  Remark 2.6.

$\mathcal{B}$ is defined so that

$\displaystyle\theta\int\left(f(y)-f(x)\right)q_{\theta}(x,dy)\to\mathcal{B}f(x)\qquad\text{as
}\theta\to\infty.$

######  Remark 2.7.

An equivalent way to describe the model would be to say that when the state of
the population is $\eta$, an individual at $x$ gives birth at rate

$\theta\gamma(x,\eta)\int r(y,\eta)q(x,dy),$

and that offspring disperse according to the kernel

$q_{\theta}^{\mathfrak{m}}(x,\eta,dy):=\frac{r(y,\eta)q_{\theta}(x,dy)}{\int
r(z,\eta)q_{\theta}(x,dz)}.$

Clearly, the random walk driven by this dispersal kernel is biased towards
regions of higher establishment probability. For comparison with future
results, it is interesting to write down the limiting generator:

$\displaystyle\lim_{\theta\to\infty}\theta\int(f(y)-f(x))q_{\theta}^{\mathfrak{m}}(x,\eta,dy)$
$\displaystyle=\frac{\mathcal{B}\left[f(\cdot)r(\cdot,\eta)\right](x)-f(x)\mathcal{B}\left[r(\cdot,\eta)\right](x)}{r(x,\eta^{N})}.$
(2.6)

In the simplest case of unbiased isotropic dispersal (i.e., $\vec{b}=0$ and
$\mathbf{C}=I$), $\mathcal{B}=\Delta$, and so (2.6) is equal to

$\displaystyle\Delta f(x)+2\nabla f(x)\cdot\nabla\log
r(\cdot,\rho_{r}\\!*\\!\eta(\cdot))(x).$

One might guess that the spatial motion described by following the ancestral
lineage of an individual back through time would be described (in the limit)
by the adjoint of this generator. However, we will see in Section 2.2 that
this is not in fact the case.

In order to pass to a scaling limit, we will need to impose some conditions on
the parameters of our model.

######  Assumptions 2.8.

We shall make the following assumptions on the parameters of our model.

Dispersal generator: We assume that

  1. 1.

$\vec{b}(x)$ and $\mathbf{C}(x)$ are $\alpha$-Hölder continuous for some
$\alpha\in(0,1]$ and uniformly bounded in each component, and

  2. 2.

the operator $\mathcal{B}$ is uniformly strictly elliptic, i.e.,
$\inf_{x}\inf_{y:\|y\|=1}\sum_{ij}y_{i}C(x)_{ij}y_{j}>0$.

Reproduction parameters: We assume that

  1. 3.

The function $F(x,m)$ satisfies

    1. (a)

$F(x,m)$ is locally Lipschitz in $m$;

    2. (b)

$F(x,m)$ is uniformly bounded above (but not necessarily below);

    3. (c)

for each fixed $m$, $\sup_{x\in\mathbb{R}^{d}}\sup_{k\leq m}|F(x,k)|<\infty$;

  2. 4.

The functions $r(x,m)$, $\gamma(x,m)$ have bounded first and second
derivatives in both arguments;

  3. 5.

$\gamma(x,m)$ is uniformly bounded;

  4. 6.

For each $f\in C_{b}^{2}(\mathbb{R}^{d})$, there is a $C_{f}$ such that

$|\gamma(x,\eta)\theta\int(r(y,\eta)f(y)-r(x,\eta)f(x))q_{\theta}(x,dy)|\leq
C_{f}(1+|f(x)|)$

for all $x\in\mathbb{R}^{d}$ and $\eta\in\mathcal{M}_{F}(\mathbb{R}^{d})$.
Furthermore, $C_{f}$ only depends on the norm of the first two derivatives of
$f$, i.e.,

$C_{f}=C(\sup_{x}\sup_{\|z\|=1}\max(\sum_{i}z_{i}\partial_{x_{i}}f(x),\sum_{ij}z_{i}z_{j}\partial_{x_{i}x_{j}}f(x))).$

To keep expressions manageable, we shall also assume that

$\mu_{\theta}(x)=r(x,\eta)\gamma(x,\eta)-\frac{1}{\theta}F(x,\eta),$

that is, this expression is non-negative so that there is no need to take the
maximum with zero in (2.2). (This is anyway true for sufficiently large
$\theta$.)

Since we take bounded $f$, for most situations the bound $C_{f}(1+|f(x)|)$ in
Condition 6 above can be safely replaced simply by $C_{f}$; however, this will
be useful in certain situations where we consider a sequence of $f$ with
increasing upper bounds. We now give two concrete situations in which
Condition 6 is satisfied. The proof is in Section 6.1.

######  Lemma 2.9.

Assume that Conditions 2.8 are satisfied, except for Condition 6. If either

  1. 1.

$|\partial_{x_{i}}r(x,\eta)|$ and $|\partial_{x_{i}x_{j}}r(x,\eta)|$ are
uniformly bounded for $x\in\mathbb{R}^{d}$,
$\eta\in\mathcal{M}_{F}(\mathbb{R}^{d})$;

  2. 2.

or, $m^{2}\gamma(x,m)$ is uniformly bounded and there exists $C<\infty$ such
that for $\theta$ sufficiently large, and all $x\in\mathbb{R}^{d}$,
$\eta\in\mathcal{M}_{F}(\mathbb{R}^{d})$,

$\theta\int\big{(}\rho_{r}\\!*\\!\eta(y)-\rho_{r}\\!*\\!\eta(x)\big{)}q_{\theta}(x,dy)\leq
C\rho_{\gamma}\\!*\\!\eta(x),$

and

$\theta\int\big{(}\rho_{r}\\!*\\!\eta(y)-\rho_{r}\\!*\\!\eta(x)\big{)}^{2}q_{\theta}(x,dy)\leq
C(\rho_{\gamma}\\!*\\!\eta(x))^{2},$

then Condition 6 is also satisfied.

The purpose of the conditions that we have placed on the reproduction
parameters is to ensure that the net per capita reproduction rate (before time
scaling) is order $1/\theta$. As remarked above, because of the non-local
reproduction mechanism, it no longer suffices to assume that
$r(x,\eta)\gamma(x,\eta)-\mu_{\theta}(x)$ is of order $1/\theta$. Perhaps the
simplest example in which we can see this is the case where $\gamma\equiv 1$
and $F\equiv 0$, so that $\mu_{\theta}=r$, and $\eta=\delta_{x}$ (i.e., the
population has all individuals at a single location), so that
$\rho_{r}\\!*\\!\eta(y)=\rho_{r}(y)$. In this case, the mean rate of change of
the total population size is
$\int(r(y,\rho_{r}(y))-r(x,\rho_{r}(x)))q_{\theta}(x,dy)$; the first condition
of Lemma 2.9 would ensure this is of order $1/\theta$.

If $r(x,m)$ is independent of $m$, then the conditions are easy to satisfy;
they just require some regularity of $r$ as a function of $x$. Condition 1 of
Lemma 2.9 is also satisfied if for example $\|\nabla\rho_{r}\|\leq C\rho_{r}$
and $m\partial_{m}r(x,m)$, $m^{2}\partial_{mm}r(x,m)$ are bounded. This is the
case, for instance, if $\rho_{r}$ decays exponentially. On the other hand, it
might seem more natural to take $\rho_{r}$ to be a Gaussian density with
parameter $\sigma_{r}$, say. Then, as we check in Lemma B.1, Condition 2 of
Lemma 2.9 is satisfied if $\rho_{\gamma}$ is also Gaussian with parameter
$\sigma_{\gamma}$ and $\sigma_{\gamma}>\sigma_{r}$. For large enough $\theta$,
this condition guarantees that $\sigma_{r}+1/\theta<\sigma_{\gamma}$, so that
the establishment probability of a juvenile is controlled by individuals that
are already ‘felt’ by the fecundity-regulating kernel $\rho_{\gamma}$ at the
location of their parent.

###  2.1 Scaling limits of the population process

Our main results depend on two dichotomies: Is the limiting process
deterministic or a (generalized) superprocess? And, are interactions pointwise
in the limit or nonlocal? See Figure 1 for snapshots of the population from
direct simulation of the process using SLiM (Haller and Messer, 2019)
illustrating this first dichotomy. Below we have results for deterministic
limits with pointwise and nonlocal interactions, and for superprocess limits
with nonlocal interactions.

Figure 1:  Snapshots of two simulations, with small $\alpha=\theta/N$ (left)
and large $\alpha=\theta/N$ (right). Simulations are run with a Fisher-KPP-
like parameterization: birth and establishment are constant, while death
increases linearly with density, at slope $1/\theta$. Left: $\alpha=0.1$.
Right: $\alpha=10$. Other parameters were the same: dispersal and interactions
distance were set to 1, and the equilibrium density is 10 individuals per unit
area.

##### Scaling limits with nonlocal interactions:

Recall that the process $(\eta_{t}^{N})_{t\geq 0}$ takes its values in the
space ${\mathcal{D}}_{[0,\infty)}(\mathcal{M}_{F}(\mathbb{R}^{d}))$ of càdlàg
paths on $\mathcal{M}_{F}(\mathbb{R}^{d})$. We endow
$\mathcal{M}_{F}(\mathbb{R}^{d})$ with the topology of weak convergence and
${\mathcal{D}}_{[0,\infty)}(\mathcal{M}_{F}(\mathbb{R}^{d}))$ with the
Skorohod topology. A sequence of processes taking values in
${\mathcal{D}}_{[0,\infty)}(\mathcal{M}_{F}(\mathbb{R}^{d}))$ is said to be
tight if the corresponding sequence of distributions is tight, i.e., that any
infinite subsequence has a weakly convergent subsubsequence. Our first main
result establishes tightness of our rescaled population processes in the case
in which interactions remain nonlocal under the scaling, and characterises
limit points as solutions to a martingale problem.

######  Theorem 2.10.

Let $(\eta^{N}_{t})_{t\geq 0}$ be as defined in Definition 2.4 and assume that
as $N\to\infty$, $\theta(N)\to\infty$ in such a way that
$\theta(N)/N\to\alpha$. (However, the kernels $\rho_{r}$, $\rho_{\gamma}$, and
$\rho_{F}$ remain fixed.) Suppose that Assumptions 2.8 hold and, further, that
$\\{\eta_{0}^{N}\\}_{N\geq 1}$ is a sequence of purely atomic measures, with
$\eta_{0}^{N}$ comprised of atoms of mass $1/N$, which is tight in
$\mathcal{M}_{F}(\mathbb{R}^{d})$. Also assume there exists a nonnegative
$f_{0}\in C(\mathbb{R}^{d})$ with uniformly bounded first and second
derivatives (i.e., with
$\sup_{x}\sup_{\|z\|=1}\sum_{i}\partial_{x_{i}}f_{0}(x)z_{i}$ and
$\sup_{x}\sup_{\|z\|=1}\sum_{ij}\partial_{x_{i}x_{j}}f_{0}(x)z_{i}z_{j}$ both
finite) and $f_{0}(x)\to\infty$ as $|x|\to\infty$ for which $\langle
f_{0}(x),\eta_{0}^{N}(dx)\rangle<C<\infty$ for some $C$ independent of $N$.
Then the sequence of processes $(\eta^{N}_{t})_{t\geq 0}$ is tight, and for
any limit point $(\eta_{t})_{t\geq 0}$, for every $f\in
C^{\infty}_{b}(\mathbb{R}^{d})$,

$\displaystyle\begin{split}M_{t}(f)&:=\langle f(x),\eta_{t}(dx)\rangle-\langle
f(x),\eta_{0}(dx)\rangle\\\
&\qquad-\int_{0}^{t}\big{\langle}\gamma(x,\eta_{s})\mathcal{B}\left(f(\cdot)r(\cdot,\eta_{s})\right)(x)+f(x)F(x,\eta_{s}),\eta_{s}(dx)\big{\rangle}ds\end{split}$
(2.7)

is a martingale (with respect to the natural filtration), with angle bracket
process

$\displaystyle\left\langle
M(f)\right\rangle_{t}=\alpha\int_{0}^{t}\big{\langle}2\gamma\left(x,\eta_{s}\right)r\left(x,\eta_{s}\right)f^{2}(x),\eta_{s}(dx)\big{\rangle}ds.$
(2.8)

If $\alpha=0$ the limit is deterministic.

Recall when interpreting (2.7) that, for instance,
$r(x,\eta_{s})=r(x,\rho_{r}\\!*\\!\eta_{s}(x))$, and so
$\mathcal{B}(fr)(x)=\mathcal{B}(f(\cdot)r(\cdot,\rho_{r}\\!*\\!\eta_{s}(\cdot)))(x)$.
The proof of this theorem appears in Section 6.2.

Theorem 2.10 provides tightness of the rescaled processes. If the limit points
are unique, then this is enough to guarantee convergence.

######  Corollary 2.11.

Under the assumptions of Theorem 2.10, if the martingale problem defined by
equations (2.7) and (2.8) has a unique solution, then $(\eta^{N}_{t})_{t\geq
0}$ converges weakly to that solution as $N\to\infty$.

When $\alpha>0$, the limit points can be thought of as interacting
superprocesses. For example, when $r$ and $\gamma$ are constant, and
$F(x,\eta_{s})=1-\rho_{F}\\!*\\!\eta_{s}(x)$, we recover a superprocess with
nonlinear death rates corresponding to logistic growth (Etheridge, 2004) that
is a continuous limit of the Bolker-Pacala model (Bolker and Pacala, 1997,
1999). We are not aware of a general result to determine when we will have
uniqueness of solutions to the martingale problem of Theorem 2.10 when
$\alpha>0$. However, the Dawson–Girsanov transform tells us that we have
uniqueness in this special case of the superprocess with nonlinear death
rates, and Perkins stochastic calculus (and its adaptation to a lookdown
setting) provides uniqueness for cases with interactions in the dispersal
mechanism of the superprocess. We refer to Dawson (1993); Perkins (1992), and
Donnelly and Kurtz (1999) for approaches to showing that these sorts of
martingale problems are well-posed.

For the deterministic case of $\alpha=0$, the limiting process is a weak
solution to a nonlocal PDE. We next describe some situations in which more is
known about uniqueness and whether the solution is close to the corresponding
local PDE. First, recall the following notion of solution to a PDE.

######  Definition 2.12 (Weak solutions).

We say that $(\eta_{t})_{t\geq 0}$, with
$\eta_{t}\in\mathcal{M}_{F}(\mathbb{R}^{d})$, is a _weak solution_ to the PDE

$\partial_{t}\varphi=r\mathcal{B}^{*}(\gamma\varphi)+\varphi F$ (2.9)

(where $r$, $\gamma$ and $F$ can all be functions of $\varphi$) if, for all
$f\in C_{b}^{\infty}(\mathbb{R}^{d})$,

$\frac{d}{dt}\langle
f,\eta_{t}\rangle=\langle\gamma\mathcal{B}(rf)+fF,\eta_{t}\rangle.$

The notation $\varphi$ is meant to be suggestive of a density, and recall that
equation (2.9) has made dependencies on $x$ and $\varphi$ implicit; written
out more explicitly, (2.9) is

$\partial_{t}\varphi_{t}(x)=r\left(x,\rho_{r}\\!*\\!\varphi_{t}(x)\right)\mathcal{B}^{*}\left[\varphi_{t}(\cdot)\gamma\big{(}x,\rho_{\gamma}\\!*\\!\varphi_{t}(\cdot)\big{)}\right](x)+\varphi_{t}(x)F\left(x,\rho_{F}\\!*\\!\varphi_{t}(x)\right).$

Because Theorem 2.10 only tells us about weak convergence, in the case
$\alpha=0$ we can only deduce that any limit point $\eta_{t}$ is a weak
solution to this nonlocal PDE.

Specialising the results of Kurtz and Xiong (1999) to the deterministic
setting provides general conditions under which we have existence and
uniqueness of solutions to (2.9) which have an $L^{2}$-density with respect to
Lebesgue measure. Recall that the Wasserstein metric, defined by

$\rho(\nu_{1},\nu_{2})=\sup\Big{\\{}\Big{|}\int fd\nu_{1}-\int
fd\nu_{2}\Big{|}:\sup_{x}|f(x)|\leq 1,|f(x)-f(y)|\leq\|x-y\|\Big{\\}},$

determines the topology of weak convergence on
$\mathcal{M}_{F}(\mathbb{R}^{d})$. We write
$r(x,\eta)\gamma(x,\eta)\mathbf{C}(x)=J(x,\eta)J(x,\eta)^{T}$, and
$\beta(x,\eta)=r(x,\eta)\gamma(x,\eta)\big{(}\vec{b}(x)+2\mathbf{C}(x)\nabla\log
r(x,\eta)\big{)}$ (quantities that will appear in Proposition 5.6). If $J$,
$\beta$, and $F$ are bounded and Lipschitz in the sense that

$|J(x_{1},\nu_{1})-J(x_{2},\nu_{2})|,|\beta(x_{1},\nu_{1})-\beta(x_{2},\nu_{2})|,|F(x_{1},\nu_{1})-F(x_{2},\nu_{2})|\leq
C(\|x_{1}-x_{2}\|+\rho(\nu_{1},\nu_{2}))$ (2.10)

for some $C>0$, the methods of Kurtz and Xiong (1999) show that if the initial
condition $\eta_{0}$ for our population process has an $L^{2}$ density, then
so does $\eta_{t}$ for $t>0$. Although the necessary estimates (for which we
refer to the original paper) are highly nontrivial, the idea of the proof is
simple. Take a solution to the equation and use it to calculate the
coefficients $r$, $\gamma$ and $F$ that depend on local population density.
Then $\eta$ solves the linear equation obtained by regarding those values of
$r$, $\gamma$ and $F$ as given. It remains to prove that the solution to the
linear equation has a density which is achieved by obtaining $L^{2}$ bounds on
its convolution with the heat semigroup at time $\delta$ and letting
$\delta\to 0$. We also have the following uniqueness result.

######  Theorem 2.13 (Special case of Kurtz and Xiong (1999), Theorem 3.5).

Suppose $J$, $\beta$, and $F$ are bounded and Lipschitz in the sense of
(2.10). If $\eta_{0}$ has an $L^{2}({\mathbb{R}}^{d})$-density, then there
exists a unique $L^{2}({\mathbb{R}}^{d})$-valued solution of (2.9) in the
sense of Definition 2.12.

######  Remark 2.14.

Kurtz and Xiong (1999) considers an infinite system of stochastic differential
equations for the locations and weights of a collection of particles that
interact through their weighted empirical measure, which is shown to be the
unique solution to a stochastic PDE. As we shall see through our lookdown
representation in Section 5, the solution to our deterministic equation can be
seen as the empirical measure of a countable number of particles (all with the
same weight) which, in the notation above, evolve according to

$X(t)=X(0)+\int_{0}^{t}\beta\big{(}X(s),\eta_{s}\big{)}ds+\int_{0}^{t}J\big{(}X(s),\eta_{s}\big{)}dW(s)$

(with an independent Brownian motion $W$ for each particle).

##### Two-step convergence to PDE:

Although the coefficients at $x$ in (2.9) are nonlocal, we can choose our
kernels $\rho_{\gamma}$, $\rho_{r}$, and $\rho_{F}$ in such a way that they
depend only on the population in a region close to $x$, and so we expect that
under rather general conditions solutions of the nonlocal PDE will be close to
the corresponding classical PDE. The following propositions provide two
concrete situations in which this is true. In the first, the PDE is a
reaction-diffusion equation, and in the proof in Section 6.3.1 we borrow an
idea from Penington (2017) to express the solutions to both the nonlocal
equation and the classical PDE through a Feynman-Kac formula.

######  Proposition 2.15.

Let $\rho^{\epsilon}_{F}(x)=\rho_{F}\big{(}x/\epsilon)/\epsilon^{d}$. Assume
$\varphi_{0}\in L^{2}(\mathbb{R})$ is a positive, uniformly Lipschitz, and
uniformly bounded function. Suppose that $\varphi^{\epsilon}\in
L^{2}(\mathbb{R}^{d})$ is a weak solution to the equation

$\partial_{t}\varphi^{\epsilon}=\mathcal{B}^{*}\varphi^{\epsilon}+\varphi^{\epsilon}F(\rho_{F}^{\epsilon}*\varphi^{\epsilon}),\qquad
x\in\mathbb{R}^{d},\,t>0,$ (2.11)

with initial condition $\varphi_{0}(\cdot)$, and that $\varphi$ is a weak
solution to the equation

$\partial_{t}\varphi=\mathcal{B}^{*}\varphi+\varphi F(\varphi),\qquad
x\in\mathbb{R}^{d},\,t>0,$ (2.12)

also with initial condition $\varphi_{0}(\cdot)$. Suppose further that $F$ is
a Lipschitz function which is bounded above, and that $\vec{b}(x)$ and
$\mathbf{C}(x)$, the drift and covariance matrix of $\mathcal{B}$, satisfy the
conditions of Assumptions 2.8. Then, for all $T>0$ there exists a constant
$K=K(T,\|\varphi_{0}\|_{\infty})<\infty$ and a function $\delta(\epsilon)$
(dependent on $\rho_{F}$) with $\delta(\epsilon)\to 0$ as $\epsilon\to 0$,
such that, for all $0\leq t\leq T$, and $\epsilon$ small enough,

$\|\varphi_{t}(\cdot)-\varphi_{t}^{\epsilon}(\cdot)\|_{\infty}\leq
K\delta(\epsilon).$

In particular, as $\epsilon\to 0$, we have that $\varphi^{\epsilon}$ converges
uniformly in compact intervals of time to $\varphi$.

######  Remark 2.16.

Note that Theorem 2.13 guarantees uniqueness of solutions to equation (2.11).

Our second example in which we know solutions to the nonlocal PDE converge to
solutions of the local PDE as interaction distances go to zero is a nonlocal
version of a porous medium equation with logistic growth. That is, we consider
non-negative solutions to the equation

$\partial_{t}\psi^{\epsilon}=\Delta\left(\psi^{\epsilon}\,\rho^{\epsilon}_{\gamma}*\psi^{\epsilon}\right)+\psi^{\epsilon}\left(1-\rho^{\epsilon}_{\gamma}*\psi^{\epsilon}\right).$
(2.13)

The case without the reaction term (and with $\mathbb{R}^{d}$ replaced by a
torus) is considered by Lions and Mas-Gallic (2001) who use it as a basis for
a particle method for numerical solution of the porous medium equation. Of
course this does not quite fit into our framework, since in the notation of
our population models this would necessitate
$\gamma(x,m)=\rho_{\epsilon}\\!*\\!m$ which is not bounded. However, this can
be overcome by an additional layer of approximation (c.f. our numerical
experiments of Section 3.1) and we do not allow this to detain us here.
Existence and uniqueness of solutions to (2.13) can be obtained using the
approach of Lions and Mas-Gallic (2001), so we should like to prove that as
$\epsilon\to 0$ we have convergence to the solution to the porous medium
equation with logistic growth:

$\partial_{t}\psi=\Delta\left(\psi^{2}\right)+\psi\left(1-\psi\right).$ (2.14)

######  Notation 2.17.

We use $\rightharpoonup$ to denote weak convergence in the sense of analysts;
that is, $\psi^{\epsilon}\rightharpoonup\psi$ in $L^{1}$ means
$\int\psi^{\epsilon}vdx\rightarrow\int\psi vdx$ for all $v\in L^{\infty}$.

We write $L_{t}^{2}(H^{1})$ for functions for which the $H^{1}$ norm in space
is in $L^{2}$ with respect to time, i.e.

$\int_{0}^{T}\int\left\\{\psi_{t}(x)^{2}+\left\|\nabla\psi_{t}(x)\right\|^{2}\right\\}dxdt<\infty,$

and $C_{t}(L^{1})$ will denote functions for which the $L^{1}$ norm in space
is continuous in time.

######  Proposition 2.18.

Suppose that we can write $\rho_{\gamma}=\zeta*\check{\zeta}$, where
$\check{\zeta}(x)=\zeta(-x)$ and $\zeta\in\mathcal{S}(\mathbb{R}^{d})$ (the
Schwartz space of rapidly decreasing functions). Furthermore, suppose that
$\psi_{0}^{\epsilon}\geq 0$ is such that there exist $\lambda\in(0,1)$ and
$C\in(0,\infty)$ (independent of $\epsilon$) such that

$\int\exp(\lambda\|x\|)\psi_{0}^{\epsilon}(x)dx<C,\quad\mbox{ and
}\sup_{\epsilon}\int\psi_{0}^{\epsilon}|\log\psi_{0}^{\epsilon}|dx<\infty,$

with $\psi_{0}^{\epsilon}\rightharpoonup\psi_{0}$ as $\epsilon\to 0$. Then
writing $\psi^{\epsilon}$ for the solution to (2.13) on
$[0,T]\times\mathbb{R}^{d}$ with initial condition $\psi_{0}^{\epsilon}$,
$\psi^{\epsilon}\rightharpoonup\psi$ as $\epsilon\to 0$ where $\psi\in
L_{t}^{2}(H^{1})\cap C_{t}(L^{1})$, $\int\psi|\log\psi|dx<\infty$, and $\psi$
solves (2.14) on $[0,T]\times\mathbb{R}^{d}$.

The example that we have in mind for the kernel $\rho_{\gamma}$ is a Gaussian
kernel. For the proof, see Section 6.3.2.

######  Remark 2.19.

Although it seems hard to formulate an all-encompassing result, Propositions
2.15 and 2.18 are by no means exhaustive. When the scaling limit is
deterministic, one can expect analogous results under rather general
conditions. However, when the limit points are stochastic, they resemble
“nonlinear superprocesses” and so one cannot expect a density with respect to
Lebesgue measure in $d\geq 2$. It is then not reasonable to expect to be able
to make sense of the limit if we scale the kernels in this way. Moreover, in
one dimension, where the classical superprocess does have a density with
respect to Lebesgue measure, the form of (2.7) suggests that even if one can
remove the local averaging from $\gamma$, it will be necessary to retain
averaging of $r$ in order to obtain a well-defined limit.

##### One-step convergence to PDE:

Theorem 2.10, combined with Proposition 2.15 or 2.18 implies that we can take
the limit $N\to\infty$ followed by the limit $\epsilon\to 0$ to obtain
solutions to the PDE (2.12). However, it is of substantial interest to know
whether we can take those two limits simultaneously. The general case seems
difficult, but we prove such “diagonal” convergence in the following
situation. The proof is provided in Section 7.

######  Theorem 2.20 (Convergence to a PDE).

Let $(\eta^{N}_{t})_{t\geq 0}$ be as defined in Definition 2.4 with
$r(x,m)\equiv 1\equiv\gamma(r,m)$, $F(x,m)\equiv F(m)$, $\rho^{\epsilon}_{F}$
a symmetric Gaussian density with variance parameter $\epsilon^{2}$, and
$\mathcal{B}=\Delta/2$. Further suppose that $F(m)$ is a polynomial with
$F(m)\mathbf{1}_{m\geq 0}$ bounded above. Assume that $\langle
1,\eta_{0}^{N}\rangle$ is uniformly bounded, and that for all
$x\in\mathbb{R}^{d}$ and $k\in\mathbb{N}$,

$\limsup_{\epsilon\to
0}\mathbb{E}\big{[}\rho_{F}^{\epsilon}*\eta_{0}(x)^{k}\big{]}<\infty,$

and

$\limsup_{\epsilon\to
0}\int\mathbb{E}\big{[}\rho_{F}^{\epsilon}*\eta_{0}(x)^{k}\big{]}dx<\infty.$

Finally assume that $N\to\infty$, $\theta\to\infty$ and $\epsilon\to 0$ in
such a way that

$\frac{1}{\theta\epsilon^{2}}+\frac{\theta}{N\epsilon^{d}}\to 0.$ (2.15)

Then the sequence of
${\mathcal{D}}_{[0,\infty)}(\mathcal{M}_{F}(\mathbb{R}^{d}))$-valued
stochastic processes
$\big{(}\rho_{F}^{\epsilon}*\eta^{N}_{t}(x)dx\big{)}_{t\geq 0}$ converges
weakly to a measure-valued process with a density $\varphi(t,x)$ that solves

$\displaystyle\partial_{t}\varphi(t,x)=\frac{1}{2}\Delta\varphi(t,x)+\varphi(t,x)F(\varphi(t,x)).$
(2.16)

######  Remark 2.21.

In fact, our proof goes through without significant change under the
conditions that $F(m)\mathbf{1}_{m\geq 0}$ is bounded above (but not
necessarily below), and that for all $m,n\in[0,\infty)$

$|F(m)|\leq\sum_{j=1}^{k}a_{j}m^{j},\quad\mbox{ and
}|F(n)-F(m)|\leq|n-m|\sum_{j=1}^{k^{\prime}}b_{j}\Big{(}n^{j}+m^{j}\Big{)},$

for some non-negative constants $\\{a_{j}\\}_{j=0}^{k}$,
$\\{b_{j}\\}_{j=0}^{k^{\prime}}$. We take $F$ to be polynomial to somewhat
simplify notation in the proof.

###  2.2 Ancestral lineages in the scaling limit

Now that we have established what we can say about how population density
changes with time, we turn to results on ancestral lineages, i.e., how
genealogical ancestry can be traced back across the landscape. Informally, a
_lineage_ $(L_{t}^{N})_{t\geq 0}$, begun at a spatial location $L_{0}^{N}=x$
where there is a focal individual in the present day can be obtained by, for
each time $t$, setting $L_{t}^{N}$ to be the spatial location of the
individual alive at time $t$ before the present from whom the focal individual
is descended. Since in our model individuals have only one parent, this is
unambiguous. Although we did not explicitly retain such information, it is
clear that for finite $N$, since individuals are born one at a time, one could
construct the lineage $(L_{t}^{N})_{t=0}^{T}$ given the history of the
population $(\eta^{N}_{t})_{t=0}^{T}$, for each starting location to which
$\eta^{N}_{T}$ assigns positive mass. It is less clear, however, how to
formally retain such information when we pass to the scaling limit.

The _lookdown construction_ in Section 5 will enable us to recover information
about ancestry in the infinite population limit. Roughly speaking, each
particle is assigned a unique “level” from $[0,\infty)$ that functions as a
label and thus allows reconstruction of lineages. The key to the approach is
that levels are assigned in such a way as to be exchangeable, so that sampling
a finite number, $k$ say, of individuals from a given region is equivalent to
looking at the individuals in that region with the $k$ lowest levels.
Moreover, as we pass to the infinite population limit, the collection of
(individual, level) pairs converges, as we show in Theorem 5.4. See Etheridge
and Kurtz (2019) for an introduction to these ideas. In particular, even in
the infinite population limit, we can sample an individual from a region (it
will be the individual in that region with the lowest level) and trace its
line of descent. This will allow us to calculate, for each $x$ and
$y\in\mathbb{R}^{d}$, the proportion of the population at location $x$ in the
present day population that are descended from a parent who was at location
$y$ at time $t$ in the past. To make sense of this in our framework, in
Section 8.2, we justify a weak reformulation of this idea.

We are interested in two questions. First, when is the motion of an ancestral
lineage, given complete knowledge of the population process, a well-defined
process? In other words, is knowledge of the process $(\eta_{t})_{t=0}^{T}$
that records numbers of individuals but not their ancestry sufficient to
define the distribution of $(L_{t})_{t=0}^{T}$? Second, does the process have
a tractable description?

We focus on the simplest situation, that in which the population process is
deterministic. However, the results here apply when the population process
solves either a nonlocal or a classical PDE. There will be no coalescence of
ancestral lineages in the deterministic limit, but understanding motion of
single lineages is useful in practice, and our results can be seen as a first
step towards understanding genealogies for high population densities.

Proofs of results in this section are found in Section 8.

######  Definition 2.22 (Ancestral lineage).

Let $(\varphi_{t}(x))_{0\leq t\leq T}$ denote the density of the scaling limit
of our population model, solving (2.9), and let $y$ be a point with
$\varphi_{T}(y)>0$. We define $(L_{s})_{s=0}^{T}$, the ancestral lineage of an
individual sampled from the population at $y$ at time $T$, by setting
$L_{0}=y$ and $L_{s}$ to be the position of the unique ancestor of that
individual at time $T-s$. We define $(Q_{s})_{s\geq 0}$ to be the time
inhomogeneous semigroup satisfying

$\displaystyle Q_{s}f(y):=\mathbb{E}_{y}[f(L_{s})].$

Our next result identifies the ancestral lineage as a diffusion by
characterizing its generator.

######  Theorem 2.23.

For $\varphi:\mathbb{R}^{d}\to\mathbb{R}$, define

$\displaystyle\mathcal{L}_{\varphi}f$
$\displaystyle=\frac{r}{\varphi}\left[\mathcal{B}^{*}(\gamma\varphi
f)-f\mathcal{B}^{*}(\gamma\varphi)\right]$ (2.17)
$\displaystyle=r\gamma\left[\sum_{ij}\mathbf{C}_{ij}\partial_{x_{i}x_{j}}f+\sum_{j}\vec{m}_{j}\partial_{x_{j}}f\right],$
(2.18)

where $\vec{m}$ is the vector

$\vec{m}_{j}=2\sum_{i}C_{ij}\partial_{x_{i}}\log(\gamma\varphi)+2\sum_{i}\partial_{x_{i}}C_{ij}-\vec{b}_{j}.$

Then the generator of the semigroup $Q_{s}$ of Definition 2.22 is given by
$\partial_{s}Q_{s}f(y)=\mathcal{L}_{\varphi_{T-s}}Q_{s}f(y).$

######  Remark 2.24.

As usual, to make the generator readable, we’ve written it in concise
notation, omitting the dependencies on location and population density, which
itself changes with time. When interpreting this, remember that everything
depends on location and density at that location and time – for instance,
“$r$” is actually $r(x,\varphi(x))$ (in the classical case), or
$r(x,\rho_{r}\\!*\\!\eta(x))$ (in the nonlocal case).

Moreover, we haven’t proved any regularity of the population density process
$\varphi$, so, as written, the generator (2.17) may not make sense. Instead,
it should be interpreted in a weak sense which is made precise in Section 8.2.

######  Corollary 2.25.

In addition to the assumptions of Theorem 2.23, if the covariance of the
dispersal process is isotropic (i.e., $\mathbf{C}=\sigma^{2}I$), then

$\mathcal{L}_{\varphi}f=r\gamma\left(\sigma^{2}\Delta
f+\left(2\sigma^{2}\nabla\log(\gamma\varphi)-\vec{b}\right)\cdot\nabla
f\right).$ (2.19)

(However, $\vec{b}$ can still depend on location.)

In other words, the lineage behaves as a diffusion driven by Brownian motion
run at speed $\sigma^{2}$ multiplied by the local per-capita production of
mature offspring ($r\gamma$) in a potential tilted by migration bias
($(\varphi_{s}\gamma)^{2}\exp(-\vec{b}\cdot x/\sigma^{2})$, whose gradient
appears in the drift term of the generator). In particular, lineages are drawn
to regions of high fecundity (production of juveniles), but their speed is
determined by the rate of production of mature offspring. This can be compared
to Remark 2.7.

######  Corollary 2.26.

In addition to the assumptions of Corollary 2.25, if the population process is
stationary (so $\varphi_{t}\equiv\varphi$), and $\vec{b}(x)=\nabla h(x)$ for
some function $h$, then $Y$ is reversible with respect to

$\displaystyle\pi(x)=\frac{\gamma}{r}\varphi(x)^{2}e^{-h(x)/\sigma^{2}}.$
(2.20)

Long-term fitness of an individual is proportional to the fraction of lineages
from the distant future that pass through the individual, and hence the total
long-term fitness at a location is proportional to the stationary distribution
of $Y$ there, if it exists. Therefore, if $\pi$ is integrable then the per-
capita long-term fitness of an individual at $x$ is proportional to
$\pi(x)/\varphi(x)$.

######  Corollary 2.27.

In addition to the assumptions of Corollary 2.25, suppose that the population
process is described by a travelling wave with velocity $\mathfrak{c}$, i.e.,
the population has density $\varphi(t,x)=w(x-t\mathfrak{c})$ where $w$ solves

$\displaystyle r\mathcal{B}^{*}(\gamma w)+wF+\mathfrak{c}\cdot\nabla
w=r\sigma^{2}\Delta(\gamma w)-\vec{b}\cdot\nabla(\gamma
w)+\mathfrak{c}\cdot\nabla w=0.$

Then the semigroup $Q_{s}$ of the motion of a lineage in the frame that is
moving at speed $\mathfrak{c}$ is time-homogeneous with generator

$\displaystyle\mathcal{L}f$ $\displaystyle=\sigma^{2}r\gamma\left(\Delta
f+2\nabla\log(\gamma w)\cdot\nabla
f\right)+(\mathfrak{c}-r\gamma\vec{b})\cdot\nabla f.$ (2.21)

##  3 Examples and applications

We now discuss some consequences of these results.

###  3.1 Beyond linear diffusion

Equation (2.9) is a nonlocal version of a reaction-diffusion equation; the
diffusion is nonlinear if $\gamma$ depends on population density: in other
words, if the diffusivity of the population depends on the population density.
Passing to the classical limit, we recover equations like (2.14). Such
equations are widely used in a number of contexts in biology in which motility
within a population varies with population density. For example, density
dependent dispersal is a common feature in spatial models in ecology,
eukaryotic cell biology, and avascular tumour growth; see Sherratt (2010) and
references therein for further discussion. In particular, it has been
suggested as a model for the expansion of a certain type of bacteria on a thin
layer of agar in a Petri dish (Cohen et al., 1999). We shall pay particular
attention to the case in which the equation can be thought of as modelling the
density of an expanding population. We focus on the monostable reaction of
(2.14).

Comparing with (2.9), we see that to set up a limit in which the population
density $\varphi$ follows the porous medium equation with logistic growth of
(2.14), we need $r=1$, $\gamma=\varphi$, and $F=1-\varphi$. Consulting
equation (2.2), this implies that
$\mu_{\theta}=\max\big{(}0,(1+1/\theta)\varphi-1/\theta\big{)}$. In other
words, establishment is certain and birth rates increase linearly with
population density, but to compensate, death rates increase slightly faster
(also linearly). Alert readers will notice that the condition from Assumptions
2.8 that $\gamma(x,m)$ be uniformly bounded is violated. This can be corrected
by use of a cut-off, and in fact the downwards drift provided by the logistic
control of the population size prevents $m$ from getting too big. In practice
the simulations shown in Figure 2 take discrete time steps of length $dt$
(with $dt$ suitably small), and have each individual reproduce and die with
probabilities, respectively,

$\displaystyle p_{\text{birth}}(m)=\left(1-e^{-mdt}\right)\qquad
p_{\text{death}}(m)=\left(1-e^{-(m(1+1/\theta)-1/\theta)dt}\right),$

where $m$ is the local density at their location. This makes
$\gamma(x,m)=p_{\text{birth}}(m)/dt\approx m$ and

$\displaystyle F(x,m)=\theta(\gamma(x,m)-\mu(x,m))/dt\approx 1-m.$

Birth and death rates are equal at density $m=1$, corresponding to an unscaled
density of $N$ individuals per unit area.

Figure 2:  Simulated populations under a porous medium equation with logistic
growth (2.14) in $d=1$, $\theta/N$ small on the left; large on the right.
Values of $\theta$ in top and bottom figures are 1 and 100, respectively, and
both have $N$ set so that the density is roughly 100 individuals per unit of
habitat (as displayed on the vertical axis). See text for details of the
simulations.

In one dimension, equation (2.14) has an explicit travelling wave solution

$\displaystyle w^{P}(t,x):=\left(1-e^{\frac{1}{2}(x-x_{0}-t)}\right)_{+}.$
(3.1)

Notice that the wave profile has a sharp boundary at $x=x_{0}+t$. There are
also travelling wave solutions with $c>1$ (Gilding and Kersner, 2005), which
lack this property. However, for initial conditions that decay sufficiently
rapidly at infinity, such as one might use in modelling a population invading
new territory, the solution converges to (3.1) (Kamin and Rosenau, 2004). In
Figure 2 we show simulations of the individual based model described above,
which display travelling wave solutions qualitatively similar to solutions of
(2.14), with better agreement for smaller $\theta/N$ (but in both cases, $N$
is reasonably large).

###  3.2 Ancestry in different types of travelling waves

Although it remains challenging to establish the distribution of genealogical
trees relating individuals sampled from our population model, as described in
the introduction, we can gain some insight by investigating the motion of a
single ancestral lineage. Here we do that in the context of a one-dimensional
population expanding into new territory as a travelling wave. We focus on
three cases in which we have explicit information about the shape of the
travelling wave profile: the Fisher-KPP equation, a special case of the Allen-
Cahn equation with a bistable nonlinearity, and the porous media equation with
logistic growth, equation (2.14). We work here in one dimension, and take
$\sigma^{2}=1$ and $\vec{b}=0$.

##### Fisher-KPP equation:

Consider the classical Fisher-KPP equation,

$\displaystyle\partial_{t}\varphi=\partial_{xx}\varphi+\varphi(1-\varphi).$
(3.2)

Even though we do not have an explicit formula for the wave shape in this
case, our methods provide information about ancestral lineages. The equation
has non-negative travelling wave solutions of speed $c$ for all $c\geq 2$,
but, started from any compact perturbation of a Heaviside function, the
solution will converge to the profile $w^{F}$ with the minimal wavespeed,
$c=2$ (Kolomogorov et al., 1937; Fife and McLeod, 1977; Bramson, 1983). No
matter what initial condition, for any $t>0$ the support of the solution will
be the whole real line. In this case, we must have $r=\gamma=1$, and
$F(x,m)=1-m$ so $\mu_{\theta}(x,m)=1+(m-1)/\theta$. By Corollary 2.27, the
generator of the motion of an ancestral lineage is

$\mathcal{L}_{F}f=\partial_{xx}f+2\frac{\partial_{x}w^{F}}{w^{F}}\partial_{x}f+2\partial_{x}f.$
(3.3)

Near the tip of the wave (for $x$ large), $w^{F}(x)\sim e^{-x}$, so (3.3)
implies that the motion of a lineage is close to unbiased Brownian motion. On
the other hand, in the “bulk”, a lineage behaves approximately as Brownian
motion with drift at rate two to the right. This implies that ancestral
lineages are pushed into the tip of the wave, and there is no stationary
distribution, so that long-term dynamics of genetic inheritance depend on the
part of the wave not well-approximated by a smooth profile, in agreement with
the previous results referred to in the Introduction.

##### Allen-Cahn equation:

Now take the Allen-Cahn equation:

$\displaystyle\partial_{t}\varphi=\partial_{xx}\varphi+\varphi(1-\varphi)(2\varphi-1+s),$
(3.4)

for a given $s\in(0,2)$. Once again we have taken $r=\gamma=1$, but now the
reaction term $F(x,m)=(1-m)(2m-1+s)$ is bistable. This equation can be used to
model the motion of so-called hybrid zones in population genetics; see, for
example, Barton (1979), Gooding (2018), and Etheridge et al. (2022). This
equation has an explicit travelling wave solution with speed $s$ and shape

$w^{A}(x)=(1+e^{x})^{-1},$

i.e., $\phi_{t}(x)=w^{A}(x-st)$ solves (3.4). Substituting $w^{A}$ in place of
$w^{F}$ in (3.3), we find that the generator of an ancestral lineage relative
to the wavefront is now,

$\displaystyle\mathcal{L}_{A}f$
$\displaystyle=\partial_{xx}f+2\frac{\partial_{x}w^{A}}{w^{A}}\partial_{x}f+s\partial_{x}f$
$\displaystyle=\partial_{xx}f-2\frac{e^{x}}{1+e^{x}}\partial_{x}f+s\partial_{x}f,$

so lineages in the tip are pushed leftwards into the bulk of the wave at a
rate $s-2e^{x}/(1+e^{x})$. The density of the speed measure for this diffusion
is

$m_{A}(x)\propto e^{sx}(1+e^{x})^{-2},$

which is integrable, and so determines the unique stationary distribution.
Thus the position of the ancestral lineage relative to the wavefront will
converge to a stationary distribution which is maximised away from the extreme
tip of the wave. This is consistent with Etheridge and Penington (2022), who
consider an analogous stochastic population model, although the stronger
result there (that the genealogy of a sample from behind the wavefront is
approximately a Kingman coalecsent) requires the stronger condition $s<1$.

##### Porous Medium equation with logistic growth:

Finally, consider equation (2.14). Setting $x_{0}=0$ (for definiteness) and
substituting the form of $w^{P}$ from equation (3.1) into Corollary 2.27, with
$c=1$, $\gamma(x,m)=m$, $r(x,w)=1$, and $F(x,m)=(1-m)$, the generator of the
diffusion governing the position of the ancestral lineage relative to the
wavefront is, for $x<0$,

$\displaystyle\mathcal{L}_{P}f$
$\displaystyle=w^{P}\left(\partial_{xx}f+2\frac{\partial_{x}((w^{P})^{2})}{(w^{P})^{2}}\partial_{x}f\right)+\partial_{x}f$
$\displaystyle=\left(1-e^{\frac{1}{2}x}\right)\partial_{xx}f-2e^{\frac{1}{2}x}\partial_{x}f+\partial_{x}f.$

The speed measure corresponding to this diffusion has density

$\displaystyle m_{P}(\xi)$
$\displaystyle\propto\frac{1}{2(1-e^{\xi/2})}\exp\left(\int_{\eta}^{\xi}\left\\{1-\frac{e^{x/2}}{1-e^{x/2}}\right\\}dx\right)$
$\displaystyle\propto e^{\xi}\left(1-e^{\xi/2}\right),\quad\text{for }\xi<0$

and $m_{P}(\xi)=0$ for $\xi\geq 0$, which is integrable and so when suitably
normalised gives the unique stationary distribution. Notice that even though
we have the same reaction term as in the Fisher-KPP equation, with this form
of nonlinear diffusion, at stationarity the lineage will typically be
significantly behind the front, suggesting a different genealogy.

It is interesting to compare the stationary distribution we have obtained here
to the expression that we’d get by setting $\vec{b}=-\mathfrak{c}$ and using
Corollary 2.26, i.e., by giving each offspring a mean displacement that
offsets the motion of the wave. In the Fisher-KPP and Allen-Cahn cases above,
we get the same expressions, but this is only the case because
$r\equiv\gamma\equiv 1$ in both. For the PME with $\vec{b}=1$ we have
$\mathcal{B}^{*}=\Delta+\nabla$ and so the equation solved by the population
density is

$\partial_{t}\varphi=\Delta\varphi+2\varphi\nabla\varphi+\varphi(1-\varphi),$

which has a traveling wave solution of the same shape but moving at half the
speed, $\varphi(x,t)=w^{P}(x-t/2)$, and a stationary distribution of the
lineage relative to the wavefront of $\pi(\xi)\propto
e^{3\xi/2}(1-e^{\xi/2})^{3}$.

###  3.3 Clumping from nonlocal interactions

Simulating these processes and exploring parameter space, one sooner or later
comes upon a strange observation: with certain parameter combinations, the
population spontaneously forms a regular grid of stable, more or less discrete
patches, separated by areas with nearly no individuals, as shown in Figure 3.
The phenomenon is discussed in Section 16.10 of Haller and Messer (2022), and
has been described in similar models, e.g., by Britton (1990); Sasaki (1997);
Hernández-García and López (2004); Young et al. (2001), and Berestycki et al.
(2009). For example, if the density-dependent effects of individuals extend
farther (but not too much farther) than the typical dispersal distance, then
depending on the interaction kernel new offspring landing between two clumps
can effectively find themselves in competition with _both_ neighbouring
clumps, while individuals within a clump compete with only one.

More mathematically, consider the case in which $\mathcal{B}=\sigma^{2}\Delta$
and all parameters are spatially homogeneous, so that
$r(x,\eta)=r(\rho_{r}\\!*\\!\eta(x))$, and similarly for $\gamma$ and $F$. If
$\varphi_{0}$ is such that $F(\varphi_{0})=0$ and $F^{\prime}(\varphi_{0})<0$,
then the constant solution $\varphi\equiv\varphi_{0}$ is a nontrivial
equilibrium of (1.1). However, this constant solution may not be unique, it
may be unstable, and a stable solution may have oscillations on a scale
determined by the interaction distance.

To understand the stability of the constant solution
$\varphi\equiv\varphi_{0}$, we linearise (1.1) around $\varphi_{0}$: let
$\varphi_{t}(x)=\varphi_{0}+\psi_{t}(x)$, and (informally) $r(x)\approx
r(\varphi_{0})+r^{\prime}(\varphi_{0})\rho_{r}\\!*\\!\psi(x)$. Writing
$r_{0}=r(\varphi_{0})$ and $r^{\prime}_{0}=r^{\prime}(\varphi_{0})$, with
analogous expressions for $\gamma$ and $F$,

$\displaystyle\partial_{t}\psi$
$\displaystyle\approx\sigma^{2}\varphi_{0}r_{0}\gamma_{0}^{\prime}\Delta\rho_{\gamma}\\!*\\!\psi+\sigma^{2}r_{0}\gamma_{0}\Delta\psi+\varphi_{0}F^{\prime}_{0}\rho_{F}\\!*\\!\psi.$

Letting $\widehat{f}(u)=\int e^{2\pi iux}f(x)dx$ denote the Fourier transform,

$\partial_{t}\widehat{\psi}(u)\approx\left\\{-u^{2}\sigma^{2}\varphi_{0}r_{0}\gamma_{0}^{\prime}\widehat{\rho}_{\gamma}(u)-u^{2}\sigma^{2}r_{0}\gamma_{0}+\varphi_{0}F^{\prime}_{0}\widehat{\rho}_{F}(u)\right\\}\widehat{\psi}(u).$
(3.5)

In the simplest case, in which $\gamma$ is constant, so
$\gamma_{0}^{\prime}=0$, this reduces to

$\partial_{t}\widehat{\psi}(u)\approx\left(-u^{2}\sigma^{2}r_{0}\gamma_{0}+\varphi_{0}F^{\prime}_{0}\widehat{\rho}_{F}(u)\right)\widehat{\psi}(u).$
(3.6)

If we take $\rho_{F}=p_{\epsilon^{2}}$, then
$\widehat{\rho}_{F}(u)=\exp(-2\pi\epsilon^{2}u^{2})$ and (recalling that
$F_{0}^{\prime}<0$) the term in brackets is always negative, and we recover
the well-known fact that in this case the constant solution is stable. If, on
the other hand, $\widehat{\rho}_{F}$ changes sign, there may be values of $u$
for which the corresponding quantity is positive. For example, if $d=1$ and
$\rho_{F}(x)={\mathbf{1}}_{[-\epsilon,\epsilon]}(x)/2\epsilon$, then
$\widehat{\rho}_{F}(u)=\sin(2\pi\epsilon u)/(2\pi\epsilon u)$, which is
negative for $u\in(1/(2\epsilon),1/\epsilon)$ (and periodically repeating
intervals). Setting $v=\epsilon u$, the bracketed term on the right hand side
of (3.6) becomes

$\varphi_{0}F_{0}^{\prime}\frac{1}{2\pi v}\sin(2\pi
v)-\frac{\sigma^{2}}{\epsilon^{2}}v^{2}r_{0}\gamma_{0},$

and we see that if $\sigma^{2}/\epsilon^{2}$ is sufficiently small, there are
values of $v$ for which this is positive. In other words, in keeping with our
heuristic above, if dispersal is sufficiently short range relative to the
range over which individuals interact, there are unstable frequencies that
scale with the interaction distance $\epsilon$. In two dimensions, replacing
the indicator of an interval by that of a ball of radius $\epsilon$, a similar
analysis applies, except that the sine function is replaced by a Bessel
function.

Now suppose that $\gamma$ is not constant. Then, from (3.5), if we take
$\rho_{\gamma}=\rho_{F}=p_{\epsilon}^{2}$,

$\partial_{t}\widehat{\psi}(u)=e^{-2\pi^{2}\epsilon^{2}u^{2}}\left\\{-\sigma^{2}\varphi_{0}r_{0}\gamma_{0}^{\prime}u^{2}-\sigma^{2}r_{0}\gamma_{0}u^{2}e^{2\pi^{2}\epsilon^{2}u^{2}}+\varphi_{0}F_{0}^{\prime}\right\\}\widehat{\psi}(u).$

If we make the (reasonable) assumption that $\gamma_{0}^{\prime}<0$, then we
see that even when the Fourier transform of $\rho$ does not change sign, there
may be parameter values for which the constant solution is unstable. As
before, we set $v=\epsilon u$. The term in brackets becomes

$\frac{\sigma^{2}}{\epsilon^{2}}v^{2}r_{0}\left(-\varphi_{0}\gamma_{0}^{\prime}-\gamma_{0}e^{2\pi^{2}v^{2}}\right)+\varphi_{0}F_{0}^{\prime},$

and, provided $-\varphi_{0}\gamma_{0}^{\prime}/\gamma_{0}>1$, for sufficiently
small $v$ the term in round brackets is positive. We now see that if
$\sigma^{2}/\epsilon^{2}$ is sufficiently large, the equilibrium state
$\varphi\equiv\varphi_{0}$ is unstable. As before, the unstable frequencies
will scale with $\epsilon$ and for given $F$, $r$ and $\gamma$, whether or not
such unstable frequencies exist will be determined by
$\sigma^{2}/\epsilon^{2}$, but in this case of Gaussian kernels, it is
interaction distance being sufficiently small relative to dispersal that will
lead to instability.

Figure 3:  Left: A snapshot of individual locations in a two-dimensional
simulation in which the constant density is unstable and a stable, periodic
pattern forms. Right: Population density in an expanding wave in a one-
dimensional simulation forming a periodic pattern; each panel shows the
wavefront in three periods of time; within each period of time the wavefront
at earlier times is shown in blue and later times in pink. In both cases,
$\gamma(m)=3/(1+m)$, $\mu\equiv 0.3$, and $r\equiv 1$; dispersal is Gaussian
with $\sigma=0.2$ and density is measured with $\rho_{\gamma}(x)=p_{9}(x)$,
i.e., using a Gaussian kernel with standard deviation 3.

###  3.4 Lineage motion is not uniquely determined by population density

It is natural for applications to wonder about identifiability: when can the
observed quantities like population density or certain summaries of lineage
movement uniquely determine the underlying demographic parameters? Consider a
deterministic, continuous population generated by parameters $\gamma$, $r$,
and $F$, with $\vec{b}=0$ and $\mathbf{C}=I$. Suppose it has a stationary
profile $w(x)$, that must satisfy

$r\Delta(\gamma w)+Fw=0.$

It is easy to see that $w$ does not uniquely specify $\gamma$, $F$, and $r$:
let $\lambda(x)$ be a smooth, nonnegative function on $\mathbb{R}^{d}$, and
let $\widetilde{r}(x,m)=\lambda(x)r(x,m)$ and
$\widetilde{F}(x,m)=\lambda(x)F(x,m)$ (and, let $\widetilde{\gamma}=\gamma$).
Since $\mu=r\gamma-F/\theta$, this corresponds to multiplying both
establishment probabilities and death rates by $\lambda$. Then the population
with parameters $\widetilde{\gamma}$, $\widetilde{r}$, and $\widetilde{F}$ has
the same stationary profile(s) as the original population.

Can these two situations be distinguished from summaries of lineage movement?
The first has lineage generator

$f\mapsto\mathcal{L}f=r\gamma\left(\Delta f+2\nabla\log(\gamma w)\cdot\nabla
f\right),$

while the second has lineage generator $f\mapsto\lambda(x)\mathcal{L}f(x)$. In
other words, although the stationary profile of the population is unchanged
when we scale local establishment and death by $\lambda$, the motion of
lineages is sped up locally by $\lambda$. This corresponds to making areas
with $\lambda>1$ more “sink-like” and $\lambda<1$ more “source-like”: if
$\lambda(x)>1$, then at $x$ both the death rate and probability of
establishment of new individuals are higher. As a result, lineages in the
second model spend more time in areas with $\lambda<1$, i.e., those areas have
higher reproductive value, something that is, in principle, discernible from
genetic data (because, for instance, making reproductive value less evenly
distributed reduces long-term genetic diversity).

##  4 Heuristics

In this section we perform some preliminary calculations and use them to
provide heuristic arguments for our main results, to build intuition before
the formal proofs.

###  4.1 The population density

We reiterate that in our prelimiting model, the population is represented by a
point measure $\eta^{N}$ in which each individual is assigned a mass $1/N$. We
use the term “population density” for this process, as it is supposed to
measure population size relative to a nominal occupancy of $N$ individuals per
unit area, but it is not absolutely continuous with respect to Lebesgue
measure.

We write $\mathcal{P}^{N}$ for the generator of the scaled population process
$\eta^{N}$ of Definition 2.4 acting on test functions of the form $G(\langle
f,\eta\rangle)$, where $f\geq 0$ is smooth and bounded on $\mathbb{R}^{d}$ and
$G\in C^{\infty}([0,\infty))$. Recall that $\theta=\theta(N)\to\infty$ as
$N\to\infty$ in such a way that $\theta(N)/N\to\alpha$.

A Taylor expansion allows us to write

$\mathcal{P}^{N}G(\langle f,\eta\rangle)=G^{\prime}(\langle
f,\eta\rangle)\lim_{\delta t\downarrow 0}\frac{1}{\delta
t}\mathbb{E}\left[\left.\langle f,\eta_{\delta t}\rangle-\langle
f,\eta\rangle\right|\eta_{0}=\eta\right]\\\
\qquad{}+\frac{1}{2}G^{\prime\prime}(\langle f,\eta\rangle)\lim_{\delta
t\downarrow 0}\frac{1}{\delta t}\mathbb{E}\left[\left.\big{(}\langle
f,\eta_{\delta t}\rangle-\langle
f,\eta\rangle\big{)}^{2}\right|\eta_{0}=\eta\right]+\epsilon_{N}(f,G,\eta),$
(4.1)

where the terms that make up $\epsilon_{N}(f,G,\eta)$ will be negligible in
our scaling limit (at least if $G^{\prime\prime\prime}<\infty$).

#### Mean measure

Recall that in our parameterization only death rates $\mu_{\theta}$ and the
dispersal kernel $q_{\theta}$ depend on $\theta$. For a suitable test function
$f$, we find

$\begin{split}\mathcal{P}^{N}\langle f,\eta\rangle&=\lim_{\delta t\downarrow
0}\frac{1}{\delta t}\mathbb{E}\left[\left.\langle f,\eta_{\delta
t}\rangle-\langle f,\eta\rangle\right|\eta_{0}=\eta\right]\\\ &=\theta\int\int
f(z)r(z,\eta)q_{\theta}(x,dz)\gamma(x,\eta)\eta(dx)-\theta\int
f(x)\mu_{\theta}(x,\eta)\eta(dx).\end{split}$ (4.2)

The first term is the increment in $\langle f,\eta\rangle$ resulting from a
birth event (recalling that we don’t kill the parent) integrated against the
rate of such events, and the second reflects death events. The factor of
$\theta$ appears from the time rescaling. In both terms, the rate of events
has a factor of $N$ (because events happen at a rate proportional to the
number of individuals, whereas $\eta$ has mass $1/N$ for each individual)
which is offset by the fact that the birth or loss of a single individual at
the point $y$, say, changes $\langle f,\eta\rangle$ by $f(y)/N$.

We use the fact that $\int q_{\theta}(x,dz)=1$ to rewrite (4.2) as

$\begin{split}\int\left(\int\theta\left(f(z)r(z,\eta)-f(x)r(x,\eta)\right)q_{\theta}(x,dz)\right)\gamma(x,\eta)\eta(dx)\\\
+\int
f(x)\theta\Big{(}r(x,\eta)\gamma(x,\eta)-\mu_{\theta}(x,\eta)\Big{)}\eta(dx).\end{split}$
(4.3)

We have defined $\mu_{\theta}$ so that the second term is simple:

$\displaystyle\theta\Big{(}r(x,\eta)\gamma(x,\eta)-\mu_{\theta}(x,\eta)\Big{)}=F(x,\eta).$

Furthermore, recall from Remark 2.6 that

$\displaystyle\int\theta\Big{(}r(z,\eta)f(z)-r(x,\eta)f(x)\Big{)}q_{\theta}(x,dz)\qquad\stackrel{{\scriptstyle\theta\to\infty}}{{\longrightarrow}}\qquad\mathcal{B}\big{(}r(\cdot,\eta)f(\cdot)\big{)}(x).$
(4.4)

In particular, if dispersal is determined by a standard multivariate Gaussian
with mean zero and covariance $\sigma^{2}I/\theta$, then
$\mathcal{B}=\sigma^{2}\Delta$, where $\Delta$ denotes the Laplacian.

In summary, equation (4.3) converges to

$\int\gamma(x,\eta)\mathcal{B}\big{(}f(\cdot)r(\cdot,\eta)\big{)}(x)\eta(dx)+\int
f(x)F(x,\eta)\eta(dx),$ (4.5)

which explains the form of the martingale of Theorem 2.10.

#### Quadratic variation

We now look at the second order term in (4.1), which will converge to the
quadratic variation of the limiting process. An individual at location $x$
gives birth to a surviving offspring at $y$ at rate

$\gamma(x,\eta)r(y,\eta)q_{\theta}(x,dy),$

and since this increments $\langle f,\eta\rangle$ by $f(y)/N$, the
contribution to the quadratic variation from birth events, which occur at rate
$\theta$ per individual (so, rate $N\theta|\eta|$ overall), is

$\int
N\theta\gamma(x,\eta)\int\frac{1}{N^{2}}f^{2}(y)r(y,\eta)q_{\theta}(x,dy)\eta(dx).$

Similarly, the increment in $\langle f,\eta\rangle$ resulting from the death
of an individual at $x$ is $-f(x)/N$, and so combining with the above, the
second order term in the generator takes the form

$\displaystyle G^{\prime\prime}(\langle
f,\eta\rangle)\frac{1}{2}N\theta\left\\{\int\gamma(x,\eta)\int\frac{1}{N^{2}}f^{2}(y)r(y,\eta)q_{\theta}(x,dy)\eta(dx)+\int\mu_{\theta}(x,\eta)\frac{1}{N^{2}}f^{2}(x)\eta(dx)\right\\}$
$\displaystyle\qquad=\frac{1}{2}G^{\prime\prime}(\langle
f,\eta\rangle)\frac{\theta}{N}\int\left\\{\gamma(x,\eta)\int
f^{2}(y)r(y,\eta)q_{\theta}(x,dy)+f^{2}(x)\mu_{\theta}(x,\eta)\right\\}\eta(dx).$

Since $\int f^{2}(y)r(y,\eta)q_{\theta}(x,dy)\to f^{2}(x)r(x,\eta)$ and
$r\gamma+\mu_{\theta}=2r\gamma-F/\theta\to 2r\gamma$ as $\theta\to\infty$,
this converges to

$\displaystyle\frac{\alpha}{2}G^{\prime\prime}(\langle
f,\eta\rangle)\big{\langle}2r(x,\eta)\gamma(x,\eta),\eta(dx)\big{\rangle}.$

An entirely analogous argument shows that if $G^{\prime\prime\prime}$ is
bounded, then the term $\epsilon_{\theta,N}(f,G,\eta)$ in (4.1) will be
$\mathcal{O}(\theta/N^{2})$.

If we hold $\rho_{\gamma}$, $\rho_{r}$, $\rho_{F}$ fixed, then by taking
$\theta/N\rightarrow 0$, the second order term in the generator will vanish
and we expect a deterministic limit, for which $\partial_{t}\langle
f,\eta_{t}\rangle$ is equal to (4.5). In other words, the limit is a weak
solution to the deterministic equation

$\partial_{t}\varphi_{t}(x)=r(x,\varphi_{t})\mathcal{B}\big{(}\gamma(\cdot,\varphi_{t})\varphi_{t}(\cdot)\big{)}(x)+F(x,\varphi_{t})\varphi_{t}(x)$
(4.6)

in the sense of Definition 2.12, where $\varphi_{t}$ is the density of
$\eta_{t}$, if it has a density. On the other hand, if $N=\alpha\theta$ for
some $\alpha>0$, the second order term remains, and we expect a “generalised
superprocess” limit. The limiting quadratic variation is exactly as seen in
Theorem 2.10.

##### One-step convergence:

In order to pass directly to a classical PDE limit in Theorem 2.20 we impose
the stronger condition that $\theta/N\epsilon^{d}\to 0$ and also require that
$\theta\epsilon^{2}\to\infty$. Recall that in this case, we take
$\rho_{F}^{\epsilon}$ to be a symmetric Gaussian density with variance
$\epsilon^{2}$. The condition $\theta\epsilon^{2}\to\infty$ ensures that
$\epsilon^{2}$ is large enough relative to $1/\theta$ that the regularity
gained by smoothing our population density by convolution with
$\rho_{\epsilon}$ is preserved under the dynamics dictated by $q_{\theta}$. To
understand the first condition, note that we are aiming to obtain a
deterministic expression for the limiting population density. It is helpful to
think about a classical Wright-Fisher model (with no spatial structure and
just two types, say). We know then that if the timescale $\theta$ is on the
same order as population size $N$, we see stochastic fluctuations in the
frequencies of the two types in the limit as $N\to\infty$; to obtain a
deterministic limit, we look over timescales that are short relative to
population size. In our setting, the total population size is replaced by the
local population size, as measured by convolution with $\rho_{\epsilon}$,
which we expect to be of order $N\epsilon^{d}$, and so in order to ensure a
deterministic limit we take $\theta/(N\epsilon^{d})\to 0$.

###  4.2 Motion of ancestral lineages

Although our proof of Theorem 2.23 uses an explicit representation in terms of
the lookdown process, the result can be understood through informal
calculations. Suppose that we have traced a lineage back to an individual at
location $y$ at time $t$. Looking further back through time, at the time of
the birth of that individual, the lineage will jump to the location of the
parent of the individual. Now, the rate at which new individuals are born to
parents at $x$ and establish at $y$ is

$\theta
N\eta^{N}_{t}(dx)\gamma(x,\eta^{N}_{t})q_{\theta}(x,dy)r(y,\eta^{N}_{t}).$

Suppose that $\eta^{N}$ did have a density (in the prelimit it does not), say
$\eta^{N}_{t}(dx)=\varphi^{N}_{t}(x)dx$. Informally, since the number of
individuals near $y$ is $N\varphi^{N}_{t}(y)dy$, the probability that a
randomly chosen individual near $y$ is a new offspring from a parent at $x$ in
$[t,t+dt)$ is

$\frac{\theta\varphi^{N}_{t}(x)\gamma(x,\eta^{N}_{t})r(y,\eta^{N}_{t})}{\varphi^{N}_{t}(y)}\frac{q_{\theta}(x,dy)}{dy}dxdt.$
(4.7)

Leaving aside questions of whether a lineage can be treated as a randomly
chosen individual, we define a continuous-time jump process whose transition
rates, conditional on $(\varphi^{N}_{t})_{t=0}^{T}$, are given by (4.7).
Because we are tracing the lineage backwards in time we make the substitution
$s=T-t$ and write $(L^{N}_{s})_{s=0}^{T}$ for the location of a lineage that
moves according to these jump rates. Then, abusing notation to write
$q_{\theta}(x,y)$ for the density of $q_{\theta}(x,dy)$,

$\displaystyle\begin{split}&\mathbb{E}[f(L^{N}_{s+ds})-f(y)\;|\;L^{N}_{s}=y]\\\
&\qquad=ds\,\theta\int\left(f(x)-f(y)\right)\frac{\varphi^{N}_{T-s}(x)\gamma(x,\eta^{N}_{T-s})r(y,\eta^{N}_{T-s})}{\varphi^{N}_{T-s}(y)}q_{\theta}(x,y)dx.\end{split}$
(4.8)

(Note that this integral is with respect to $x$.) Referring back to Remark
2.6, a quick calculation shows that as $N\to\infty$,

$\displaystyle\theta\int\big{(}f(x)-f(y)\big{)}g(x)q_{\theta}(x,y)dx$
$\displaystyle\qquad=\theta\int\big{\\{}(f(x)g(x)-f(y)g(y))-f(y)(g(x)-g(y))\big{\\}}q_{\theta}(x,y)dx$
$\displaystyle\qquad\to\mathcal{B}^{*}(fg)(y)-f(y)\mathcal{B}^{*}g(y).$

Applying this to (4.8) with $g=\varphi_{T-s}\gamma$, this suggests that the
generator of the limiting process is

$\displaystyle\mathcal{L}_{s}f$
$\displaystyle=\frac{r}{\varphi_{T-s}}\left\\{\mathcal{B}^{*}(\gamma\varphi_{T-s}f)-f\mathcal{B}^{*}(\gamma\varphi_{T-s})\right\\}.$
(4.9)

This agrees with Theorem 2.23.

##  5 The lookdown process

Our characterisation of the motion of lines of descent (from which we
establish that of ancestral lineages) when we pass to the scaling limit in our
model will be justified via a lookdown construction. In this section we
present such a construction for the general population model of Definition
2.4. It will be in the spirit of Kurtz and Rodrigues (2011). The general set-
up is as follows. Each individual will be labelled with a “level”, a number in
$[0,N]$. We will still encode the process embellished by these levels as a
point measure: if the $i^{\mathrm{th}}$ individual’s spatial location is
$x_{i}$ and level is $u_{i}$, then we will write

$\xi^{N}=\sum_{i}\delta_{x_{i},u_{i}},$

which is a measure on $\mathbb{R}^{d}\times[0,N]$. Note that each individual
contributes mass 1 to the measure, not $1/N$ as above. If we assign mass $1/N$
to each individual and ignore the levels we will recover our population model.
Moreover, at any time, the levels of individuals in a given spatial region
will be exchangeable and conditionally uniform on $[0,N]$: in particular,
choosing the $k$ individuals with the lowest levels in that region is
equivalent to taking a uniform random sample of size $k$ from the population
in the region. However, this exchangeability is only as regards the _past_ :
an individual’s level encodes information about their future reproductive
output, since individuals with lower levels tend to live longer, and have more
offspring. For more explanation of the set-up and how this is possible, see
Kurtz and Rodrigues (2011) and Etheridge and Kurtz (2019) (and note that our
$N$ corresponds to the $\lambda$ of those papers). The power of this approach
is that we can pass to a limit under the same scalings as described in Theorem
2.10, and the limiting “spatial-level” process will still be a point measure,
and so we explicitly retain the notion of individuals and lineages in the
infinite-population limit.

###  5.1 Lookdown representation of the model of Definition 2.4

For the remainder of this section, when there is no risk of ambiguity we shall
suppress the superscript $N$ on the processes $\eta$ and $\xi$.

In this subsection, we’ll define the process $(\xi_{t})_{t\geq 0}$ in terms of
the dynamics of labelled particles, and write down its generator. The dynamics
depend on the spatial locations of particles, and in this section $\eta_{t}$
is the corresponding spatial measure, i.e.,

$\eta_{t}(\cdot)=\frac{1}{N}\xi_{t}(\cdot\times[0,N]).$

A nontrivial consequence of the way we define $\xi_{t}$ will be that this
process has the same distribution as the process $(\eta_{t})_{t\geq 0}$ of
Definition 2.4. (This provides our justification for using the same notation
for both.)

Following Etheridge and Kurtz (2019), we build the generator step by step from
its component parts. Suppose that the initial population is composed of $O(N)$
particles with levels uniformly distributed on $[0,N]$, and that the current
state of the population is $\xi$, with spatial projection $\eta$.

An individual at spatial location $x$ with level $u$ produces one juvenile
offspring at rate

$2\theta\left(1-\frac{u}{N}\right)\gamma(x,\eta),$

which disperses to a location relative to $x$ drawn from the kernel
$q_{\theta}(x,\cdot)$. Averaging over the uniform distribution of the level
$u$, we recover the birth rate $\theta\gamma(x,\eta)$. This juvenile – suppose
its location is $y$ – either survives, with probability $r(y,\eta)$, or
immediately dies. (As before, “maturity” is instantaneous.) If it survives, a
new level $u_{1}$ is sampled independently and uniformly from $[u,N]$, and the
parent and the offspring are assigned in random order to the levels
$\\{u,u_{1}\\}$. This random assignment of levels to parent and offspring will
ensure that assignment of individuals to levels remains exchangeable.

Evidently this mechanism increases the proportion of individuals with higher
levels. To restore the property that the distribution of levels is
conditionally uniform given $\eta$, we impose that the level $v$ of an
individual at location $x$ evolves according to the differential equation

$\dot{v}=-\theta\frac{v}{N}\left(N-v\right)\gamma(x,\eta)\int_{\mathbb{R}^{d}}r(y,\eta)q_{\theta}(x,dy).$

Since $v\in[0,N]$, this moves levels down; see Etheridge and Kurtz (2019),
Section 3.4 for a detailed explanation.

Levels never cross below 0, while particles whose levels move above $N$ are
regarded as dead (and are removed from the population). Therefore, in order to
incorporate death, the level of the individual at location $x$ with level $u$
moves upwards at an additional rate $\theta\mu_{\theta}(x,\eta)u$. Since
levels are uniform, it is easy to check that if $\mu_{\theta}$ were constant,
this would imply an exponential lifetime for each individual; see Etheridge
and Kurtz (2019), Section 3.1 for more general justification.

Putting these together, the level $u$ of an individual at $x$ evolves
according to:

$\dot{u}=-\theta\frac{u}{N}\left(N-u\right)\gamma(x,\eta)\int_{\mathbb{R}^{d}}r(y,\eta)q_{\theta}(x,dy)+\theta\mu_{\theta}(x,\eta)u.$
(5.1)

We shall write

$\displaystyle
b_{\theta}(x,\eta):=\theta\left(\gamma(x,\eta)\int_{\mathbb{R}^{d}}r(y,\eta)q_{\theta}(x,dy)-\mu_{\theta}(x,\eta)\right),$

which captures the local net difference between reproduction and death, and

$\displaystyle
c_{\theta}(x,\eta):=\frac{\theta}{N}\gamma(x,\eta)\int_{\mathbb{R}^{d}}r(y,\eta)q_{\theta}(x,dy),$
(5.2)

which captures the local rate of production of successful offspring. Recall
from equation (2.2) that
$F(x,\eta)=\theta(r(x,\eta)\gamma(x,\eta)-\mu_{\theta}(x,\eta))$, and so

$b_{\theta}(x,\eta)=\theta\gamma(x,\eta)\int_{\mathbb{R}^{d}}\left(r(y,\eta)-r(x,\eta)\right)q_{\theta}(x,dy)+F(x,\eta).$
(5.3)

Under Assumptions 2.8, as $\theta\to\infty$ and $c_{\theta}(x,\eta)$ will tend
to $\alpha\gamma(x,\eta)r(x,\eta)$

$b_{\theta}(x,\eta)\to\gamma(x,\eta)\mathcal{B}r(x,\eta)+F(x,\eta).$ (5.4)

We can then rewrite the differential equation governing the dynamics of the
level of each individual as

$\displaystyle\dot{u}$
$\displaystyle=\theta\gamma(x,\eta)\int_{\mathbb{R}^{d}}r(y,\eta)q_{\theta}(x,dy)\left\\{-\frac{u}{N}\left(N-u\right)+u\right\\}-b_{\theta}(x,\eta)u$
$\displaystyle=c_{\theta}(x,\eta)u^{2}-b_{\theta}(x,\eta)u.$ (5.5)

Now, we can write down the generator for $(\xi_{t})_{t\geq 0}$, the lookdown
process. In what follows, we will write sums (and, products) over
“$(x,u)\in\xi$” to mean a sum over the (location, level) pairs of each
individual in the population. Test functions for $\xi$ will take the form

$f(\xi)=\prod_{(x,u)\in\xi}g(x,u)=\exp\left(\int\log g(x,u)\xi(dx,du)\right),$
(5.6)

where $g(x,u)$ is differentiable in $u$ and smooth in $x$. We will also assume
that $0\leq g(x,u)\leq 1$ for all $u\in[0,N]$, and $g(x,u)\equiv 1$ for $u\geq
N$. In the expressions that follow, we shall often see one or more factor of
$1/g(x,u)$; it should be understood that if $g(x,u)=0$, then it simply cancels
the corresponding factor in $f(\xi)$.

First consider the terms in the generator that come from birth events. When a
birth successfully establishes, a new level is generated above the parent’s
level, and this new level is assigned to either the offspring or the parent.
Since the probability of each is 1/2, the contribution of birth to the
generator maps $f(\xi)$ to

$\displaystyle
f(\xi)\sum_{(x,u)\in\xi}2\frac{\theta}{N}\gamma(x,\eta)\int_{u}^{N}\int_{\mathbb{R}^{d}}\left(\frac{1}{2}\bigg{\\{}g(y,u_{1})+\frac{g(y,u)g(x,u_{1})}{g(x,u)}\bigg{\\}}-1\right)r(y,\eta)q_{\theta}(x,dy)du_{1}$
(5.7)
$\displaystyle\begin{split}&=f(\xi)\sum_{(x,u)\in\xi}2\gamma(x,\eta)\bigg{\\{}\frac{1}{2N}\int_{u}^{N}g(x,u_{1})du_{1}\frac{\theta\int_{\mathbb{R}^{d}}(g(y,u)-g(x,u))r(y,\eta)q_{\theta}(x,dy)}{g(x,u)}\\\
&\qquad\qquad\qquad{}+\frac{\theta}{N}\int_{u}^{N}\int_{\mathbb{R}^{d}}\left(\frac{g(y,u_{1})+g(x,u_{1})}{2}-1\right)r(y,\eta)q_{\theta}(x,dy)du_{1}\bigg{\\}}.\end{split}$
(5.8)

In (5.7), $u_{1}$ is the new level and $y$ is the offspring’s location, and so
the two terms in the integral correspond to the two situations: in the first,
we have added an individual at $(y,u_{1})$, while in the second, we replace an
individual at $(x,u)$ by one at $(x,u_{1})$ and another at $(y,u)$. We’ve
rewritten it in the form (LABEL:eqn:birth_generator) because each of the two
pieces naturally converges to a separate term in the limit.

The remaining term in the generator is due to the motion of particles’ levels.
Reading off from (5.5), it takes the form

$\displaystyle
f(\xi)\sum_{(x,u)\in\xi}\left(c_{\theta}(x,\eta)u^{2}-b_{\theta}(x,\eta)u\right)\frac{\partial_{u}g(x,u)}{g(x,u)}.$
(5.9)

We can now define the spatial-level process explicitly as a solution to a
martingale problem, whose generator is just the sum of
(LABEL:eqn:birth_generator) and (5.9). We need some notation. Write
${\mathcal{C}}={\mathcal{C}}(\mathbb{R}^{d}\times[0,\infty))$ for the counting
measures on $\mathbb{R}^{d}\times[0,\infty)$ and ${\mathcal{C}}_{N}$ for the
subset consisting of counting measures on $\mathbb{R}^{d}\times[0,N]$.

######  Definition 5.1 (Martingale Problem Characterisation).

For given positive values of $N$ and $\theta$, define the generator $A^{N}$ by

$\begin{split}&A^{N}f(\xi)\\\
&\quad=f(\xi)\,\sum_{(x,u)\in\xi}\,2\gamma(x,\eta)\Bigg{\\{}\frac{1}{2N}\int_{u}^{N}g(x,u_{1})du_{1}\frac{\theta\int_{\mathbb{R}^{d}}(g(y,u)-g(x,u))r(y,\eta)q_{\theta}(x,dy)}{g(x,u)}\\\
&\qquad\qquad\qquad\qquad\qquad\qquad\qquad{}+\frac{\theta}{N}\int_{u}^{N}\int_{\mathbb{R}^{d}}\left(\frac{g(y,u_{1})+g(x,u_{1})}{2}-1\right)r(y,\eta)q_{\theta}(x,dy)du_{1}\Bigg{\\}}\\\
&\qquad\qquad{}+f(\xi)\sum_{(x,u)\in\xi}\,\left(c_{\theta}(x,\eta)u^{2}-b_{\theta}(x,\eta)u\right)\frac{\partial_{u}g(x,u)}{g(x,u)},\end{split}$
(5.10)

where $f(\xi)=\prod_{(x,u)\in\xi}g(x,u)$ is as defined in (5.6), and
$\eta(\cdot)=\xi(\cdot\times[0,N])/N$ as before. Given
$\xi_{0}\in{\mathcal{C}}_{N}$, we say that a
${\mathcal{D}}_{[0,\infty)}({\mathcal{C}}_{N})$-valued process
$(\xi_{t})_{t\geq 0}$ is a solution to the $(A^{N},\xi_{0})$ martingale
problem if $f(\xi_{t})-f(\xi_{0})-\int_{0}^{t}A^{N}f(\xi_{s})ds$ is a
martingale (with respect to the natural filtration) for all test functions $f$
as defined above.

The martingale problem for finite $N$ has a unique solution. Next we state the
limiting martingale problem, for which we do not necessarily have uniqueness.
As before, the parameter $\alpha$ will correspond to
$\lim_{N\to\infty}\theta(N)/N$. Whereas for finite $N$, conditional on the
population process $\eta_{t}^{N}$, the levels of particles are independent and
uniformly distributed on $[0,N]$, in the infinite population limit,
conditional on $\eta_{t}$, the process $\xi_{t}$ is Poisson distributed on
$\mathbb{R}^{d}\times[0,\infty)$ with mean measure $\eta_{t}\times\lambda$,
where $\lambda$ is Lebesgue measure.

######  Definition 5.2 (Martingale Problem Characterisation, scaling limit).

Fix $\alpha\in[0,\infty)$, and define test functions $f$ by
$f(\xi)=\prod_{(x,u)\in\xi}g(x,u)$ with $g$ differentiable in $u$, smooth in
$x$, satisfying $0\leq g(x,u)\leq 1$ and such that there exists a $u_{0}$ with
$g(x,u)=1$ for all $u>u_{0}$. Then, define the operator $A$ on such test
functions by

$\displaystyle\begin{split}Af(\xi)&=f(\xi)\sum_{(x,u)\in\xi}\gamma(x,\eta)\frac{\mathcal{B}(g(\cdot,u)r(\cdot,\eta))(x)-g(x,u)\mathcal{B}r(x,\eta)}{g(x,u)}\\\
&\qquad{}+f(\xi)\sum_{(x,u)\in\xi}2\alpha\gamma(x,\eta)r(x,\eta)\int_{u}^{\infty}(g(x,u_{1})-1)du_{1}\\\
&\qquad{}+f(\xi)\sum_{(x,u)\in\xi}\left(\alpha\gamma(x,\eta)r(x,\eta)u^{2}-\left\\{\gamma(x,\eta)\mathcal{B}r(x,\eta)+F(x,\eta)\right\\}u\right)\frac{\partial_{u}g(x,u)}{g(x,u)}.\end{split}$
(5.11)

We say that a ${\mathcal{D}}_{[0,\infty)}(\mathcal{C})$-valued process
$(\xi_{t})_{t\geq 0}$ is a solution to the $(A,\xi_{0})$ martingale problem if
it has initial distribution $\xi_{0}$ and
$f(\xi_{t})-f(\xi_{0})-\int_{0}^{t}Af(\xi_{s})ds$ is a martingale (with
respect to the natural filtration) for all test functions $f$ as defined
above.

The lookdown processes have been carefully constructed so that observations
about the past spatial positions of individuals in the population do not give
us any information about the assignment of individuals to levels. In other
words, the dynamics of the lookdown process preserve the conditionally uniform
(or in the limit, conditionally Poisson) structure – if started with uniform
levels, levels are uniform at all future times. Moreover, if we average over
levels in the expression for the generator (equation (5.10) or (5.11)) we
recover the generator for the population process. Once this is verified (along
with some boundedness conditions) the Markov Mapping Theorem (Theorem A.1;
also see Etheridge and Kurtz (2019)) tells us that by “removing labels” from
the lookdown process $\xi$ we recover the population process $\eta$.

To make this precise, define the spatial projection maps
$\kappa^{N}:\mathcal{M}(\mathbb{R}^{d}\times[0,N])\to\mathcal{M}(\mathbb{R}^{d})$
by $\kappa^{N}(\xi^{N})(\cdot)=\xi^{N}(\cdot\times[0,N])/N$, and
$\kappa:\mathcal{M}(\mathbb{R}^{d}\times[0,\infty))\to\mathcal{M}(\mathbb{R}^{d})$
by $\kappa(\xi)(\cdot)=\lim_{u_{0}\to\infty}\xi(\cdot\times[0,u_{0}])/u_{0}$.
We will also need an inverse notion: for a measure $\xi^{N}$ on
$\mathbb{R}^{d}\times[0,N]$ and a $\sigma$-field $\mathcal{F}$, we say that
_$\xi^{N}$ is conditionally uniform given $\mathcal{F}$_ if $\kappa^{N}(\xi)$
is $\mathcal{F}$-measurable and for all compactly supported $f$,

$\displaystyle\mathbb{E}[e^{-\langle f,\xi\rangle}\;|\;\mathcal{F}]$
$\displaystyle=e^{-\langle H_{f}^{N},\kappa^{N}(\xi)\rangle},$ (5.12)

where

$H_{f}^{N}(x)=-N\log\frac{1}{N}\int_{0}^{N}e^{-f(x,u)}du.$

In other words, the $[0,N]$ components of $\xi$ are independent, uniformly
distributed on $[0,N]$, and independent of $\kappa^{N}(\xi)$. Similarly, for a
measure $\xi$ on $\mathbb{R}^{d}\times[0,\infty)$ we say that _$\xi$ is a
conditionally Poisson random measure given $\mathcal{F}$_ if $\kappa(\xi)$ is
$\mathcal{F}$-measurable and for all compactly supported $f$,

$\displaystyle\mathbb{E}[e^{-\langle f,\xi\rangle}\;|\;\mathcal{F}]$
$\displaystyle=e^{-\langle\int_{0}^{\infty}(1-e^{-f(x,u)})du,\kappa(\xi)(dx)\rangle}.$
(5.13)

In other words, $\xi$ is conditionally Poisson with Cox measure
$\kappa(\xi)\times\lambda$, where $\lambda$ is Lesbegue measure.

######  Proposition 5.3.

If $\widetilde{\eta}^{N}$ is a solution of the martingale problem of
Definition 2.4 with initial distribution $\eta^{N}_{0}$ then there exists a
solution $\xi^{N}$ of the $(A^{N},\xi^{N}_{0})$-martingale problem of
Definition 5.1 such that $\eta^{N}=\kappa^{N}\circ\xi^{N}$ has the same
distribution on $D_{\mathcal{M}_{F}(\mathbb{R}^{d})}[0,\infty)$ as
$\widetilde{\eta}^{N}$. Furthermore, for each $t$, $\xi^{N}_{t}$ is
conditionally uniform given $\mathcal{F}^{\eta^{N}}_{t}$ in the sense of
(5.12).

Similarly, if $\widetilde{\eta}$ is a solution of the limiting martingale
problem of Theorem 2.10 with initial distribution $\eta_{0}$ then there exists
a solution $\xi$ of the martingale problem of of Definition 5.2 such that
$\eta=\kappa\circ\xi$ has the same distribution on
$D_{\mathcal{M}_{F}(\mathbb{R}^{d})}[0,\infty)$ as $\widetilde{\eta}$.
Furthermore, $\xi_{t}$ is conditionally Poisson given $\mathcal{F}^{\eta}_{t}$
in the sense of (5.13).

Now we can present the main convergence theorem that is analogous to Theorem
2.10 for the population process.

######  Theorem 5.4.

Let $(\xi_{t}^{N})$ satisfy Definition 5.1 and assume that as $N\to\infty$,
$\theta\to\infty$ in such a way that $\theta/N\to\alpha$. Let
$\eta_{0}^{N}=\kappa(\xi_{0}^{N})$ and suppose also that
$\eta^{N}_{0}\to\eta_{0}$ in $\mathcal{M}_{F}(\mathbb{R}^{d})$, and that for
each $N$, $\xi^{N}_{0}$ is conditionally uniform given $\eta_{0}^{N}$ in the
sense of (5.12). Then, $(\xi_{t}^{N})_{t\geq 0}$ has a subsequence which
converges in distribution as $N\to\infty$ to a measure-valued process
$(\xi_{t})_{t\geq 0}$ with $\xi_{t}$ conditionally Poisson given
$\eta_{t}=\kappa(\xi_{t})$ for each $t$ in the sense of (5.13), that is a
solution to the martingale problem of Definition 5.2.

Both results are proved in Section 8.

###  5.2 Explicit construction of lines of descent

The main interest in using a lookdown construction for our population
processes is that it allows us to retain information about the relatedness of
individuals as we pass to the infinite population limit. In order to exploit
this, in this section we write down stochastic equations for the locations and
levels of individuals in the prelimiting lookdown model. We will then be able
to pass to the scaling limit. This provides an explicit description of the
solution to the limiting martingale problem of Definition 5.2 which will
enable us to identify all individuals in the current population that are
descendants of a given ancestor at time zero. In theory at least, this allows
us to recover all the information about genealogies relating individuals
sampled from the present day population. This idea draws on the notion of
“tracers”, popular in statistical physics and used in population genetics by a
number of authors including Hallatschek and Nelson (2008), Durrett and Fan
(2016), and Biswas et al. (2021).

We will construct the process using a Ulam-Harris indexing scheme. First, we
assign each individual alive at time 0 a unique label from $\mathbb{N}$.
Suppose an individual with label $a$ and level $u$ reproduces, and as a result
there are two individuals, one with level $u$ and one with a new level
$u_{1}>u$. The parent individual, previously labeled $a$, might be assigned
either level. We will track chains of descendant individuals forwards through
time by following levels, rather than individuals, and will call this a _line
of descent_. So, after reproduction, we give a new label to _only_ the
individual that is given the new level $u_{1}$, retaining the label $a$ for
the individual with the old level $u$. In this way, at each birth event, a
unique label is assigned to the resulting individual with the higher level,
and the label of an individual may change throughout its lifetime.

Concretely, then: for each label $a$ in $\mathcal{I}=\bigcup_{k\geq
1}\mathbb{N}^{k}$, let $\Pi_{a}$ be an independent Poisson process on
$[0,\infty)^{2}\times\mathbb{R}^{d}\times\\{0,1\\}$. The mean measure of each
$\Pi_{a}$ is a product of Lebesgue measure on $[0,\infty)^{2}$, the density of
the standard Gaussian on $\mathbb{R}^{d}$, and $(\delta_{0}+\delta_{1})/2$ on
$\\{0,1\\}$. It will also be convenient to suppose that for each label $a$ we
have an enumeration of the points in $\Pi_{a}$, so we may refer to “the
$j^{\text{th}}$ point in $\Pi_{a}$”, although the precise order of this
enumeration is irrelevant. If $(\tau,v,z,\kappa)$ is the $j^{\text{th}}$ point
in $\Pi_{a}$, then $\tau$ will determine a possible birth time, $v$ will
determine the level of the offspring, $z$ will determine the spatial
displacement of the offspring relative to the parent, $\kappa$ will be used to
determine whether parent or offspring is assigned the new level, and the new
label produced will be $a\oplus j$, i.e., the label $a$ with $j$ appended (so,
if $a=(a_{1},\ldots,a_{k})$ then $a\oplus j=(a_{1},\ldots,a_{k},j)$). Each
label $a$ has a birth time $\tau_{a}$, when it is first assigned, and a
(possibly infinite) death time $\sigma_{a}$, when its level first hits $N$.
For any $\tau_{a}\leq t\leq\sigma_{a}$ we denote by $X_{a}(t)$ and $U_{a}(t)$
the spatial location and level of the individual carrying label $a$ at time
$t$, respectively. Furthermore, define

$\displaystyle\eta^{N}_{t}=\frac{1}{N}\sum_{a:\tau_{a}\leq
t<\sigma_{a}}\delta_{X_{a}(t)}\qquad\text{and}\qquad\xi^{N}_{t}=\sum_{a:\tau_{a}\leq
t<\sigma_{a}}\delta_{(X_{a}(t),U_{a}(t))}.$

Now, since we have defined labels so that the level does not jump, $U_{a}$
satisfies (5.5) for $\tau_{a}\leq t\leq\sigma_{a}$, i.e.,

$\begin{split}&U_{a}(t)=U_{a}(\tau_{a})\\\
&\qquad{}+\int_{\tau_{a}}^{t}\left(c_{\theta}(X_{a}(s),\eta_{s})U_{a}(s)^{2}-b_{\theta}(X_{a}(s),\eta_{s})U_{a}(s)\right)ds,\end{split}$
(5.14)

and, of course, $\sigma_{a}=\inf\\{t\geq\tau_{a}:U_{a}(t)>N\\}$.

Potential reproduction events occur at times $\tau$ for each point
$(\tau,v,z,\kappa)\in\Pi_{a}$ with $\tau_{a}\leq\tau<\sigma_{a}$. (We say
“potential” since if the level of the resulting offspring is greater than $N$,
the event does not happen.) If this is the $j^{\text{th}}$ point in $\Pi_{a}$,
the potential new label is $a\oplus j$, the birth time is $\tau_{a\oplus
j}=\tau$, and the spatial displacement of the potential offspring is
$y(X(\tau-),z)$, where

$y(x,z):=\frac{1}{\theta}\vec{b}(x)+\frac{1}{\sqrt{\theta}}K(x)z,$

and $K(x)K^{T}(x)=\mathbf{C}(x)$.

Next we must choose the new level created at the birth event. We would like an
individual with level $u$ and at spatial position $x$ to produce offspring at
$y$ at instantaneous rate

$2\Big{(}1-\frac{u}{N}\Big{)}\theta\gamma(x,\eta)r(x+y,\eta).$ (5.15)

To do this we will associate the point $(\tau,v,z,\kappa)\in\Pi_{a}$ with
level $u+v\ell$, where $\ell$ is chosen so that the rate of appearance of
points in $\Pi_{a}$ with level below $N$, that is points with $v\ell<N-u$, is
given by (5.15). Since the mean measure of $\Pi_{a}$ is Lebesgue measure in
the $t$ and $v$ directions, we must take

$\ell(x,y,\eta)=\frac{N-u}{2(1-u/N)\theta\gamma(x,\eta)r(x+y,\eta)}=\frac{1}{2N^{-1}\theta\gamma(x,\eta)r(x+y,\eta)},$
(5.16)

and, using this, the (potential) new level is

$U_{a\oplus
j}(\tau)=U_{a}(\tau)+v\ell\big{(}X_{a}(\tau-),y(X_{a}(\tau-),z),\eta_{\tau-}\big{)}.$

If $U_{a\oplus j}(\tau)<N$, the new individual labeled $a\oplus j$ is
produced, and $\kappa$ determines which label, $a$ or $a\oplus j$, is
associated with the new location, so

$\displaystyle X_{a\oplus j}(\tau)$
$\displaystyle=X_{a}(\tau-)+(1-\kappa)y\big{(}X_{a}(\tau-),z\big{)}.$

On the other hand if $U_{a\oplus j}(\tau)\geq N$, then $X_{a}$ is unchanged
and $X_{a\oplus j}$ is undefined, so

$\displaystyle X_{a}(\tau)$ $\displaystyle=X_{a}(\tau-)+\kappa
y(X_{a}(\tau-),z)\mathbf{1}_{U_{a\oplus j}(\tau)<N}.$ (5.17)

Recall that the parental _individual_ always retains their spatial location,
so that $\kappa=0$ corresponds to the parent being assigned a new level, and
our line of descent switching to the offspring. Combining these observations,
$X_{a}$, for $\tau_{a}\leq t<\sigma_{a}$, solves the equation

$\displaystyle X_{a}(t)$
$\displaystyle=X_{a}(\tau_{a})+\int_{[\tau_{a},t)\times[0,\infty)\times\mathbb{R}\times[0,1]}y(X_{a}(\tau-),z)\kappa\mathbf{1}_{U_{a}(\tau)+v\ell(X_{a}(\tau-),y(X_{a}(\tau-),z),\eta_{\tau-})<N}d\Pi_{a}(\tau,v,z,\kappa).$

Although we have described the evolution of a line of descent only for a given
label (i.e., for $\tau_{a}\leq t<\sigma_{a}$), we can extend the definition to
times $0\leq t<\sigma_{a}$ by setting $X_{a}(t)$ equal to $X_{[a]_{t}}(t)$,
where $[a]_{t}$ is the label of the ancestor of label $a$ alive at time $t$,
and similarly for $U_{a}(t)$. It is then straightforward, albeit tedious, to
write down the time evolution of $(X_{a}(t),U_{a}(t))$ for all time back to
$t=0$ in terms of the driving Poisson processes.

######  Remark 5.5.

Although we have a single construction that couples the processes across all
$N$, unlike in Kurtz and Rodrigues (2011) the actual trajectories,
$X_{a}(\cdot)$, do not necessarily coincide for different values of $N$, since
they are affected by the whole population process. However, this does suggest
approximating the genealogies in the infinite density limit by simulating up
until a sufficiently high level that we have a good approximation to the
population process.

###  5.3 Limiting processes for lines of descent

The previous section constructed the lookdown process using the same
underlying Poisson processes $\\{\Pi_{a}\\}_{a\in\mathcal{I}}$ for different
values of $N$. As a result, if the spatial projections $\eta$ converge, then
individual lines of descent converge pointwise as $N\to\infty$. To see this,
first note that if the Poisson processes are fixed then the set of events with
which a given label $a\in\mathcal{I}$ is associated is also fixed – this is
the sequence $(\tau_{k},v_{k},z_{k},\kappa_{k})$ associated with the label
$a$. To conclude that the lines of descent converge, first, we clearly need
that the spatial projections $\eta$ converge. Supposing that they do, consider
how a line of descent $(X_{a}(t),U_{a}(t))$ evolves. It throws off a new line
of descent at a higher level when there is a point $(\tau,v,z,\kappa)$ in
$\Pi_{a}$ with $\tau>\tau_{a}$ and

$v<2\frac{\big{(}N-U_{a}(\tau)\big{)}}{N}\theta\gamma(X_{a}(\tau-),\eta_{\tau-})r\Big{(}X_{a}(\tau-)+y\big{(}X_{a}(\tau-),z\big{)},\eta_{\tau-}\Big{)}.$
(5.18)

Since the mean measure of the $v$ coordinate is Lebesgue measure,
$\theta/N\to\alpha$, and $q_{\theta}(x,dy)\to\delta_{x}(dy)$, this corresponds
in the limit to new lines of descent being thrown off according to a Poisson
process with intensity

$2\alpha\gamma(X_{a}(t),\eta_{t})r(X_{a}(t),\eta_{t})dt\times du.$

Now consider the location of the line of descent: at each birth event, with
probability one half the line of descent jumps to $X_{a}(t)+y$. Taking $g$ to
be a suitable test function on $\mathbb{R}^{d}$, and rewriting (5.18), when
the level is $u$ and the state of the population is $\eta$, the generator of
the spatial motion of the line of descent applied to $g(x)$ is

$\displaystyle\left(1-\frac{u}{N}\right)\gamma(x,\eta)\theta\int_{\mathbb{R}^{d}}r(x+y,\eta)(g(x+y)-g(x))q_{\theta}(x,dy)$
$\displaystyle\qquad{}=\left(1-\frac{u}{N}\right)\gamma(x,\eta)\bigg{\\{}\theta\int_{\mathbb{R}^{d}}(r(x+y,\eta)g(x+y)-r(x,\eta)g(x))q_{\theta}(x,dy)$
$\displaystyle\qquad\qquad{}-\theta\int_{\mathbb{R}^{d}}(r(x+y,\eta)-r(x,\eta))g(x)q_{\theta}(x,dy)\bigg{\\}}$
$\displaystyle\qquad{}\to\gamma(x,\eta)\left(\mathcal{B}(rg)(x)-g(x)\mathcal{B}(r)(x)\right),\qquad\text{as
}N,\theta\to\infty.$

Notice that the factors of 2 have cancelled, and that the result is
independent of $u$. Also recall that $r(x,\eta)$ depends on $\eta$ only
through $\rho_{r}\\!*\\!\eta(x)$, which is guaranteed to be smooth, so that
$\mathcal{B}(r)$ and $\mathcal{B}(gr)$ are well-defined.

We write out the differential operator above in more detail. Recall that
$\mathcal{B}g(x)=\sum_{i}\vec{b}_{i}\partial_{i}g(x)+\sum_{ij}\mathbf{C}_{ij}\partial_{ij}g(x)$,
and for the moment write $r(x)$ for $r(x,\eta)$, $\vec{b}(x)=\vec{b}$, and
$\mathbf{C}(x)=\mathbf{C}$ so that

$\displaystyle\mathcal{B}(rg)(x)-g(x)\mathcal{B}(r)(x)$
$\displaystyle=r(x)\sum_{i}\vec{b}_{i}\partial_{i}g(x)+2\sum_{ij}\partial_{i}r(x)\mathbf{C}_{ij}\partial_{j}g(x)+r(x)\sum_{ij}\mathbf{C}_{ij}\partial_{ij}g(x)$
$\displaystyle=r(x)\left\\{\left(\vec{b}+2\mathbf{C}\nabla\log
r(x)\right)\cdot\nabla
g(x)+\sum_{ij}\mathbf{C}_{ij}\partial_{ij}g(x)\right\\}.$ (5.19)

The only thing that remains is to describe how the levels change, but this is
immediate from applying limit (5.4) to equation (5.5).

We summarize the results in a proposition.

######  Proposition 5.6 (Line of descent construction).

Define $J(x,\eta)$ and $\beta(x,\eta)$ by

$\displaystyle r(x,\eta)\gamma(x,\eta)\mathbf{C}(x)=J(x,\eta)J(x,\eta)^{T}$
$\displaystyle\beta(x,\eta)=r(x,\eta)\gamma(x,\eta)\big{(}\vec{b}(x)+2\mathbf{C}(x)\nabla\log
r(x,\eta)\big{)}.$

Associate with each label $a\in\mathcal{I}=\cup_{k\geq 1}\mathbb{N}^{k}$ an
independent $d$-dimensional Brownian motion $W_{a}$ and an independent Poisson
process $R_{a}$ on $[0,\infty)^{2}$ with Lebesgue mean measure, and with
points ordered in some way. Given
$\eta_{0}\in\mathcal{M}_{F}(\mathbb{R}^{d})$, let $(x_{i},u_{i})$ be the
points of a Poisson process on $\mathbb{R}^{d}\times[0,\infty)$ with mean
measure $\eta_{0}\times\lambda$ (the product of $\eta_{0}$ and Lebesgue
measure). For each $i$, begin a line of descent with label $i$, location
$X_{i}(0)=x_{i}$, level $U_{i}(0)=u_{i}$, and birth time $\tau_{i}=0$.

Write $\tau_{a}$ for the birth time of the label $a$ and
$\sigma_{a}=\lim_{u_{0}\to\infty}\inf\\{t\geq 0:U_{a}(t)>u_{0}\\}$ the time
the level hits $\infty$. Suppose that the spatial locations and level of each
line of descent $a$ solve, for $\tau_{a}\leq t<\sigma_{a}$,

$\displaystyle\begin{split}X_{a}(t)&=X_{a}(\tau_{a})+\int_{\tau_{a}}^{t}\beta(X_{a}(s),\eta_{s})ds+\int_{\tau_{a}}^{t}J(X_{a}(s),\eta_{s})dW_{a}(s)\\\
U_{a}(t)&=U_{a}(\tau_{a})+\int_{\tau_{a}}^{t}\bigg{(}\alpha\gamma(X_{a}(s),\eta_{s})r(X_{a}(s),\eta_{s})U_{a}(s)^{2}\\\
&\qquad\qquad\qquad{}-\big{\\{}\gamma(X_{a}(s),\eta_{s})\mathcal{B}r(X_{a}(s),\eta_{s})+F(X_{a}(s),\eta_{s})\big{\\}}U_{a}(s)\bigg{)}ds.\end{split}$
(5.20)

Each point in each $R_{a}$ denotes a potential birth time for $a$: if the
$j^{\text{th}}$ point in $R_{a}$ is $(\tau,v)$, with
$\tau_{a}\leq\tau<\sigma_{a}$, then a new line of descent with label $a\oplus
j$ is produced, with birth time $\tau_{a\oplus j}=\tau$, location $X_{a\oplus
j}(\tau)=X_{a}(\tau)$, and level

$U_{a\oplus
j}(\tau)=U_{a}(\tau)+\frac{v}{2\alpha\gamma(X_{a}(\tau),\eta_{\tau})r(X_{a}(\tau),\eta_{\tau})},$

if this is finite. For any solution to the equations above, the processes
defined by

$\displaystyle\eta_{t}=\lim_{u_{0}\to\infty}\frac{1}{u_{0}}\sum_{a:\tau_{a}\leq
t<\sigma_{a}\;;\;U_{a}(t)<u_{0}}\delta_{X_{a}(t)}\qquad\text{and}\qquad\xi_{t}=\sum_{a:\tau_{a}\leq
t<\sigma_{a}}\delta_{(X_{a}(t),U_{a}(t))}$

are solutions to the martingale problems of Theorems 2.10 and 5.4,
respectively.

In particular, note that if $\alpha=0$, no new lines of descent are produced.
More precisely, comparing with (5.16), they are produced, but “at infinity”,
and their trace is seen in the spatial motion of the line of descent which
results from the production of these lineages.

###### of Proposition 5.6.

The fact that a solution to the system of equations (5.20) is a solution to
the martingale problem of Theorem 5.4 is an application of Itô’s theorem.
Furthermore, in Proposition 5.3 we showed that the conditional Poisson
property of $\xi_{0}$ is preserved (i.e., holds for $\xi_{t}$ for all $t$),
and so $(\eta_{t})_{t\geq 0}$ is well-defined, and furthermore that $\eta_{t}$
is a solution to the martingale problem of Theorem 2.10.

Proofs of the remainder of these results are in Section 8.

######  Remark 5.7.

The process we consider is similar to the state-dependent branching processes
of Kurtz and Rodrigues (2011), so one might expect that the proofs there would
carry over with little change. However, there is an important difference:
Recall that the level $U_{a}(t)$ of a line of descent evolves as

$\dot{u}=c_{\theta}(x,\eta)u^{2}-b_{\theta}(x,\eta)u,$ (5.21)

where $b_{\theta}(x,\eta)$ and $c_{\theta}(x,\eta)$ are defined in (5.3) and
(5.2) respectively. Note that $c_{\theta}(x,\eta)\geq 0$, while
$b_{\theta}(x,\eta)$ may take either sign. Assumptions 2.8 imply that
$c_{\theta}(x,\eta)$ is bounded, while $b_{\theta}(x,\eta)$, because of
$F(x,\eta)$, is bounded above but not necessarily below. In Kurtz and
Rodrigues (2011), $b_{\theta}$ was bounded above and $c_{\theta}$ was bounded
away from zero, so they noted that if $U_{a}(t)\geq b_{\theta}/c_{\theta}$ for
some label $a$, that line of descent would only move upwards from that time
onwards. Furthermore, coefficients did not depend on the state of the process
(i.e., on $\eta$), thus allowing the processes to be jointly and
simultaneously constructed for all values of $N$, with a pointwise embedding
of $(\xi^{N}_{t})_{t\geq 0}$ within $(\xi^{M})_{t\geq 0}$ for
$b_{\theta}/c_{\theta}<N<M$. In other words, individuals with levels above
$N>b_{\theta}/c_{\theta}$ at time $t_{0}$ do not affect $(\xi^{N}_{t})_{t\geq
t_{0}}$, thus allowing a comparison of the number of lines of descent below
level $u_{0}$ to a branching process. Although we have provided a joint
construction of $\xi^{N}$ for all $N$ in Section 5.2, it does not have this
monotonicity: for one thing, $b_{\theta}$ and $c_{\theta}$ depend on the
population process $\eta$ and so all individuals can affect all other ones
(even those with lower levels). Furthermore, in the deterministic case
$\theta/N$, and hence $c$, converges to zero, and so lines of descent with
arbitrarily high level may drift back downwards. Indeed, this must be the case
if the population persists, since in the deterministic case there is no
branching.

##  6 Proofs of convergence for nonlocal models

In this section we present formal proofs of the first two of our three scaling
limits. In Subsection 6.2 we prove Theorem 2.10, to obtain (both stochastic
and deterministic) limits in which interactions between individuals in the
population are nonlocal. In Subsection 6.3 we show how, in two important
examples in which the nonlocal limit is respectively a deterministic solution
to a non-local equation of reaction-diffusion type and a deterministic
solution to a nonlocal porous medium equation with an additional logistic
growth term, one can pass to a further limit to obtain a classical PDE.

###  6.1 Preliminaries

Below we will have frequent use for the quantity

$B^{\theta}_{f}(x,\eta)=\theta\int_{\mathbb{R}^{d}}(f(y)r(y,\eta)-f(x)r(x,\eta))q_{\theta}(x,dy).$
(6.1)

First, we prove Lemma 2.9.

###### of Lemma 2.9.

Here, we need to prove that $|\gamma(x,\eta)B^{\theta}_{f}(x,\eta)|$ is
bounded, uniformly over $x$ and $\eta$. First suppose that Condition 1 of
Lemma 2.9 is satisfied. We write

$\displaystyle r(y,\eta)f(y)-r(x,\eta)f(x)$
$\displaystyle=r(y,\eta)(f(y)-f(x))+(r(y,\eta)-r(x,\eta))f(x)$
$\displaystyle=r(y,\eta)\left(\sum_{i}(y-x)_{i}\partial_{x_{i}}f(x)+\sum_{ij}(y-x)_{i}(y-x)_{j}\partial_{x_{i}x_{j}}f(z_{1})\right)$
$\displaystyle\qquad{}+f(x)\left(\sum_{i}(y-x)_{i}\partial_{x_{i}}r(x,\eta)+\sum_{ij}(y-x)_{i}(y-x)_{j}\partial_{x_{i}x_{j}}r(z_{2},\eta)\right)$
$\displaystyle=\left(r(x,\eta)+\sum_{j}(y-x)_{j}\partial_{x_{j}}r(z_{3},\eta)\right)\left(\sum_{i}(y-x)_{i}\partial_{x_{i}}f(x)\right)$
$\displaystyle\qquad{}+r(y,\eta)\sum_{ij}(y-x)_{i}(y-x)_{j}\partial_{x_{i}x_{j}}f(z_{1})$
$\displaystyle\qquad{}+f(x)\left(\sum_{i}(y-x)_{i}\partial_{x_{i}}r(x,\eta)+\sum_{ij}(y-x)_{i}(y-x)_{j}\partial_{x_{i}x_{j}}r(z_{2},\eta)\right),$

for some $z_{i}=\kappa_{i}x+(1-\kappa_{i})y$. Integrating this against
$q(x,dy)$, we get that

$\displaystyle\bigg{|}\theta\int\left(r(y,\eta)f(y)-r(x,\eta)f(x)\right)q_{\theta}(x,dy)\bigg{|}$
$\displaystyle\qquad\leq\bigg{|}\sum_{i}\left(r(x,\eta)\partial_{x_{i}}f(x)+f(x)\partial_{x_{i}}r(x,\eta)\right)\theta\int(y-x)_{i}q_{\theta}(x,dy)\bigg{|}$
$\displaystyle\qquad\qquad{}+\bigg{|}f(x)\theta\int\sum_{ij}\partial_{x_{i}x_{j}}r(z_{2},\eta)(y-x)_{i}(y-x)_{j}q_{\theta}(x,dy)\bigg{|}$