# Weak lensing in the blue:
a counter-intuitive strategy for stratospheric observations

Mohamed M. Shaaban Dunlap Institute for Astronomy and Astrophysics, University
of Toronto, 50 St. George Street, Toronto, ON, Canada M5S 3H4  Department of
Physics, University of Toronto, 60 St. George Street, Toronto, ON, Canada M5R
2M8  Ajay S. Gill Dunlap Institute for Astronomy and Astrophysics, University
of Toronto, 50 St. George Street, Toronto, ON, Canada M5S 3H4  David A. Dunlap
Dept. of Astronomy and Astrophysics, University of Toronto, 50 St. George
Street, Toronto, ON, Canada M5S 3H4  Jacqueline McCleary Department of
Physics, Northeastern University, 360 Huntington Ave, Boston, MA  Richard J.
Massey Institute for Computational Cosmology, Durham University, South Road,
Durham DH1 3LE, UK  Centre for Extragalactic Astronomy, Department of Physics,
Durham University, Durham DH1 3LE, UK  Centre for Advanced Instrumentation
(CfAI), Durham University, South Road, Durham DH1 3LE, UK  Steven J. Benton
Department of Physics, Princeton University, Jadwin Hall, Princeton, NJ, USA
08544  Anthony M. Brown  Centre for Advanced Instrumentation (CfAI), Durham
University, South Road, Durham DH1 3LE, UK  Centre for Extragalactic
Astronomy, Department of Physics, Durham University, Durham DH1 3LE, UK
Christopher J. Damaren  University of Toronto Institute for Aerospace Studies
(UTIAS), 4925 Dufferin Street, Toronto, ON, Canada M3H 5T6  Tim Eifler
Department of Astronomy/Steward Observatory, University of Arizona, 933 North
Cherry Avenue, Tucson, AZ, USA 85721  Aurelien A. Fraisse  Department of
Physics, Princeton University, Jadwin Hall, Princeton, NJ, USA 08544  Spencer
Everett  Jet Propulsion Laboratory (JPL), California Institute of Technology,
4800 Oak Grove Drive, Pasadena, CA, USA 91109  Mathew N. Galloway  Institute
of Theoretical Astrophysics, University of Oslo, Blindern, Oslo 0315, Norway
Michael Henderson  University of Toronto Institute for Aerospace Studies
(UTIAS), 4925 Dufferin Street, Toronto, ON, Canada M3H 5T6  Bradley Holder
University of Toronto Institute for Aerospace Studies (UTIAS), 4925 Dufferin
Street, Toronto, ON, Canada M3H 5T6  Dunlap Institute for Astronomy and
Astrophysics, University of Toronto, 50 St. George Street, Toronto, ON, Canada
M5S 3H4  Eric M. Huff  Jet Propulsion Laboratory (JPL), California Institute
of Technology, 4800 Oak Grove Drive, Pasadena, CA, USA 91109  Mathilde Jauzac
Centre for Extragalactic Astronomy, Department of Physics, Durham University,
Durham DH1 3LE, UK  Institute for Computational Cosmology, Durham University,
South Road, Durham DH1 3LE, UK  William C. Jones  Department of Physics,
Princeton University, Jadwin Hall, Princeton, NJ, USA 08544  David Lagattuta
Centre for Extragalactic Astronomy, Department of Physics, Durham University,
Durham DH1 3LE, UK  Jason S.-Y. Leung  David A. Dunlap Dept. of Astronomy and
Astrophysics, University of Toronto, 50 St. George Street, Toronto, ON, Canada
M5S 3H4  Dunlap Institute for Astronomy and Astrophysics, University of
Toronto, 50 St. George Street, Toronto, ON, Canada M5S 3H4  Lun Li  Department
of Physics, Princeton University, Jadwin Hall, Princeton, NJ, USA 08544  Thuy
Vy T. Luu  Department of Physics, Princeton University, Jadwin Hall,
Princeton, NJ, USA 08544  Johanna M. Nagy  Department of Physics, Washington
University in St. Louis, 1 Brookings Drive, St. Louis, MO, USA, 63130
McDonnell Center for the Space Sciences, Washington University in St. Louis, 1
Brookings Dr., St. Louis, MO USA 63130  C. Barth Netterfield  David A. Dunlap
Dept. of Astronomy and Astrophysics, University of Toronto, 50 St. George
Street, Toronto, ON, Canada M5S 3H4  Dunlap Institute for Astronomy and
Astrophysics, University of Toronto, 50 St. George Street, Toronto, ON, Canada
M5S 3H4  Department of Physics, University of Toronto, 60 St. George Street,
Toronto, ON, Canada M5R 2M8  Susan F. Redmond Department of Physics, Princeton
University, Jadwin Hall, Princeton, NJ, USA 08544  Jason D. Rhodes  Jet
Propulsion Laboratory (JPL), California Institute of Technology, 4800 Oak
Grove Drive, Pasadena, CA, USA 91109  Andrew Robertson  Jet Propulsion
Laboratory (JPL), California Institute of Technology, 4800 Oak Grove Drive,
Pasadena, CA, USA 91109  Jürgen Schmoll  Centre for Advanced Instrumentation
(CfAI), Durham University, South Road, Durham DH1 3LE, UK  Ellen Sirks
Institute for Computational Cosmology, Durham University, South Road, Durham
DH1 3LE, UK  Suresh Sivanandam  David A. Dunlap Dept. of Astronomy and
Astrophysics, University of Toronto, 50 St. George Street, Toronto, ON, Canada
M5S 3H4  Dunlap Institute for Astronomy and Astrophysics, University of
Toronto, 50 St. George Street, Toronto, ON, Canada M5S 3H4

###### Abstract

The statistical power of weak lensing measurements is principally driven by
the number of high-redshift galaxies whose shapes are resolved. Conventional
wisdom and physical intuition suggest this is optimised by deep imaging at
long (red or near-IR) wavelengths, to avoid losing redshifted Balmer-break and
Lyman-break galaxies. We use the synthetic Emission Line (‘EL’)-COSMOS
catalogue to simulate lensing observations using different filters, from
various altitudes. Here were predict the number of exposures to achieve a
target $z\gtrsim 0.3$ source density, using off-the-shelf and custom filters.
Ground-based observations are easily better at red wavelengths, as (more
narrowly) are space-based observations. However, we find that SuperBIT, a
diffraction-limited observatory operating in the stratosphere, should instead
perform its lensing-quality observations at blue wavelengths.

##  1 Introduction

The standard model of cosmology predicts that most of the energy density of
the universe is dominated by dark energy and dark matter, both of which
possess exotic physical properties that have only been weakly constrained. As
a result, constraining the nature of dark energy and dark matter are
considered two of the greatest challenges in contemporary science. Weak
gravitational lensing is a powerful tool for probing the physics of dark
matter and dark energy (Weinberg et al., 2013). Weak lensing occurs when light
radiated from distant galaxies changes path as it passes through a foreground
gravitational field. This change of path results in an apparent distortion of
the galaxy image. Measuring this change in shape allows us to reconstruct the
gravitational field and thus the distribution of (all) mass (Mandelbaum,
2018). Mapping the large-scale structure of dark matter can constrain
cosmological parameters (Haiman et al., 2001; Holder et al., 2001); mapping
its distribution around specific entities such as galaxy clusters can also
constrain cosmological parameters or test dark matter physics (Massey et al.,
2010; Harvey et al., 2015; Robertson et al., 2019). In practice, the original
shape of the distorted galaxy is unknown, so weak lensing signals can only be
inferred via statistical methods applied to ensembles of galaxies. With the
weak lensing signal of a galaxy cluster constrained inside a small patch of
sky, the signal-to-noise of cluster mass maps is driven by the spatial density
of galaxies whose shapes can be measured and averaged.

Conventional wisdom dictates that one can maximize the density of high-redshft
galaxies by observing at the red end of the spectrum. The physical intuition
behind this is the desire to avoid losing redshifted Balmer-break and Lyman-
break galaxies. However, the intuition can be violated by specific observing
conditions, environment, and instrumental configurations. One example of an
instrument for which this convention is violated is the Super-pressure
Balloon-borne Imaging Telescope (SuperBIT). SuperBIT is a balloon-borne,
diffraction-limited, wide-field, near-ultraviolet to near-infrared observatory
designed to exploit the stratosphere’s space-like conditions to perform deep
lensing observations of galaxy clusters (Gill et al., 2020; Romualdez et al.,
2019). The telescope has a $15\arcmin\times 23\arcmin$ field of view that is
designed to capture a 4 Mpc diameter galaxy cluster at $z\geq 0.2$ in a single
exposure. SuperBIT has 4 broadband filters and 2 ultra-broad band filters
covering a wavelength range from 300 nm to 900 nm. SuperBIT plans to use
either its broad blue band B or ultra-broad luminous band $L$ for its deep
lensing observations. The purpose of this paper is to provide justification
for this seemingly unconventional choice by simulating artificial SuperBIT
observations of the COSMOS survey (Scoville et al., 2007) and comparing the
relative performance of the bands.

##  2 Method

Any figure of merit to compare the performance of different filters is likely
to track some aspect of image depth as a function of exposure time. For a weak
lensing measurement, obtained statistically from ensembles of galaxies, a
natural choice is the density of high redshift galaxies that are both bright
enough to be detected and large enough to be resolved. We shall equivalently
adopt the observing time required to reach a target threshold in the density
of such galaxies.

To compute this figure of merit, we simulate observations in a range of
instrumental configurations and observing environments. The general procedure
is as follows. First we fix instrument-specific parameters for the telescope
optics and science detector (subsection 2.1). Then we specify the observing
environment and determine the flux from sky background that reaches the
detector in all of the instrument’s bands (subsection 2.2). Finally we model
the flux reaching the detector from distant galaxies (subsection 2.3), using
galaxies’ spectral energy distributions (SEDs) from (Saito et al., 2020)’s
synthetic Emission Line COSMOS (EL-COSMOS) catalogue and morphologies from
Mandelbaum et al. (2012)’s COSMOS-Galsim catalogue. We apply cuts similar to
those in a real cluster lensing analysis, requiring

  * •

galaxy redshift $z>0.3$

  * •

galaxy flux signal-to-noise ratio SNR$\geq 5$ inside a fixed $3\times 3$ pixel
aperture

  * •

galaxy half-light radius $r_{e}$ greater than the Full Width at Half Maximum
(FWHM) of the instrument’s Point Spread Function (PSF).

Our figure of merit is calculated from the density of galaxies remaining after
these cuts, as a function of exposure time.

###  2.1 Instrument Specifications

In this analysis we consider one set of system optics, with two different
configurations of detector throughput and optical transmission. Our optics
model matches the SuperBIT telescope (Romualdez et al., 2020): a diffraction-
limited 500 mm diameter telescope with 38% of the collecting area obscured by
a secondary mirror, and a focal plane Gaussian jitter with
$\sigma=0.05\arcmin\arcmin$. The first configuration of detector throughput
(SuperBIT configuration) comprises the nominal SuperBIT science detector: the
Sony IMX455 CMOS camera, with dark current of 0.002 electrons/s/pixel, read
noise of 1.5 electrons/pixel, and sensitivity in the 300-900 nm range measured
and presented in Gill et al. (2022). The second configuration (flat
configuration) is a fictional configuration with noise specifications
identical to the SuperBIT configuration, but with a total transmission
multiplied by quantum efficiency of 0.7 $e^{-}$/photon, irrespective of
wavelength. This is to enable us to quantify the extent to which our results
are dependent on the integrated QE curve shape. Both detector configurations
adopt a plate scale of 0.141$\arcmin\arcmin$ per pixel. For both
configurations we model the system’s Point Spread Function (PSF) as the
convolution of a diffraction limited airy disk and a Gaussian kernel
representing the system jitter. A summary of the relevant instrument
specifications is presented in Table 1.

The final instrument specification of relevance to this analysis is the choice
of photometric bandpasses. We perform this analysis in 6 bands. Four are
SuperBIT’s nominal photometric bands: U(330–430 nm), B(370–570 nm), G(520–700
nm), and R(640–800 nm). The remaining two are broad bands, one of which is a
near-red band that was the originally proposed deep lensing observation band:
S(530–830 nm) that has been scrapped as a result of this analysis. The other
is the near-blue band currently considered the deep lensing observation band
L(400–700 nm). This band was determined using a band optimization tool that
will be publicly released in an upcoming publication. The fully integrated
bandpass curves for the SuperBIT configuration are presented in Figure 1. Note
that at the time of this analysis, the SuperBIT science filters were still in
the manufacturing process so only the median in-band throughput was known. As
a result, the filters are taken to be boxcar with amplitude equal to the
median throughput as provided by the vendor ($\geq$0.9) and an out-of-band
transmission of 0.

Aperture | 500 mm ($f$/# = 11)
---|---
Obscuring Fraction | 0.38
Collecting Area ($A$) | 0.121 m2
Focal Plane Jitter (1$\sigma$) | 0.05$\arcmin\arcmin$
Bandwidth | 300-900 nm (See Figure 1)
Plate Scale | 0.141$\arcmin\arcmin$/pixel
Dark Current at -20 C∘ ($N_{d}$) | 0.002 $e^{-}/s$/pixel
Read Noise ($N_{r}$) | 1.5 $e^{-}$
Table 1: Summary of SuperBIT’s relevant specifications. Figure 1: The overall
SuperBIT configuration bandpass which is constructed by taking the product of
the throughput of all of SuperBIT’s optical elements, filter transmission, and
science camera quantum efficiency. The top panel shows the bandpass for the
photometric bands U(330-430 nm), B(370-570 nm), G(520-700 nm), and R(640-800
nm). The bottom panel shows the bandpass for the broad bands L(400-700 nm) and
S(530-830 nm). The flat configuration bandpass would be simple a series of
horizontal lines at 0.5.

###  2.2 Observing Environment

The next step in the analysis is to model the observing environment.
Specifically, we need to model the sky background level and atmospheric
transmission for both instrument configurations in three distinct
environments: the ground, the stratosphere, and low-earth orbit. While
SuperBIT will only collect data in the stratosphere, we include the other more
familiar environments to serve as a reference for the reader.

For the model of sky background in the stratosphere, we use the measured
average values from Gill et al. (2020). Specifically, we take the flux values
at the pivot wavelength $\lambda_{p}$ of each band in which the sky background
was measured, defined as

$\lambda_{p}^{2}=\frac{\int(e_{\lambda}\lambda)~{}\mathrm{d}\lambda}{\int(e_{\lambda}/\lambda)~{}\mathrm{d}\lambda}$
(1)

where $e_{\lambda}$ is the product of the total in-band transmission and
quantum efficiency as a function of wavelength, in units of $e^{-}/$photon. We
then linearly interpolate the flux as a function of wavelength between the
data points to generate an approximate model of the background SED in the
stratosphere. While this method has its limitations, it suffices for this
analysis as the bands of interest are similar in bandwidth and pivot
wavelength to the bands in which the measurements were performed. For the sky
background model in low earth orbit, we use the high sky background data
publicly available in the WFC3 manual for the Hubble Space Telescope (HST)
(Dressel, 2022). Finally, for the model of sky background from the ground, we
use the publicly available optical background data for the Gemini
Observatory111The Gemini Observatory background data used can be found at
http://www.gemini.edu/observing/telescopes-and-sites/sites.

The atmospheric transmission in the bands of interest is taken to be order
unity in the stratosphere and low-earth orbit. For the ground, we calculated
the atmospheric transmission using the MODTRAN webapp 222The MODTRAN webapp
can be found here http://modtran.spectral.com/modtran_home for an altitude of
5 km (Figure 2).

Figure 2: Atmospheric transmission at an altitude of 5 km above sea level in
desert conditions at mid-latitude winter, as calculated by the MODTRAN webapp.

All three background models are presented in Figure 3. Note that the ground
model is an empirical one and therefore the atmospheric transmission is
already folded in. We only need to use our modelled transmission curve for
simulating the observations of astronomical sources from the ground.

Figure 3: The background flux density in units of $\frac{\text{photons}}{s\
nm\ m^{2}\ \text{arcsec}^{2}}$ at the aperture for the three environments
considered.

Once the environments are modeled, we calculate the total number of electrons
per second, $N_{b}$, generated in each pixel of the detector by the sky
background, as

$N_{b}=AP^{2}\int f_{\lambda}e_{\lambda}~{}\mathrm{d}\lambda~{},$ (2)

where $A$ is the total collecting area, $P$ is the detector’s pixel scale in
arcseconds/pixel, $f_{\lambda}$ is the source flux density as a function of
wavelength, and $e_{\lambda}$ is the product of the total in band transmission
and quantum efficiency as a function of wavelength in units of $e^{-}/$photon.
The calculated background electron flux for all configurations is presented in
Figure 4. We note that the background flux per pixel for SuperBIT in the
stratosphere is significantly higher in $S$ than $L$ despite the comparable
bandwidth. This long-wavelength sky background plays a major role in
SuperBIT’s differential performance at red and blue wavelengths, as we will
show in section 3.

Figure 4: The modelled number of electrons per second per pixel, generated on
a detector by the sky background, for the SuperBIT (blue) and flat (orange)
instrument configurations, in all three observing environments.

###  2.3 Source Galaxy Population

####  2.3.1 Which galaxies are bright enough to be detected?

The publicly available EL-COSMOS catalogue (Saito et al., 2020) provides
spectral energy distribution (SED) models of more than 500,000 galaxies from
the COSMOS2015 catalogue (Laigle et al., 2016). For each galaxy in the central
region333The EL-COSMOS catalogue is divided into 9 area coverage zones for
ease of data management. We use the central zone (zone 5), which is already
$\sim$1.9 times larger than the field of view of SuperBIT. Parts of the
analysis were repeated on the other zones with no qualitative change in our
results, indicating no zone-selection bias. of EL-COSMOS, we convert the
provided SEDs from energy units to photon units, then calculate the total
number of electrons per second, $N_{s}$ that will be generated on the detector

$N_{s}=A\int f_{\lambda}T_{\lambda}e_{\lambda}~{}\mathrm{d}\lambda~{},$ (3)

where $A$ is the total collecting area, $f_{\lambda}$ is the source flux
density as a function of wavelength, $T_{\lambda}$ is the atmospheric
transmission as a function of wavelength and $e_{\lambda}$ is the product of
the total in-band transmission and quantum efficiency as a function of
wavelength, in units of $e^{-}/$photon. We shall distribute the photons
between pixels in section 2.3.2. As intuition suggested, for the flat
instrument configuration, the median source flux is higher in the redder S
band than the bluer L band for all three observing environments (see Figure
5). However, once the wavelength-dependence of SuperBIT’s quantum efficiency
is taken into account, the median on-detector flux is higher in L. As will be
shown in section 3, this will further separate the performance of SuperBIT
between red and blue wavelengths.

Figure 5: The median number of electrons per second from galaxies in the EL-
COSMOS catalogue for the SuperBIT (blue) and flat (orange) configurations for
all three observing environments. Note that due to identical atmospheric
transmission curves, the values in orbit and stratosphere are identical.

To predict the detection SNR of each galaxy in a $3\times 3$ pixel
aperture444To check that a 3$\times$3 pixel aperture is a reasonable choice,
note that galaxies’ median half-light radii are within 0.190$\arcmin\arcmin$
and 0.220$\arcmin\arcmin$ for all relevant redshifts, the SuperBIT PSF ranges
from 0.200$\arcmin\arcmin$-0.400$\arcmin\arcmin$, and that the SuperBIT
detector pixel scale is 0.141$\arcmin\arcmin$., we calculate

$\text{SNR}=\frac{nt\mathcal{F}_{p}N_{s}}{\sqrt{nt(\mathcal{F}_{p}N_{s}+pN_{b}+pN_{d})+nN_{r}^{2}}}~{},$
(4)

where $n$ is the total number of stacked exposures, $t$ is the duration of a
single exposure in seconds, $\mathcal{F}_{p}$ is the fraction of the source
intensity that lands inside $p=3\times 3=9$ pixels on the focal plane, $N_{s}$
is the total number of electrons per second from the observed source, $N_{b}$
is the number of electrons per second per pixel from the sky background,
$N_{d}$ is the dark current in electrons per pixel, and $N_{r}$ is the read
noise in electrons per pixel. We determined {$N_{s},N_{b},N_{d},N_{r}$} in
previous sections, and the duration of each exposure is dictated by SuperBIT’s
pointing and tracking specifications to be $t=300$s (see Romualdez et al.,
2020). The number of such exposures, $n$, will be our independent variable.
Our figure of merit will be the value of $n$ that guarantees a density of
galaxies with SNR$\geq$5 above a threshold (with low $n$ being better).

####  2.3.2 Assumptions about galaxy morphologies and sizes

To evaluate the figure of merit from equation 4, we still need to estimate the
parameter $\mathcal{F}_{p}$ for each galaxy. This require knowledge of galaxy
morphology that the EL-COSMOS catalog does not provide. We instead exploit
fits by Mandelbaum et al. (2012) of each COSMOS galaxy’s surface brightness
$I(r)$ to a Sérsic (1963) profile

$I(r)=I_{e}\exp\left\\{-b_{n}\left[\left({\frac{r}{r_{e}}}\right)^{1/n}-1\right]\right\\}~{},$
(5)

where $r_{e}$ is the radius inside which half of the total light of a galaxy
is emitted (also known as the half-light radius), $I_{e}$ is the source
intensity at the half-light radius, $n$ is the Sérsic index parameter, and
$b_{n}$ is a constant determined by the source index. Since our analysis will
use all galaxies behind a galaxy cluster at fiducial redshift $z=0.3$, we
determine the median values of $r_{e}$ and $n$ as a function of redshift (see
Figure 6).

Figure 6: The median Sérsic half light radius (top panel) and index (bottom
panel) of galaxies behind a fixed redshift, as a function of that redshift.

We calculate the median value of $\mathcal{F}_{p}$ as a function of cluster
redshift in all SuperBIT bands, assuming that galaxies have Sérsic profiles
with parameters drawn from Figure 6. We convolve these with a PSF, which is
modelled as a diffraction-limited Airy disk centered at the filter’s pivot
wavelength. In the stratospheric environment, the PSF is also convolved with a
Gaussian kernel of width $\sigma=0.05\arcsec$ to represent pointing jitter
(Romualdez et al., 2020). We draw both in an 11$\times$11 pixel window, with
each pixel over-sampled by a factor of 141$\times$141 (so each subsampled
pixel is conveniently 1 milliarcsecond on a side). We then convolve the galaxy
and PSF models and integrate the fraction of flux inside the central
3$\times$3 pixel aperture (see Figure 7). As one would expect for a fixed
aperture, the resulting value of $\mathcal{F}_{9}$ inversely scales with
wavelength (see Figure 8), because bluer bands will have a smaller diffraction
limited PSF. This effect also contributes to the red-blue performance
discrepancy. 555It is important to note that this step induces a non-trivial
selection bias in our simulations. In an actual observation the value of p
would be dynamically allocated on a source by source basis. However, this step
is unavoidable as there currently does not exist public catalogues with
sufficient depths to simulate SuperBIT observations that also provide
morphological fits. As a result, we make the most conservative choices to
ensure that the selection bias is in favour of larger PSFs (the redder bands)
and is therefore dampening the boast in performance gained by the blue bands.
The aforementioned choices are: (1) Taking the median source shape of all
galaxies behind the cluster rather than in a redshift bin thus putting more
weight on smaller sources. (2) Fixing $p=9$ rather than dynamically allocating
p on a source by source basis thus favouring PSFs with a size closer to 0.423
″.

Figure 7: Schematic illustration of the process used to calculate the fraction
$\mathcal{F}_{9}$ of each galaxy’s flux that falls within a $3\times 3$ pixel
aperture on the focal plane. Figure 8: The fraction $\mathcal{F}_{9}$ of a
galaxy’s flux that falls within a $3\times 3$ pixel aperture, as a function of
cluster redshift for median Sersíc morphology.

We can now use this information to calculate the number of background sources
detected by SuperBIT when observing a cluster at $z=0.3$. The result is shown
in Figure 9.

Figure 9: The number of sources with $z>0.3$ that SuperBIT can expect to
detect if it was to observe the EL-COSMOS zone5 as a function of exposure
time. The results are presented in all three relevant environments for each
one of the bands presented in Figure 1.

###  2.4 Which galaxies are large enough to be resolved?

Measurements of weak lensing are only possible for galaxies whose shapes are
well resolved. Typically, weak lensing analyses impose size cuts of the form
$({R_{\mathrm{gal}}}/{R_{\mathrm{PSF}}})>r$, where $R_{\mathrm{gal}}$ is some
measure of galaxy size, $R_{\mathrm{PSF}}$ is some measure of the PSF size,
and $r$ is a threshold ratio that depends upon mission systematic errors and
knowledge of the PSF (Leauthaud et al., 2007; Massey et al., 2013; Zuntz et
al., 2018). We perform a version of this size cut using the galaxy half-light
radii from Mandelbaum et al. (2012) as $R_{\mathrm{gal}}$, and the PSF FWHM as
$R_{\mathrm{PSF}}$. Since this analysis is focused on relative performance of
bands, the value of $r$ can be arbitrary, so we set it to be unity. This is
equivalent to labeling a source as resolved if its half light radius is larger
than the FWHM of the in-band PSF. Using these assumptions, we calculate the
ratio of resolved to detected sources and present the results in Figure 10.

There is an important consideration when evaluating the ratios of resolved to
detected galaxies in Figure 10. The Mandelbaum et al. (2012) catalog is
shallower than EL-COSMOS, and the relationship between galaxy size and
brightness is complex. Therefore, the intrinsic ratio of resolved to detected
galaxies would likely differ between the Mandelbaum et al. (2012) and EL-
COSMOS catalogs, had the latter provided galaxy size parameters. This means
that the ratio is not directly transferable between catalogues. However, we
can reasonably expect that the relative relationship of the ratio of resolved
to detected galaxies between bandpasses to change very slowly as a function of
catalogue depth (this assumption is partially supported by Figure 10).

Therefore, we can transfer the relative ratios to EL-COSMOS catalogue by
folding them into the curves in Figure 9. It is important to note that this
transfer renders the absolute scale of the y-axis meaningless. Instead of a
robust physical prediction of what SuperBIT could expect to observe in flight,
the $y$-axis is abstracted to a metric that only possesses comparative
information. However, we did calculate the galaxy densities predicted for HST
observations, and check that they roughly match the density (as a function of
magnitude) in the Leauthaud et al. (2017) catalogue. That suggests the
normalisation can be meaningful.

To avoid the misinterpretation of the results as a forecast we perform the
following normalization. Since SuperBIT expects to observe each cluster for
order 3 hours (36 exposures), our normalization will be chosen such that 3
hours of SuperBIT observations in the stratosphere in the L band will have a
metric value of 1.

Figure 10: The ratio of resolvable to detectable sources in the Mandelbaum et
al. (2012) catalogue as a function of redshift for each of SuperBIT’s bands.
Figure 11: The value of the performance metric for each SuperBIT band as a
function of time, in all three relevant observing environments.

##  3 Results

We now have everything we need to evaluate the relative performance of the
bands. To do so, we iterate over exposure time by varying the number $n$ of
300s exposures and calculating the SNR of each EL-COSMOS source at each
iteration. We simulate an observation of a cluster at $z=0.3$ by discarding
all sources with redshift $z\leq 0.3$, SNR$\leq$5, and
$R_{\mathrm{gal}}<R_{\mathrm{PSF}}$, and count the total number of remaining
sources in each band as a function of $n$. This simulates the number of
detected, resolved sources in each band as a function of exposure time (Figure
11).

We recover conventional wisdom that lensing measurements from the ground are
most efficient in red bands. For observations from the stratosphere however,
we find that bluer L (400–700 nm) and B (370–570 nm) bands are more efficient,
because of the diffraction-limited PSF size compared to the ground, and the
near-IR sky background compared to low-Earth orbit. As a result of this
analysis, the SuperBIT team has decided to avoid using the redder S band all
together and will primarily perform its deep cluster observations in the B
band, which has simultaneous utility for estimating photometric redshifts.

###  3.1 Systematic Errors

We recognise that the redshift uncertainty in the EL-COSMOS catalogue
increases for dimmer objects. In order to quantify the extent to which this
increased uncertainty affects our results, we repeat our analysis on the
catalogue with various HSC i-band magnitude cuts. We find that the performance
metric ratios between bands is mostly unaffected, as shown in Figure 12.

Figure 12: The ratio of the metric value for the $B$ and $S$ bands in the
stratospheric environment as a function number of 300 second exposures with
varying brightness cuts applied on the catalogue. This figure demonstrates
that the bias introduced into our result due to the increasing uncertainty in
source redshift is negligible.

Additionally, this analysis does not simulate two systematic errors commonly
associated with weak lensing observations of galaxy clusters: source blending
and background noise due to intracluster light (ICL). We believe that these
effects do not appreciably alter our conclusions for the following reasons:

  1. 1.

Blending effects become more important for observations with larger PSFs.
Inclusion of blending would only strengthen the case for observations in the
bluer $B$ band, as its smaller PSF enhances the ability to resolve galaxies.

  2. 2.

Weak lensing measurements are most sensitive to mass distributions projected
at 1 Mpc from the cluster center. With SuperBIT’s large field of view, the
majority of sources used in WL analysis will lie even farther than this from
the cluster’s centre. At this radius from the cluster center, the ICL is
neglible ($<$ 28 mag arcsec-2, see Chu et al., 2022), and should not
appreciably decrease the observation depth or increase photometric error.

##  4 Discussion

This analysis has involved simulating a variety of observing environments,
galaxy sizes, detector configurations, and optical bandpasses. These provide
sanity checks in familiar observing environments. They are also physically
informative, to test how much each effect contributes to our observed red-blue
performance discrepancy. To investigate these effects’ relative contributions,
we toggle the effects creating a variety of effect-free configurations. The
configurations considered are the following:

  * •

No size cut ($r=0$): In this configuration, we set the size cutoff ratio $r$
to 0 which is equivalent to treating all detected sources as resolved and
therefore bluer bands lose the performance advantage gained from having a
smaller PSF.

  * •

No QE shape (Flat): In this configuration, the integrated system’s quantum
efficiency curve is taken to be flat at 0.5 $e^{-}$/photon eliminating the
performance advantage blue bands gain due to SuperBIT s superior sensitivity
in the blue.

  * •

No background flux ($N_{b}=0$): In this configuration, the background flux
$N_{b}$ is set to 0 thus eliminating the performance advantage gained by blue
bands due to the bright red background in the stratosphere.

  * •

No source smearing ($\mathcal{F}_{9}=1$): In this configuration, we enforce
that the entire flux from the source lands within the aperture thus
eliminating the advantage gained by bluer bands due to lower smearing
resulting from smaller PSFs.

  * •

Raw flux density: In this configuration, all the effects mentioned above are
turned off (i.e. $r=0$ with a flat quantum efficiency, no background flux, and
no image smearing).

Each of the above effects is toggled on and off, and the metric value at three
hours (36 nominal exposures) is calculated in the stratospheric environment.
The results of this procedure are presented in Table 2 and Table 3.

| SuperBIT | r = 0 | Flat |  $N_{b}$ = 0 |  $\mathcal{F}_{9}$ = 1 | Raw
---|---|---|---|---|---|---
U | 0.4 | 0.56 | 0.61 | 0.69 | 0.83 | 1.73
B | 0.98 | 1.59 | 0.96 | 1.69 | 1.59 | 3.43
G | 0.67 | 1.41 | 0.63 | 1.31 | 1.16 | 3.58
R | 0.32 | 0.87 | 0.34 | 0.75 | 0.62 | 3.2
S | 0.6 | 1.4 | 0.57 | 1.34 | 1.03 | 3.8
L | 1 | 1.86 | 0.97 | 1.8 | 1.57 | 3.73
Table 2: The performance metric value after 36 exposures for each SuperBIT band with each effect independently toggled off in the stratosphere.  | r = 0 | Flat |  $N_{b}$ = 0 |  $\mathcal{F}_{9}$ = 1 | Raw
---|---|---|---|---|---
U | 1.40 | 1.53 | 1.73 | 2.08 | 4.33
B | 1.62 | 0.98 | 1.72 | 1.62 | 3.50
G | 2.10 | 0.94 | 1.96 | 1.73 | 5.34
R | 2.72 | 1.06 | 2.34 | 1.94 | 10.00
S | 2.33 | 0.95 | 2.23 | 1.72 | 6.33
L | 1.86 | 0.97 | 1.80 | 1.57 | 3.73
Table 3: This table is simply the result of dividing columns 2-6 in Table 2
with the the corresponding column 1 entry. This highlights how each each
effect contributes to the final result.

This result shows that while all of these effects favour bluer bands, the
dominant effects resulting in the superior blue performance are the
background, and the PSF size difference at the half meter diffraction limit.
Surprisingly,SuperBIT’s superior QE in the blue plays only a minor role.

The SuperBIT forecasting analysis presented in the upcoming publication
McCleary et al. 2022, strongly supports the central thesis of this paper. The
analysis found that when weak lensing pipeline was run on simulated B images
it was able to use 7447 background sources as opposed to the 2578 it was able
to use in S images. This result is summarised in Figure 13.

Figure 13: Median tangential shear profiles over 30 simulated observations of
a $4.1\times 10^{14}$ Msol galaxy cluster at $z=0.45$. Error bars represent
the standard deviation of the mean tangential shear in a radial bin. The red
line is the theoretical NFW shear input to all simulations. The variance is
more than three times larger in S than B or L. Reproduced from the upcoming
publication McCleary et al. 2022

##  5 Acknowledgements

Support for the development of SuperBIT is provided by NASA through APRA grant
NNX16AF65G. Launch and operational support for the sequence of test flights
from Palestine, Texas are provided by the Columbia Scientific Balloon Facility
(CSBF) under contract from NASA’s Balloon Program Office (BPO). Launch and
operational support for test flights from Timmins, Ontario are provided by the
Centre National d’Études Spatiales (CNES) and the Canadian Space Agency (CSA).

JR, EH, and SE are supported by JPL, which is run under a contract by Caltech
for NASA. Canadian coauthors acknowledge support from the Canadian Institute
for Advanced Research (CIFAR) as well as the Natural Science and Engineering
Research Council (NSERC). The Dunlap Institute is funded through an endowment
established by the David Dunlap family and the University of Toronto. UK
coauthors acknowledge funding from the Durham University Astronomy Projects
Award, STFC [grant ST/P000541/1], and the Royal Society [grants UF150687 and
RGF/EA/180026].

## References

  * Chu et al. (2022) Chu, A., Sarron, F., Durret, F., & Márquez, I. 2022, A&A, 666, A54, doi: 10.1051/0004-6361/202243504
  * Dressel (2022) Dressel, L. 2022, “Wide Field Camera 3 Instrument Handbook, Version 14.0” (Baltimore: STScI)
  * Gill et al. (2020) Gill, A., Benton, S. J., Brown, A. M., et al. 2020, The Astronomical Journal, 160, 266
  * Gill et al. (2022) Gill, A. S., Shaaban, M. M., Tohuvavohu, A., et al. 2022, arXiv preprint arXiv:2207.13052
  * Haiman et al. (2001) Haiman, Z., Mohr, J. J., & Holder, G. P. 2001, ApJ, 553, 545, doi: 10.1086/320939
  * Harvey et al. (2015) Harvey, D., Massey, R., Kitching, T., Taylor, A., & Tittley, E. 2015, Science, 347, 1462, doi: 10.1126/science.1261381
  * Holder et al. (2001) Holder, G., Haiman, Z., & Mohr, J. J. 2001, ApJ, 560, L111, doi: 10.1086/324309
  * Laigle et al. (2016) Laigle, C., McCracken, H. J., Ilbert, O., et al. 2016, ApJS, 224, 24, doi: 10.3847/0067-0049/224/2/24
  * Leauthaud et al. (2007) Leauthaud, A., Massey, R., Kneib, J.-P., et al. 2007, The Astrophysical Journal Supplement Series, 172, 219, doi: 10.1086/516598
  * Leauthaud et al. (2017) Leauthaud, A., Saito, S., Hilbert, S., et al. 2017, Monthly Notices of the Royal Astronomical Society, 467, 3024
  * Mandelbaum (2018) Mandelbaum, R. 2018, Annual Review of Astronomy and Astrophysics, 56, 393, doi: 10.1146/annurev-astro-081817-051928
  * Mandelbaum et al. (2012) Mandelbaum, R., Hirata, C. M., Leauthaud, A., Massey, R. J., & Rhodes, J. 2012, MNRAS, 420, 1518, doi: 10.1111/j.1365-2966.2011.20138.x
  * Mandelbaum et al. (2012) Mandelbaum, R., Lackner, C., Leauthaud, A., & Rowe, B. 2012, COSMOS real galaxy dataset, Zenodo, doi: 10.5281/zenodo.3242143
  * Massey et al. (2010) Massey, R., Kitching, T., & Richard, J. 2010, Reports on Progress in Physics, 73, 086901, doi: 10.1088/0034-4885/73/8/086901
  * Massey et al. (2013) Massey, R., Hoekstra, H., Kitching, T., et al. 2013, MNRAS, 429, 661, doi: 10.1093/mnras/sts371
  * Robertson et al. (2019) Robertson, A., Harvey, D., Massey, R., et al. 2019, MNRAS, 488, 3646, doi: 10.1093/mnras/stz1815
  * Romualdez et al. (2019) Romualdez, L. J., Benton, S. J., Brown, A. M., et al. 2019, arXiv preprint arXiv:1911.11210
  * Romualdez et al. (2020) —. 2020, Review of Scientific Instruments, 91, 034501, doi: 10.1063/1.5139711
  * Saito et al. (2020) Saito, S., de la Torre, S., Ilbert, O., et al. 2020, Monthly Notices of the Royal Astronomical Society, 494, 199, doi: 10.1093/mnras/staa727
  * Scoville et al. (2007) Scoville, N., Abraham, R. G., Aussel, H., et al. 2007, ApJS, 172, 38, doi: 10.1086/516580
  * Sérsic (1963) Sérsic, J. L. 1963, Boletin de la Asociacion Argentina de Astronomia La Plata Argentina, 6, 41
  * Weinberg et al. (2013) Weinberg, D. H., Mortonson, M. J., Eisenstein, D. J., et al. 2013, Physics Reports, 530, 87, doi: 10.1016/j.physrep.2013.05.001
  * Zuntz et al. (2018) Zuntz, J., Sheldon, E., Samuroff, S., et al. 2018, Monthly Notices of the Royal Astronomical Society, 481, 1149, doi: 10.1093/mnras/sty2219


Noé Tits1, Zoé Broisson1

# Flowchase: a Mobile Application for Pronunciation Training

###### Abstract

In this paper, we present a solution for providing personalized and instant
feedback to English learners through a mobile application, called Flowchase,
that is connected to a speech technology able to segment and analyze speech
segmental and supra-segmental features. The speech processing pipeline
receives linguistic information corresponding to an utterance to analyze along
with a speech sample. After validation of the speech sample, a joint forced-
alignment and phonetic recognition is performed thanks to a combination of
machine learning models based on speech representation learning that provides
necessary information for designing a feedback on a series of segmental and
supra-segmental pronunciation aspects.

Index Terms: pronunciation training, language learning, speech analysis,
machine learning, transfer learning, human-computer interaction

##  1 Introduction

In the field of Computer-Assisted Language Learning (CALL), there are nowadays
still very few solutions focusing on oral skills, and specifically on
pronunciation. Computer-Assisted Pronunciation Training (CAPT) is an important
research discipline, but there is a lack of concrete applications, although
explicit focus on pronunciation, when combined with the use of technologies,
has a significant impact on L2 learners pronunciation [1, 2]. A reason for
this situation is the gap of complexity between developing feedback on
written, reading or listening skills compared to spoken skills. Indeed for the
first three skill sets, implementing simple heuristics based on multiple
answer exercises, or matching a user answer to a gold standard is
straightforward. On the contrary, providing feedback on spoken skills is not.
A speech technology tailored to analyzing segmental and supra-segmental
patterns is necessary.

The techniques of mispronunciation errors have been close to the findings of
speech recognition area, from HMM-GMM [3], to DNN-HMM [4] and more recently,
transformers [5]. Indeed the tasks share a strong common characteristic, which
is extracting information from audio, a representation of human speech, be it
text or phonetics.

Transfer Learning [6] is today a widely used technique in Deep Learning for
leveraging models trained on related tasks for which there exist abundant
datasets towards tasks for which few data exist. This principle has been
applied successfully for speech technology application [7] with few available
data such as speech recognition for low resource languages, emotion
recognition in speech [8], emotional or expressive speech synthesis [9, 10] or
voice conversion [11], and also to pronunciation assessment [12]. A specific
form of Transfer Learning that was shown to be very efficient is self-
supervised learning where a model is trained to learn representations of input
data without the need for explicit supervision.

In this paper, we present a complete system able to provide a pronunciation
training based on a speech technology built on top of a wav2vec2 [13] model
adapted for mispronunciation detection, integrated in a mobile application.
Although the application contains a mix of tutorials, listening activities and
speaking activities, we focus here on the description of the speaking
activities that involves the speech processing pipeline for analyzing English
learners' pronunciation and providing feedback.

##  2 System

Figure 1 describes the main steps of the user experience inside a speaking
exercise of a learning program. First, the exercise data is shown to the user.
Specifically, it shows an English utterance that the user is expected to say,
with a pronunciation guide to help him understand how it has to be pronounced.
The pronunciation of the sentence can also be heard thanks to a set of
different actor recordings with different English variations. On this screen,
the user can record himself. Then the audio recording is sent to the speech
technology backend along with the exercise information in order to perform
segmentation and analysis of the speech sample. From this analysis, a number
of information are extracted depending on the pronunciation aspect analyzed.

In the second screen, feedback cards are shown to the user in order to
communicate the result of the analysis, and advice in order to improve.

Figure 1: Sequence of steps in a Speaking Activity

Figure 2 details the processing steps happening in the second step explained
above. The speech analysis takes as inputs exercise information such as the
phonetic content and the English learner's speech sample. The user recording
has first to be validated thanks to a series of test on audio that checks is a
valid speech sample, including:

  * •

the duration of the audio is plausible to have a human speech rate compared to
the expected utterance

  * •

the speech sample contains voiced content

  * •

the phonetic content in speech is sufficiently close from the phonetic content

Figure 2: Description of the processing steps of a user recording towards
pronunciation feedback

If the speech sample is validated, a combination of machine learning models
based on speech representation learning is used for performing a forced-
alignment between the speech sample and the phonetic transcription in order to
extract the start and end timings of each phoneme of the sequence. The machine
learning model also analyzes the phonetic content of the audio and allows us
to extract information related to set of different pronunciation aspects such
as analysis of vowels or consonants, and specifically analyzing minimal pairs,
as shown in Figure 3, analysis of intonation such as word stress or sentence
stress, and other supra-segmental aspects like an analysis of pauses between
breath groups in an utterance.

An example of analysis results on a word from a sentence is shown in Figure 3.
Expected phonemes and predicted phonemes are extracted along with the start
and end timings, as well as the respective posterior probabilities according
to the statistical model.

Figure 3: Pronunciation feedback on a word

##  3 Conclusions

In this paper, we presented Flowchase, a mobile application for personalized
pronunciation training that utilizes a speech technology pipeline for
analyzing English learners' pronunciation and providing instant feedback. We
employed transfer learning and self-supervised learning techniques to build a
speech technology model for detecting mispronunciations based on the wav2vec2
architecture.

The system provides feedback on both segmental and supra-segmental aspects of
pronunciation. Our solution addresses the gap in current computer-assisted
language learning applications, which mostly focus on written, reading, or
listening skills. Flowchase provides a tool for improving oral language
skills, particularly pronunciation, which is crucial for effective
communication. Future work includes testing the effectiveness of the
application and the speech technology pipeline in real-world settings and
extending the system to support other languages.

##  4 Acknowledgements

This work is part of the project REDCALL that is partially funded by a FIRST
Entreprise Docteur program from SPW Recherche111https://recherche.wallonie.be/

## References

  * [1] N. Cordier and F. Meunier, ``Digital tools to improve the english pronunciation of l2 learners: a case study with flowchase.''
  * [2] X. Anguera and V. Van, ``English language speech assistant.'' in _INTERSPEECH_ , 2016, pp. 1962–1963.
  * [3] S. M. Witt and S. J. Young, ``Phone-level pronunciation scoring and assessment for interactive language learning,'' _Speech communication_ , vol. 30, no. 2-3, pp. 95–108, 2000.
  * [4] W. Hu, Y. Qian, F. K. Soong, and Y. Wang, ``Improved mispronunciation detection with deep neural network trained acoustic models and transfer learning based logistic regression classifiers,'' _Speech Communication_ , vol. 67, pp. 154–166, 2015.
  * [5] Y. Gong, Z. Chen, I.-H. Chu, P. Chang, and J. Glass, ``Transformer-based multi-aspect multi-granularity non-native english speaker pronunciation assessment,'' in _ICASSP 2022-2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2022, pp. 7262–7266.
  * [6] C. Tan, F. Sun, T. Kong, W. Zhang, C. Yang, and C. Liu, ``A survey on deep transfer learning,'' in _International conference on artificial neural networks_. Springer, 2018, pp. 270–279.
  * [7] D. Wang and T. F. Zheng, ``Transfer learning for speech and language processing,'' in _2015 Asia-Pacific Signal and Information Processing Association Annual Summit and Conference (APSIPA)_. IEEE, 2015, pp. 1225–1237.
  * [8] N. Tits, K. El Haddad, and T. Dutoit, ``Asr-based features for emotion recognition: A transfer learning approach,'' in _Proceedings of Grand Challenge and Workshop on Human Multimodal Language (Challenge-HML)_. Association for Computational Linguistics, 2018, pp. 48–52. [Online]. Available: http://aclweb.org/anthology/W18-3307
  * [9] ——, ``Exploring Transfer Learning for Low Resource Emotional TTS,'' in _Intelligent Systems and Applications_ , Y. Bi, R. Bhatia, and S. Kapoor, Eds. Cham: Springer International Publishing, 2020, pp. 52–60.
  * [10] N. Tits, F. Wang, K. E. Haddad, V. Pagel, and T. Dutoit, ``Visualization and Interpretation of Latent Spaces for Controlling Expressive Speech Synthesis through Audio Analysis,'' in _Proc. Interspeech 2019_ , 2019, pp. 4475–4479. [Online]. Available: http://dx.doi.org/10.21437/Interspeech.2019-1426
  * [11] K. Zhou, B. Sisman, R. Liu, and H. Li, ``Emotional voice conversion: Theory, databases and esd,'' _Speech Communication_ , vol. 137, pp. 1–18, 2022.
  * [12] B. Lin and L. Wang, ``Deep feature transfer learning for automatic pronunciation assessment.'' in _Interspeech_ , 2021, pp. 4438–4442.
  * [13] A. Baevski, Y. Zhou, A. Mohamed, and M. Auli, ``wav2vec 2.0: A framework for self-supervised learning of speech representations,'' _Advances in neural information processing systems_ , vol. 33, pp. 12 449–12 460, 2020.


# Time evolution of a mean-field generalized contact process

by Logan Chariker1 and Joel L. Lebowitz2

1School of Natural Sciences, Institute for Advanced Study, Princeton, New
Jersey.
2Departments of Mathematics and Physics, Rutgers University, New Brunswick,
New Jersey.

Keywords: generalized contact process, neurons with discrete voltage, mean
field, exact solution, spatial dependence

Abstract. We investigate the macroscopic time evolution and stationary states
of a mean field generalized contact process in $\mathbb{R}^{d}$. The model is
described by a coupled set of nonlinear integral-differential equations. It
was inspired by a model of neurons with discrete voltages evolving by a
stochastic integrate and fire mechanism. We obtain a complete solution in the
spatially uniform case and partial solutions in the general case. The system
has one or more fixed points and also traveling wave solutions.

##  1\. Introduction

We consider the mean-field (hydrodynamic) limit of a (novel) stochastic
lattice system inspired by neuronal integrate-and-fire models [7]. On the
microscopic level the system consists of variables $S(\boldsymbol{z},t)$
associated to lattice site $\boldsymbol{z}\in\Omega\subset\mathbb{Z}^{d}$, at
time $t\in\mathbb{R}_{+}$. $S(\boldsymbol{z},t)$ can take integer values
$0,1,\ldots,k$. These can be thought of as discrete values of the voltage of a
neuron at site $\boldsymbol{z}$, or as the state of infection of the
individual at $\boldsymbol{z}$. There are $k$ stages of the infection with
$S(\boldsymbol{z},t)=0$ corresponding to the healthy state.

A sketch of the microscopic dynamics is as follows: when none of the
$S(\boldsymbol{z},t)$ are in state $k$, the system is in a static state. When
$S(\boldsymbol{z},t)=k$ it has a probability $dt$ of firing (healing, dying)
during the time interval $(t,t+dt)$. This is independent of the values
$S(\boldsymbol{w},t)$ for $\boldsymbol{w}\neq\boldsymbol{z}$. When
$S(\boldsymbol{z},t)$ fires it is instantaneously reset to
$S^{\prime}(\boldsymbol{z},t)=0$. A neuron at site $\boldsymbol{w}$, with
$S(\boldsymbol{w},t)=j$, $j\leq k-1$, will jump from $j$ to $j+1$ with rate
$k\lambda J_{\gamma}(\boldsymbol{z},\boldsymbol{w})$, $\lambda>0$, when
$S(\boldsymbol{z},t)=k$. The function
$J_{\gamma}(\boldsymbol{x},\boldsymbol{y})$ has the Kac form [8],

(1) $\displaystyle
J_{\gamma}(\boldsymbol{z},\boldsymbol{w})=\gamma^{d}J(\gamma(\boldsymbol{z}-\boldsymbol{w}))=\gamma^{d}J(\gamma(\boldsymbol{w}-\boldsymbol{z}))\geq
0$

with $\int_{\mathbb{R}^{d}}J(\boldsymbol{r})d\boldsymbol{r}=1$.

When $k=1$, each site can exist in one of only two states, one an inactive
state and the other active, and in this respect the model becomes similar to a
contact process and also to a popular two-state neuron system known as the
stochastic Wilson-Cowan model [9, 15]. The case $k>1$ introduces multiple
inactive states $0,\ldots,k-1$, in which a site exerts no influence on its
neighbors, and it must traverse sequentially through the states to reach the
active state $k$, which is followed by a reset at rate 1 to the state $0$.
This is inspired by integrate and fire neuron models, where the firing of
other neurons is required to drive the membrane potential of a particular
neuron from a resting potential to a threshold level in order for it to become
active itself. There is an abundance of integrate and fire neuron models used
throughout neuroscience. Less tractable than two-state neuron models, many
results for integrate and fire models are given by numerical simulation (e.g.,
[11], [3], [12]), although there are some results proved, in particular for
non-spatially-dependent networks in the mean field limit[2].

We are interested here in the inclusion of a spatial structure. [5] and [6]
have shown rigorous results in the hydrodynamic limit of a network of neurons,
which have continuous membrane potentials and a stochastic firing threshold.
Individual neuron firings produce infinitesimal jumps in the membrane
potential of all other neurons, weighted by a spatially dependent connectivity
function. In contrast, neurons in the model studied in this paper have a fixed
firing threshold.

We will be interested in the macroscopic equations on the scale $x=\gamma z$
in the mean field limit $\gamma\rightarrow 0$. We shall not discuss here the
derivation of the macroscopic equations from the microscopic models. That will
be done in a different publication [4]. Here we shall discuss the solution of
the resulting macroscopic equations for different values of $k$.

##  2\. Macroscopic Equations

The $\gamma\rightarrow 0$ limit of the microscopic model yields the following
equations on the macroscopic spatial scale for the $v_{j}(\boldsymbol{x},t)$,
the fraction of the population density at position $x$ in state $j$,
$j=\\{0,1,\ldots,k\\}$,

(2) $\displaystyle\frac{\partial v_{0}(\boldsymbol{x},t)}{\partial t}$
$\displaystyle=v_{k}(\boldsymbol{x},t)-v_{0}(\boldsymbol{x},t)\lambda
kR_{k}(\boldsymbol{x},t)$ (3) $\displaystyle\frac{\partial
v_{j}(\boldsymbol{x},t)}{\partial t}$
$\displaystyle=[v_{j-1}(\boldsymbol{x},t)-v_{j}(\boldsymbol{x},t)]\lambda
kR_{k}(\boldsymbol{x},t),\quad j=1,\ldots,k-1$ (4)
$\displaystyle\frac{\partial v_{k}(\boldsymbol{x},t)}{\partial t}$
$\displaystyle=-v_{k}(\boldsymbol{x},t)+v_{k-1}(\boldsymbol{x},t)\lambda
kR_{k}(\boldsymbol{x},t),$

where $\boldsymbol{x}\in\Lambda\subset\mathbb{R}^{d}$ is a cubical box of
sides $L$, with periodic boundary conditions and uniform density equal to one,

(5) $\displaystyle
R_{k}(\boldsymbol{x},t)=\int_{\Lambda}J(\boldsymbol{x}-\boldsymbol{y})v_{k}(\boldsymbol{y},t)d\boldsymbol{y}.$

We assume that $J(r)$ has a range less than $L/2$. In the spatially uniform
state, when $v_{j}(\boldsymbol{y},t)$ is independent of $\boldsymbol{y}$,
then, by (5) and (1), $R_{k}=v_{k}$.

It follows from equations (2)-(4) that starting with
$v_{j}(\boldsymbol{x},0)\geq 0$, $\sum_{j=0}^{k}v_{j}(\boldsymbol{x},0)=1$,
then

(6)
$\displaystyle\sum_{i=0}^{k}v_{i}(\boldsymbol{x},t)=\sum_{i=0}^{k}v_{i}(\boldsymbol{x},0)=1,\quad
v_{i}(\boldsymbol{x},t)\geq 0,\qquad\text{for all $t\geq 0$}.$

It is clear from the above equations that if $v_{k}(\boldsymbol{x},0)=0$ for
all $\boldsymbol{x}$ then $v_{i}(\boldsymbol{x},t)=v_{i}(\boldsymbol{x},0)$
and the system remains in its initial state forever. This is not very
interesting and we shall assume from now on that $v_{k}(\boldsymbol{x},0)>0$
for some values of $\boldsymbol{x}$. It is still possible however that
$v_{k}(\boldsymbol{x},t)\rightarrow 0$ as $t\rightarrow\infty$. We shall call
that “extinction”.

Using (6) we can replace $v_{k}(x,t)$ in (2) and (3) by

$v_{k}(x,t)=1-\sum_{j=0}^{k-1}v_{j}(x,t)$

leading to a closed set of equations for $v_{j}(x,t)$, $j=0,\ldots,k-1$ with

$R_{k}(x,t)=\int_{\Lambda}J(x-y)\Big{[}1-\sum_{j=0}^{k-1}v_{j}(y,t)\Big{]}dy.$

##  3\. Stationary States

Consider now the stationary solutions of (2)-(4). We see that in addition to
the $R_{k}=0$ solution corresponding to $v_{k}=0$, there is a stationary
solution $R_{k}(\boldsymbol{x},t)=\bar{R}(\boldsymbol{x})>0$ of the form

(7) $\displaystyle\bar{v}_{j}(\boldsymbol{x})=\frac{1}{\lambda
k\bar{R}(\boldsymbol{x})+k},\quad j\leq k-1$

while

(8) $\displaystyle\bar{v}_{k}(\boldsymbol{x})$
$\displaystyle=\frac{\lambda\bar{R}(\boldsymbol{x})}{1+\lambda\bar{R}(\boldsymbol{x})},$

and $\bar{R}(\boldsymbol{x})$ satisfies the eq.

(9) $\displaystyle\bar{R}(\boldsymbol{x})$ $\displaystyle=\lambda\int
J(\boldsymbol{x}-\boldsymbol{y})\frac{\bar{R}(\boldsymbol{y})}{\lambda\bar{R}(\boldsymbol{y})+1}d\boldsymbol{y}.$

Equations (7)-(9) have the spatially uniform solution

(10) $\displaystyle\bar{R}=\bar{v}_{k}=\frac{\lambda-1}{\lambda}\quad\text{
and }\quad\bar{v}_{j}=\frac{1}{\lambda k},\quad j<k,$

which is a physical sustaining stationary solution with $\bar{v}_{k}>0$, for
$\lambda>1$. For $\lambda\leq 1$ the only physical stationary solution is
$\bar{v}_{k}=0$.

Note that equation (9) for the stationary $\bar{R}(\boldsymbol{x})$ is
independent of $k$. Thus if there exists a non-vanishing, spatially dependent
stationary state $\bar{R}(\boldsymbol{x})$, then it will be so for all $k$. We
will show in section 6 that for $k=1$, the only stationary nonzero $v_{k}$ is
the spatially uniform one. Hence this will be true for all $k$.

The stationary state $\bar{\boldsymbol{v}}=(\bar{v}_{0},\ldots,\bar{v}_{k})$
given in (10) is linearly stable (as shown in section 5), so if we start close
enough to $\bar{\boldsymbol{v}}$ then the system will always approach
$\bar{\boldsymbol{v}}$ as $t\rightarrow\infty$. On the other hand there are,
as shown in section 4 for $k>1$, initial uniform states $v_{j}(0)$ with
$v_{k}(0)>0$ such that $v_{k}(t)\rightarrow 0$ as $t\rightarrow\infty$.
Initial states $v_{j}(\boldsymbol{x},0)$ close to $v_{j}(0)$ would also have
$v_{k}(\boldsymbol{x},t)\rightarrow 0$ as $t\rightarrow\infty$ for all
$\boldsymbol{x}$; see figure 1 and section 7.

##  4\. Solution of the macroscopic equations in the spatially uniform case

The macroscopic equations in the spatially uniform case take the form

(11) $\displaystyle\frac{dv_{0}}{dt}$ $\displaystyle=v_{k}-(k\lambda
v_{k})v_{0}$ (12) $\displaystyle\frac{dv_{j}}{dt}$ $\displaystyle=(k\lambda
v_{k})[v_{j-1}-v_{j}],\qquad j=1,\ldots,k-1$ (13)
$\displaystyle\frac{dv_{k}}{dt}$ $\displaystyle=-v_{k}+(k\lambda
v_{k})v_{k-1}$

where $v_{j}(t)=$ fraction of neurons in state $j=0,1,\ldots,k$.

To simplify equations (11)-(13), we introduce the variable $r(t)$ defined by
the equations

$\frac{dr}{dt}=k\lambda v_{k},\qquad r(0)=0.$

Writing $v_{j}(t)=\tilde{v}_{j}(r(t))$, valid as long as $v_{k}>0$, we get

(14) $\displaystyle\frac{d\tilde{v}_{0}}{dr}$
$\displaystyle=-\tilde{v}_{0}+q_{0},\qquad q_{0}=\frac{1}{k\lambda},$ (15)
$\displaystyle\frac{d\tilde{v}_{j}}{dr}$
$\displaystyle=\tilde{v}_{j-1}-\tilde{v}_{j},\qquad j=1,\ldots,k-1,$ (16)
$\displaystyle\frac{d\tilde{v}_{k}}{dr}$
$\displaystyle=\tilde{v}_{k-1}+q_{k},\qquad q_{k}=-q_{0}=-\frac{1}{\lambda
k}.$

The autonomous equations for $j=0,1,\ldots,k-1$ can be rewritten in vector
form

(17) $\displaystyle\frac{d\tilde{\boldsymbol{v}}}{dr}$
$\displaystyle=\mathbf{A}\tilde{\boldsymbol{v}}+\boldsymbol{q},\qquad\tilde{\boldsymbol{v}}(0)=\boldsymbol{v}(0)$

where
$\tilde{\boldsymbol{v}}=(\tilde{v}_{0},\tilde{v}_{1},\ldots,\tilde{v}_{k-1})^{T}$,
$\boldsymbol{q}=(q_{0},0,\ldots,0)^{T}$, and $\mathbf{A}$ is a $k\times k$
square matrix with $i,j\in\\{0,1,\ldots,k-1\\}$. $A$ has $-1$’s along the
diagonal and $+1$’s along the first subdiagonal:

(18)
$\displaystyle\boldsymbol{\boldsymbol{A}}=-\boldsymbol{\boldsymbol{I}}+\boldsymbol{\boldsymbol{B}}_{k},\qquad\boldsymbol{\boldsymbol{B}}_{k}=\begin{pmatrix}0&&&&&\\\
1&0&&&&\\\ &1&&&&\\\ &&\ddots&&&\\\ &&&1&0&\\\ &&&&1&0\\\ \end{pmatrix}.$

Note that in successive powers of $\mathbf{B}_{k}$, the $+1$’s move to lower
subdiagonals:

$(\boldsymbol{B}_{k})_{i,j}=\delta_{i-1,j},(\boldsymbol{B}_{k}^{2})_{i,j}=\delta_{i-2,j},\ldots,(\boldsymbol{B}_{k}^{k-1})_{i,j}=\delta_{i-k+1,j},\text{
and }\boldsymbol{B}_{k}^{l}=\boldsymbol{0}\text{ for $l\geq k$}.$

Equation (17) has the solution

(19)
$\displaystyle\tilde{\boldsymbol{v}}(r)=e^{\mathbf{A}r}\tilde{\boldsymbol{v}}(0)+\int_{0}^{r}e^{s\mathbf{A}}\boldsymbol{q}ds.$

Using equation (18), and the fact that $\mathbf{B}_{k}^{k}=\mathbf{0}$, we
have that

$e^{Ar}=e^{-r}\mathbf{1}e^{-\mathbf{B}r}=e^{-r}\left[\mathbf{1}+r\mathbf{B}_{k}+\frac{r^{2}}{2}\mathbf{B}^{2}_{k}+\cdots+\frac{r^{k-1}}{(k-1)!}\mathbf{B}^{k-1}_{k}\right].$

Each term within the square brackets corresponds to a distinct subdiagonal, so
the matrix $e^{Ar}$ is 0 above the diagonal and constant along each
subdiagonal. Let $H_{j}$ be the value in the $j$th subdiagonal: explicitly,

$H_{j}(r)=e^{-r}\frac{r^{j}}{j!},$

for $j=0,\ldots,k-1$. Then (19) can be expanded to give an explicit solution
for $\tilde{v}_{j}$ in terms of $r$.

(20)
$\displaystyle\tilde{v}_{j}(r)=\sum_{i=0}^{j}H_{j-i}\tilde{v}_{i}(0)+\frac{1}{k\lambda}\int_{0}^{r}H_{j}(s)ds,\qquad
j=0,\ldots,k-1,$

We also have

(21) $\displaystyle\tilde{v}_{k}(r)$
$\displaystyle=1-\sum_{j=0}^{k-1}\tilde{v}_{j}(r)$
$\displaystyle=1-\sum_{j=0}^{k-1}\sum_{i=0}^{j}H_{j-i}(r)\tilde{v}_{i}(0)-\frac{1}{k\lambda}\int_{0}^{r}\sum_{j=0}^{k-1}H_{j}(s)ds.$

Clearly $\tilde{v}_{i}\rightarrow\bar{v}_{i}=(\lambda k)^{-1}$ for $i<k$, and
$\tilde{v}_{k}\rightarrow\bar{v}_{k}=(\lambda-1)/\lambda$ as
$r\rightarrow\infty$. Recalling now that $dr/dt=k\lambda v_{k}$, we get

(22) $\displaystyle\frac{dr}{dt}$
$\displaystyle=k\Big{(}\lambda-1-\lambda\sum_{j=0}^{k-1}\sum_{i=0}^{j}H_{j-i}(r)\Big{[}v_{i}(0)-\frac{1}{\lambda
k}\Big{]}\Big{)}=\phi(r),$

which is an autonomous ODE, emphasized by introducing the notation $\phi(r)$
on the right. The behavior of $r(t)$ can then be determined by analyzing
$\phi(r)$.

Starting with $v_{k}(0)>0$, we see that $\phi(0)>0$, and so $r(t)$ is monotone
increasing, and satisfies

(23) $\displaystyle t=\int_{0}^{r}\frac{ds}{\phi(s)}$

so we have either

  1. 1.

there exists $r_{0}$ the smallest positive solution to $\phi(r)=0$, such that
$t\rightarrow\infty$ as $r\rightarrow r_{0}$, or

  2. 2.

the integral is finite, in which case $r\rightarrow\infty$ as
$t\rightarrow\infty$. This will certainly be the case if $\phi(r)>0$ for all
$r>0$.

In case 1, since $\phi(r)=dr/dt=k\lambda v_{k}$, we see that $v_{k}\rightarrow
0$ as $t\rightarrow\infty$; that is, in this case the system goes to an
inactive state with the firings dying out. In case 2, we see by equation (22)
that as $t\rightarrow\infty$,

$v_{k}(t)\rightarrow\frac{\phi(\infty)}{\lambda k}=1-\frac{k}{\lambda
k}=\bar{v}_{k},$

and by equation (20) for $j<k$, that

$v_{j}(t)\rightarrow\tilde{v}_{j}(\infty)=\frac{1}{k\lambda}=\bar{v}_{j},$

so the system goes to the unique sustaining stationary state.

Illustrative Examples

The $\mathbf{k=1}$ case. In this case,

(24) $\displaystyle\frac{dr}{dt}=\lambda
v_{1}(t)=\lambda\left(1-e^{-r}v_{0}(0)-\frac{1}{\lambda}\int_{0}^{r}e^{-s}ds\right)=\phi(r).$

If $0\leq v_{0}(0)<1$, then it can be checked that the right-hand side of
equation (24) is bounded below by a positive constant for all $r\geq 0$.
Therefore the system always approaches the sustaining steady state solution if
$v_{1}(0)>0$. More explicitly, we have in this case, $v_{0}(t)=1-v_{1}(t)$, so
$v_{1}(t)$ satisfies the autonomous equation

(25) $\displaystyle\frac{dv_{1}(t)}{dt}=(\lambda-1)v_{1}(t)-\lambda
v_{1}^{2}(t),$

whose solution is

(26) $\displaystyle
v_{1}(t)=\frac{v_{1}(0)}{\frac{\lambda}{\lambda-1}v_{1}(0)+\big{(}1-\frac{\lambda}{\lambda-1}v_{1}(0)\big{)}e^{-(\lambda-1)t}}\rightarrow\bar{v}_{1}=\frac{\lambda-1}{\lambda}\text{
as $t\rightarrow\infty$}.$

Thus for $k=1$ any initial state with $v_{1}(0)\neq 0$ will approach, as
$t\rightarrow\infty$, the sustaining stationary state exponentially, as long
as $\lambda>1$.

Thinking of the process as a mean field model of infection, with $v_{1}$
representing the infected fraction of the population the model predicts as in
the standard contact process a persistent percentage of infected individuals
for $\lambda>1$, the percentage increasing with $\lambda$. For $\lambda<1$,
$v_{1}(t)\rightarrow 0$, there is no epidemic, as everyone gets eventually
cured.

The $\mathbf{k=2}$ case. Unlike the $k=1$ case, here it is possible to start
the system with $v_{2}(0)>0$ and still have the firing die out,
$v_{2}(t)\rightarrow 0$ as $t\rightarrow\infty$.

Writing the solution (20) for $k=2$, we get,

$\displaystyle\tilde{v}_{0}(r)+\tilde{v}_{1}(r)$
$\displaystyle=\frac{1}{\lambda}[1-e^{-r}]-\frac{1}{2\lambda}re^{-r}+v_{0}(0)e^{-r}[1+r]+v_{1}(0)e^{-r}.$

This yields

(27) $\displaystyle\phi(r)$
$\displaystyle=2\lambda[1-\tilde{v}_{1}(r)-\tilde{v}_{0}(r)]$
$\displaystyle=2(\lambda-1)(1-e^{-r})+re^{-r}[1-2\lambda v_{0}(0)]+2\lambda
e^{-r}[1-v_{1}(0)-v_{0}(0)].$

Clearly if $v_{0}(0)<\frac{1}{2\lambda}$ then $\phi(r)>0$ for all $r>0$ and
the system will go to the sustaining stationary state. For $v_{0}(0)$ close to
1, there exists $r_{0}>0$ for which $\phi(r_{0})=0$, in which case the system
goes to an inactive state with the firing dying out, i.e., $v_{2}\rightarrow
0$ as $t\rightarrow\infty$. The region in the $v_{0}(0)$, $v_{1}(0)$ plane for
which this occurs shrinks as $\lambda$ increases, as demonstrated numerically
in figure 1. We expect similar behavior for $k>2$. In fact it is easy to see
from (22) that if $v_{i}(0)\leq\frac{1}{\lambda k}$ for all $i\leq k-1$ then
$\phi(r)>0$.

Figure 1. Solutions of (27) for the $k=2$ spatially uniform case, indicating
in black the initial conditions $(v_{0}(0),v_{1}(0))$ which go to the extinct
states $v_{2}=0$ as $t\rightarrow\infty$. The extinct states are indicated by
the dashed line. The sustaining stationary state
$\bar{v}_{0}=\bar{v}_{1}=1/(2\lambda)$ is shown by a red dot: it attracts all
initial states in the white region below the dotted line.

##  5\. Linear Stability of the Sustaining State

Let us consider now the linear stability of the sustaining stationary states
for general $k$

$\displaystyle\bar{v}_{j}=\frac{1}{\lambda k},\qquad
j<k,\qquad\bar{v}_{k}=1-\frac{1}{\lambda}.$

Let $v_{j}(\boldsymbol{x},0)=\bar{v}_{j}+f_{j}(\boldsymbol{x},t)$. Then
setting

$v_{j}(\boldsymbol{x},t)=\bar{v}_{j}+f_{j}(\boldsymbol{x},t),$

The conditions

$\displaystyle\sum_{j=0}^{k}f_{j}(\boldsymbol{x},t)=0,$ $\displaystyle
1-\bar{v}_{j}\geq f_{j}(\boldsymbol{x},t)\geq-\bar{v}_{j}$

are satisfied for all $t\geq 0$ if they are satisfied for $t=0$.

Setting $(\lambda-1)k=\alpha$, we obtain to first order in the $f_{j}$,

(28) $\displaystyle\frac{\partial f_{0}}{\partial t}(\boldsymbol{x},t)$
$\displaystyle=-\alpha f_{0}(\boldsymbol{x},t)+\int
J(\boldsymbol{x}-\boldsymbol{y})[f_{k}(\boldsymbol{x})-f_{k}(\boldsymbol{y})]d\boldsymbol{y}$
(29) $\displaystyle\frac{\partial f_{j}}{\partial t}(\boldsymbol{x},t)$
$\displaystyle=\alpha[f_{j-1}(\boldsymbol{x},t)-f_{j}(\boldsymbol{x},t)],\qquad
j=1,\ldots,k-1.$

Taking the Fourier series in $\boldsymbol{x}$,

(30)
$\displaystyle\hat{f}_{j}(\boldsymbol{\xi},t)=\frac{1}{L^{d}}\int_{\Lambda}f_{j}(\boldsymbol{x},t)e^{-2\pi
i\boldsymbol{\xi}\cdot\boldsymbol{x}}d\boldsymbol{x},\qquad j=0,\ldots,k-1,$

where $\boldsymbol{\xi}=(\xi_{1},\ldots,\xi_{d})\in(L^{-1}\mathbb{Z})^{d}$,
gives

(31) $\displaystyle\frac{\partial\hat{f}_{0}}{\partial t}(\boldsymbol{\xi},t)$
$\displaystyle=-\alpha\hat{f_{0}}(\boldsymbol{\xi},t)+\hat{f}_{k}(\boldsymbol{\xi},t)\big{(}1-\hat{J}(\boldsymbol{\xi})\big{)}$
(32) $\displaystyle\frac{\partial\hat{f}_{j}}{\partial t}(\boldsymbol{\xi},t)$
$\displaystyle=\alpha[\hat{f}_{j-1}(\boldsymbol{\xi},t)-\hat{f}_{j}(\boldsymbol{\xi},t)],\qquad
j=1,\ldots,k-1.$

Letting $\beta(\boldsymbol{\xi})=(1-\hat{J}(\boldsymbol{\xi}))/\alpha$, we can
write the above equations in vector form

(33) $\displaystyle\frac{\partial\hat{\boldsymbol{f}}}{\partial
t}(\boldsymbol{\xi},t)=\alpha(\boldsymbol{A}-\beta(\boldsymbol{\xi})\boldsymbol{M})\hat{\boldsymbol{f}},\qquad\boldsymbol{M}=\begin{pmatrix}1&1&\cdots&1\\\
0&0&\cdots&0\\\ \vdots&&\ddots&\vdots\\\ 0&0&\cdots&0\end{pmatrix},$

where $\boldsymbol{M}$ is the $k\times k$ matrix with ones in the top row and
zeros elsewhere.

We will show that for any fixed value of $\boldsymbol{\xi}$, each eigenvalue
of the matrix $(\boldsymbol{A}-\beta\boldsymbol{M})$ has a negative real part,
implying the convergence of $\hat{\boldsymbol{f}}(\boldsymbol{\xi},t)$ to
$\boldsymbol{0}$ as $t\rightarrow\infty$. The characteristic polynomial of the
matrix, computed by cofactor expansion along the top row, is

(34) $\displaystyle p(x)$
$\displaystyle=(-1-\beta-x)-\sum_{k=2}^{k}(-1)^{j}(-\beta)(+1)^{j-1}(-1-x)^{k-j}$
$\displaystyle=(-1)^{k}\Big{[}(1+\beta+x)(1+x)^{k-1}+\beta\sum_{j=2}^{k}(1+x)^{k-j}\Big{]}.$

Letting $p(x)=0$, $y=x+1$, and noting that $y=1$ is not a solution to the
equation, yields

(35) $\displaystyle 0$ $\displaystyle=(y+b)y^{k-1}+\beta\sum_{j=2}^{k}y^{k-j}$
$\displaystyle=y^{k}+\beta\sum_{j=1}^{k}y^{k-j}$
$\displaystyle=y^{k}+\beta\frac{1-y^{k}}{1-y},$

which simplifies to

(36) $\displaystyle y^{k}(y-(1-\beta))=\beta.$

We now note that $\beta\geq 0$. This follows from the fact that $J$ is a real-
valued function with $J(\boldsymbol{x})=J(-\boldsymbol{x})\geq 0$ and
$\int_{\Lambda}J(\boldsymbol{x})d\boldsymbol{x}=1$. Therefore applying the
absolute value to both sides of (36) gives

(37) $\displaystyle|y|^{k}|y-(1-\beta)|=\beta.$

We can deduce that $|y|\neq 1$, as otherwise the equation (37) implies $y=1$,
and we have remarked that this not a solution to the characteristic equation.
Therefore either

  1. 1.

$|y|<1$, or

  2. 2.

$|y|>1$.

In case 1, it follows easily that the real part of $y$ is less than 1, and
therefore the real part of $x$ is negative. Case 2 implies that
$|y-(1-\beta)|<\beta$, from which it also follows that the real part of $y$ is
less than 1, and that the real part of $x$ is negative. This implies that
$\hat{f}_{j}(\boldsymbol{\xi},t)\rightarrow 0$ for all fixed $\xi$, and from
this point, we can prove that $\boldsymbol{f}$ converges to $\boldsymbol{0}$
uniformly.

Linear stability of the inert state $\mathbf{v_{0}=1}$. Having shown the
linear stability of the stationary sustaining state we consider now the linear
stability of the state $v_{k}(t)=v_{k}(0)=0$. Rather than considering all
initial states $v_{j}(0)$, $\sum_{j=0}^{k-1}v_{j}(0)=1$, we consider here only
perturbations around the extreme case $v_{0}(\boldsymbol{x},t)=1$,
$v_{j}(\boldsymbol{x},t)=0$, $j>0$.

Let

(38) $\displaystyle v_{0}(\boldsymbol{x},t)$
$\displaystyle=1+f_{0}(\boldsymbol{x},t),\qquad f_{0}(\boldsymbol{x},t)<0,$
(39) $\displaystyle v_{j}(\boldsymbol{x},t)$
$\displaystyle=f_{j}(\boldsymbol{x},t),\qquad f_{j}(\boldsymbol{x},t)>0,\qquad
j=1,\ldots,k,\qquad\sum_{j=0}^{k}f_{j}=0.$

Linearizing in the $f$’s gives, for $k>1$,

(40) $\displaystyle\frac{\partial f_{0}}{\partial t}(\boldsymbol{x},t)$
$\displaystyle=f_{k}(\boldsymbol{x},t)-\lambda k\int
J(\boldsymbol{x}-\boldsymbol{y})f_{k}(\boldsymbol{y},t)d\boldsymbol{y}$ (41)
$\displaystyle\frac{\partial f_{1}}{\partial t}(\boldsymbol{x},t)$
$\displaystyle=\lambda k\int
J(\boldsymbol{x}-\boldsymbol{y})f_{k}(\boldsymbol{y},t)d\boldsymbol{y}$ (42)
$\displaystyle\frac{\partial f_{j}}{\partial t}(\boldsymbol{x},t)$
$\displaystyle=0,\quad\text{for $j\in\\{2,\ldots,k-1\\}$},$ (43)
$\displaystyle\frac{\partial f_{k}}{\partial t}(\boldsymbol{x},t)$
$\displaystyle=-f_{k}=\sum_{j=0}^{k-1}f_{j},\quad\text{ implying
}f_{k}(\boldsymbol{x},t)=f_{k}(\boldsymbol{x},0)e^{-t},$

so for $k>1$ the dead state in the vicinity of $v_{0}(\boldsymbol{x},t)=1$ is
linearly stable, see figure 1.

For $k=1$, $f_{1}(\boldsymbol{x},t)\geq 0$,

(44) $\displaystyle\frac{\partial f_{1}}{\partial
t}(\boldsymbol{x},t)=-f_{1}(\boldsymbol{x},t)+\lambda\int
J(\boldsymbol{x}-\boldsymbol{y})f_{1}(\boldsymbol{y},t)d\boldsymbol{y}.$

Taking spatial Fourier transforms yields

(45) $\displaystyle\frac{\partial\hat{f}_{1}}{\partial t}(\boldsymbol{\xi},t)$
$\displaystyle=(\lambda\hat{J}(\boldsymbol{\xi})-1)\hat{f}_{1}(\boldsymbol{\xi},t),$

which has the solution

(46) $\displaystyle\hat{f}_{1}(\boldsymbol{\xi},t)$
$\displaystyle=\hat{f}_{1}(\boldsymbol{\xi},0)e^{(\lambda\hat{J}(\boldsymbol{\xi})-1)t}.$

Since $\hat{J}(0)=1$ there will be growth at least for small values of
$\boldsymbol{\xi}$, for which $\hat{f}_{1}(\boldsymbol{\xi},0)>0$, so the
state $v_{0}=1$ is unstable for $\lambda>1$. In fact as we shall now show for
$k=1$, any perturbation of the state $v_{0}(0)=1$, will lead asymptotically to
the stable stationary state $\bar{v}_{0}=\lambda^{-1}$,
$\bar{v}_{1}=(\lambda-1)/\lambda.$

##  6\. $k=1$, general case.

In this two-level case, $v_{1}(\boldsymbol{x},t)=1-v_{0}(\boldsymbol{x},t)$ is
the only unknown function. It satisfies the equation

(47) $\displaystyle\frac{\partial v_{1}(\boldsymbol{x},t)}{\partial
t}=-v_{1}+\lambda(1-v_{1}(\boldsymbol{x},t))\int
J(\boldsymbol{x}-\boldsymbol{y})v_{1}(\boldsymbol{y},t)d\boldsymbol{y}.$

Define $f(\boldsymbol{x},t)$ by

(48) $\displaystyle
f(\boldsymbol{x},t)=v_{1}(\boldsymbol{x},t)-\bar{v}_{1},\qquad\bar{v}_{1}=\frac{\lambda-1}{\lambda}.$

Using the fact that $\int J(\boldsymbol{x}-\boldsymbol{y})d\boldsymbol{y}=\int
J(\boldsymbol{x}-\boldsymbol{y})d\boldsymbol{x}=1$ we get from (47)

(49) $\displaystyle\frac{\partial f(\boldsymbol{x},t)}{\partial t}=$
$\displaystyle-(\lambda-1)f(\boldsymbol{x})$
$\displaystyle+\int_{\Lambda}J(\boldsymbol{x}-\boldsymbol{y})[f(\boldsymbol{y})-f(\boldsymbol{x})]d\boldsymbol{y}$
$\displaystyle-\lambda\int_{\Lambda}J(\boldsymbol{x}-\boldsymbol{y})f(\boldsymbol{y})f(\boldsymbol{x})d\boldsymbol{y}.$
$\displaystyle=-\int
d\boldsymbol{y}J(\boldsymbol{x}-\boldsymbol{y})[f(\boldsymbol{x})-f(\boldsymbol{y})]-\lambda\int
d\boldsymbol{y}J(\boldsymbol{x}-\boldsymbol{y})f(x)\left[\frac{\lambda-1}{\lambda}+f(\boldsymbol{y})\right]$

Multiplying (49) by $f(\boldsymbol{x},t)$ and integrating over
$\boldsymbol{x}$ yields

(50)
$\displaystyle\frac{1}{2}\frac{d}{dt}\int_{\Lambda}f^{2}(\boldsymbol{x},t)d\boldsymbol{x}=$
$\displaystyle-\frac{1}{2}\iint
J(\boldsymbol{x}-\boldsymbol{y})[f(\boldsymbol{x})-f(\boldsymbol{y})]^{2}d\boldsymbol{x}d\boldsymbol{y}$
$\displaystyle-\lambda\iint
J(\boldsymbol{x}-\boldsymbol{y})f^{2}(\boldsymbol{x},t)v_{1}(\boldsymbol{y},t)d\boldsymbol{y}d\boldsymbol{x}\leq
0$

The inequality is strict for all initial conditions with
$v_{1}(\boldsymbol{x},0)$ not identically 0, and shows that
$v_{1}(\boldsymbol{x},t)\rightarrow\bar{v}_{1}$ as long as
$v_{1}(\boldsymbol{x},0)>0$.

Traveling wave solution. We consider the time evolution of $v_{1}(x,t)$ when
$x\in\mathbb{R}$, i.e., we let $L\rightarrow\infty$, and the initial state is
one in which $v_{1}(x,0)$ goes to the stable solution
$\bar{v}_{1}=(\lambda-1)/\lambda$ as $x\rightarrow-\infty$ and to the unstable
solution $v_{1}=0$ as $x\rightarrow\infty$. Equation (47) for $k=1$ can be
considered as a special case of the non-local KPP equation, with the diffusion
constant set equal to zero [1]. To get a feeling for the evolution of such an
initial state we first consider the limiting case when the width $J(x-y)$ goes
to 0, i.e.,

(51) $\displaystyle J(x-y)=\delta(x-y).$

Equation (47) then has the traveling wave solution

(52) $\displaystyle
u_{1}(x,t)=\frac{\lambda-1}{\lambda}\Big{[}1-\tanh[\alpha(x-Vt)]\Big{]}/2,\quad\alpha
V=(\lambda-1)/2$

Numerical solutions of (47) with

$\displaystyle
J(x)=\frac{1}{2b}\theta(b-|x|),\quad\theta(x)=\begin{cases}1,&\quad x\geq
0,\\\ 0,&\quad x<0,\end{cases}$

and initial conditions

$v_{1}(x,0)=\begin{cases}\frac{\lambda-1}{\lambda},\quad&\text{for $x<0$},\\\
0,\quad&\text{for $x>0$},\end{cases}$

show that $v_{1}(x,t)$ approaches a form close to (52) with $\alpha
V\sim(1-\lambda)/2$ as $t\rightarrow\infty$; see figure 2. Similar behavior is
found for $J(x)$ a Gaussian.

Figure 2. Numerical simulation of $k=1$ traveling front solutions. In all
plots, $\lambda=1.1$ and $J(x)=1$ for $|x|<1/2$ and 0 otherwise. A. Traveling
front arising from step function initial conditions. At $t=0$,
$v_{1}(x)=\bar{v}_{1}$ for $x<0$ and 0 otherwise. The traveling front is well
approximated by a rescaled and shifted hyperbolic tangent function, shown in
red. B. Fronts developing from initial conditions set to eq. (52) with
$\alpha=0.1$ and $\alpha=0.05$. Front velocity is approx. double for
$\alpha=0.05$ compared to $0.1$. C. Dependence of front velocity on $\alpha$.
$v_{1}(x)$ is initialized with eq. (52) for a range of $\alpha$, and velocity
is computed numerically. Comparison with velocity given in eq. (52) is given
by the red dashed line. Note that for $1/\alpha$ close to 0, the initial
condition is nearly a step function as in panel A.

##  7\. Conclusions

We have shown that when $\lambda>1$, the macroscopic equations (2)-(4) have a
unique, linearly-stable stationary state with nonzero firing rates
$v_{k}(\boldsymbol{x},t)$ given by (10), which we refer to as the sustaining
stationary state. For $k=1$, the basin of attraction of the sustaining
stationary state includes all initial conditions such that
$v_{1}(\boldsymbol{x},0)$ is not identically $0$ for all $x$. The case $k>1$
is qualitatively different: even starting with $v_{k}(0)>0$ we can have
$v_{k}(t)\rightarrow 0$. We have shown this explicitly for $k=2$ and found
that there are linearly-stable extinct states like
$\boldsymbol{v}(x,t)=(1,0,\ldots,0)$.

For $k=2$, we have seen in numerical simulations that even for initial
conditions with $\boldsymbol{v}(x,0)=\bar{\boldsymbol{v}}$ for $x$ in a small
region of $\mathbb{R}$ and $\boldsymbol{v}(x,0)=(1,0,0)$ outside of that
region ($J$ the same as in figure 2), the firing (epidemic) can die out and
approach an extinct state. As the size of the region increases, eventually a
point is reached where the firing becomes self-sustaining and spreads
throughout the system. We conjecture that there exists some $M>0$, depending
on $J(x)$, such that the initial condition with the region
$|\boldsymbol{x}|<M$ set to $\bar{\boldsymbol{v}}$ will necessarily converge
to $\bar{\boldsymbol{v}}$ pointwise on the whole domain.

In the one-dimensional case $x\in\mathbb{R}$, we showed the existence of
traveling wavefront solutions $\boldsymbol{v}(x,t)=\boldsymbol{v}(x-Vt)$ in
the $k=1$ case, with an analytic solution in the $J=\delta$ case and
numerically for other forms of $J$. Wavefront solutions have been studied
extensively in the case of neural models like the mean field Wilson-Cowan
equations [13]. It remains to be rigorously shown that traveling fronts exist
in the case $k>1$. Based on numerical solutions of the equations, we
conjecture that stable traveling wave solutions exist for $k>1$, and that
there are qualitative differences with the case $k=1$. The existence of
linearly stable extinct states like $\boldsymbol{v}=(1,0\ldots,0)$ when $k>1$
changes the properties of traveling fronts with
$\boldsymbol{v}(-\infty)=\bar{\boldsymbol{v}}$ and
$\boldsymbol{v}(\infty)=(1,0,\ldots,0)$. In particular, we conjecture that
there exists unique traveling waves, with wave velocity proportional to the
width of $J$.

As the microscopic system described in the Introduction has been inspired by
integrate and fire models from neuroscience, we are interested in extensions
of the model to include more realistic components of biological neurons and
neural networks. In particular, including sites in the model which have an
inhibitory effect on nearby sites will lead to a richer dynamical landscape.
Other models with mixed excitatory and inhibitory components have been shown
to have oscillatory activity, and such models have been used to study brain
rhythms (e.g., [2], [10], [14]). Additionally, the effects of inhibition on
the propagation of traveling fronts in the brain has been examined in [13]. An
interesting question is how brain activity is contained and localized to a
particular region when externally driven. We would like to explore this
question in our model with the addition of both inhibition and an external
drive component.

Acknowledgements. We would like to thank Lai-Sang Young, Anna De Masi, and
Errico Presutti for many helpful discussions.

## References

  * [1] Henri Berestycki, Grégoire Nadin, Benoit Perthame, and Lenya Ryzhik. The non-local Fisher-KPP equation: traveling waves and steady states. Nonlinearity, 22(12):2813-2844, 2009.
  * [2] N. Brunel. Dynamics of Sparsely Connected Networks of Excitatory and Inhibitory Spiking Neurons. J Comput Neurosci 8:183–208, 2000. https://doi.org/10.1023/A:1008925309027
  * [3] L. Chariker and L.-S. Young. Emergent spike patterns in neuronal populations J. Comp. Neurosci. 38:203–220, 2015.
  * [4] Logan Chariker, Anna De Masi, Joel Lebowitz, and Errico Presutti. In preparation, 2021.
  * [5] A. De Masi, A. Galves, E. Löcherbach, and E. Presutti. Hydrodynamic Limit for Interacting Neurons. J Stat Phys 158:866–902, 2015. https://doi.org/10.1007/s10955-014-1145-1
  * [6] A. Duarte, G. Ost, and A.A. Rodríguez. Hydrodynamic Limit for Spatially Structured Interacting Neurons. J Stat Phys, 161:1163–1202, 2015. https://doi.org/10.1007/s10955-015-1366-y
  * [7] Wulfram Gerstner and Werner Kistler. Spiking Neuron Models: Single Neurons, Populations, Plasticity. Cambridge University Press, Cambridge, UK, 2002.
  * [8] G. Giacomin, J. Lebowitz, and E. Presutti. Deterministic and Stochastic Hydrodynamic Equations Arising From Simple Microscopic Model Systems, in Stochastic Partial Differential Equations: Six Perspectives, Chapter 3. Edited by R. Carmona and B. Rozovskii, Math. Surveys Monogr., Vol. 64, p. 107, Amer. Math. Soc., Providence, RI, 1999.
  * [9] Igor Goychuk and Andriy Goychuk 2015 New J. Phys. 17 045029
  * [10] S. Keeley, À. Byrne, A. Fenton, J. Rinzel. Firing rate models for gamma oscillations. J Neurophysiol., 121(6):2181-2190, 2019.
  * [11] J. Milton et al. Spiral Waves in Integrate-and-Fire Neural Networks. NIPS, 1992.
  * [12] P.R. Protachevicz, M. Hansen, K.C. Iarosz, I.L. Caldas, A.M. Batista, J. Kurths. Emergence of Neuronal Synchronisation in Coupled Areas. Frontiers in Computational Neuroscience, 15:35, 2021.
  * [13] L.R. González-Ramírez, M.A. Kramer. The effect of inhibition on the existence of traveling wave solutions for a neural field model of human seizure termination. J Comput Neurosci 44, 393–409, 2018. https://doi.org/10.1007/s10827-018-0685-9
  * [14] E. Wallace, M. Benayoun, W. van Drongelen, J.D. Cowan. Emergent Oscillations in Networks of Stochastic Spiking Neurons. PLoS ONE 6(5):e14804, 2011. https://doi.org/10.1371/journal.pone.0014804
  * [15] C. Zankoc, T. Biancalani, D. Fanelli, and R. Livi. Diffusion approximation of the stochastic Wilson-Cowan model. Chaos, Solitons and Fractals, 103:504-512, 2017.