# Observational Inference on the Delay Time Distribution of Short Gamma-ray
Bursts

Michael Zevin<EMAIL_ADDRESS>Kavli Institute for Cosmological
Physics, The University of Chicago, 5640 South Ellis Avenue, Chicago, IL
60637, USA  Enrico Fermi Institute, The University of Chicago, 933 East 56th
Street, Chicago, IL 60637, USA  Anya E. Nugent Center for Interdisciplinary
Exploration and Research in Astrophysics (CIERA) and Department of Physics and
Astronomy, Northwestern University, Evanston, IL 60208, USA  Susmita Adhikari
Department of Physics, Department of Astronomy & Astrophysics, The University
of Chicago, 5640 South Ellis Avenue, Chicago, IL 60637, USA  Department of
Physics, Indian Institute of Science Education and Research, Homi Bhaba Road,
Pune, Maharashtra, 411045, India  Wen-fai Fong Center for Interdisciplinary
Exploration and Research in Astrophysics (CIERA) and Department of Physics and
Astronomy, Northwestern University, Evanston, IL 60208, USA  Daniel E. Holz
Kavli Institute for Cosmological Physics, The University of Chicago, 5640
South Ellis Avenue, Chicago, IL 60637, USA  Enrico Fermi Institute, The
University of Chicago, 933 East 56th Street, Chicago, IL 60637, USA
Department of Physics, Department of Astronomy & Astrophysics, The University
of Chicago, 5640 South Ellis Avenue, Chicago, IL 60637, USA  Luke Zoltan
Kelley Center for Interdisciplinary Exploration and Research in Astrophysics
(CIERA) and Department of Physics and Astronomy, Northwestern University,
Evanston, IL 60208, USA

###### Abstract

The delay time distribution of neutron star mergers provides critical insights
into binary evolution processes and the merger rate evolution of compact
object binaries. However, current observational constraints on this delay time
distribution rely on the small sample of Galactic double neutron stars (with
uncertain selection effects), a single multimessenger gravitational-wave
event, and indirect evidence of neutron star mergers based on $r$-process
enrichment. We use a sample of $67$ host galaxies of short gamma-ray bursts to
place novel constraints on the delay time distribution, and leverage this
result to infer the merger rate evolution of compact object binaries
containing neutron stars. We recover a power-law slope of
$\alpha=-1.82^{+0.40}_{-0.41}$ and a minimum delay time of
$\mathit{t_{\mathrm{min}}}=171^{+74}_{-91}~{}\mathrm{Myr}$ (median and 90%
credible interval), with the maximum delay time constrained to
$\mathit{t_{\mathrm{max}}}>7.73~{}\mathrm{Gyr}$ at 99% credibility. We find
these constraints to be broadly consistent with theoretical expectations,
although our recovered power-law slope is substantially steeper than the
conventional value of $\mathit{t_{\mathrm{d}}}^{-1}$, and our minimum delay
time is larger than the typically assumed value of $10~{}\mathrm{Myr}$.
Pairing this cosmological probe of the fate of compact object binary systems
with the Galactic population of double neutron stars will be crucial for
understanding the unique selection effects governing both of these
populations. In addition to probing a significantly larger redshift regime of
neutron star mergers than possible with current gravitational-wave detectors,
complementing our results with future multimessenger gravitational-wave events
will also help determine if short gamma-ray bursts ubiquitously result from
compact object binary mergers.

††software: Astropy (Robitaille et al., 2013; Price-Whelan et al., 2018);
Bilby (Ashton et al., 2019); iPython (Pérez & Granger, 2007); Matplotlib
(Hunter, 2007); NumPy (Oliphant, 2006; Van Der Walt et al., 2011); Pandas
(McKinney, 2010); Prospector (Johnson et al., 2021); SciPy (Virtanen et al.,
2020).††thanks: NASA Hubble Fellow

##  1 Introduction

The transient event GW170817 (Abbott et al., 2017a), which was observed in
both gravitational waves and electromagnetic (EM) waves, firmly established
the connection between binary neutron star (BNS) mergers and short gamma-ray
bursts. The host galaxy of this event, NGC4993, was subsequently identified
through broadband EM emission of the ensuing kilonova (Allam et al., 2017;
Arcavi et al., 2017; Coulter et al., 2017a; Lipunov et al., 2017; Tanvir et
al., 2017; Yang et al., 2017) and eventually, the non-thermal afterglow of the
GRB (see Margutti & Chornock, 2021, and references therein). Under the
paradigm that all short GRBs originate from compact object binary mergers,
such host associations can unveil unprecedented information regarding the
stellar populations in galaxies that host compact object binary mergers,
thereby allowing for novel constraints on the progenitors of these events.

One key aspect of compact object binary mergers that can be illuminated using
such host associations is the delay-time distribution (DTD). The DTD yields
important information regarding the birth properties of compact object
binaries, such as their typical orbital separations at birth and inspiral
times. This, in turn, provides critical constraints on evolutionary processes
of the progenitor binary stellar system and can inform modeling of massive-
star binaries, which form GW sources accessible by current and future
observatories.

GW170817 is the only GW event to date where the host galaxy has been
confidently established, and a number of fortuitous aspects of the system that
produced this signal were paramount in the discovery of its EM counterpart:
for example, it was well within the horizon of the GW detectors (Abbott et
al., 2018), the GRB jet was aligned $\approx 20~{}\mathrm{deg}$ from our line
of sight and eventually came into view (Margutti & Chornock, 2021), and it was
in close proximity to a massive galaxy that was prioritized by targeted
searches (Coulter et al., 2017b), although wide-field follow-up searches would
have found it regardless (Soares-Santos et al., 2017). While the BNS detection
rate in the forthcoming observing run of the LIGO–Virgo–Kagra interferometer
network will be $\mathcal{O}(1)$/month (Abbott et al., 2018), the rate of BNS
mergers with an EM counterpart detection precise enough for confident host
association is more uncertain given the unique observational challenges (e.g.,
Coughlin et al., 2018; Dichiara et al., 2020; Colombo et al., 2022; Perna et
al., 2022). In the next observing run, the BNS range for the LIGO detectors is
expected to be $160\mbox{--}190~{}\mathrm{Mpc}$ (Abbott et al., 2018), and GW
mergers are more likely to be detected close to the detector horizon where the
sensitive volume is the largest. Three-dimensional GW sky localization volumes
will be larger at these distance, and even with improved search capabilities
from state-of-the-art wide-field telescopes, it is unclear how rapidly
multimessenger events with definitive host associations will be accumulated in
the coming years.

On the other hand, the number of short GRBs with confident host galaxy
associations has increased substantially over the past two decades. Over the
next decade, the cumulative number of these host associations may still exceed
the number of host associations found via multimessenger observations. GRB
monitoring missions such as Swift and Fermi have observed hundreds of short
GRBs, with subsequent searches in the X-ray, optical, and radio bands for GRB
afterglows leading to the identification of probable host galaxies for dozens
of these events over a wide range of redshifts (Berger, 2014; Fong et al.,
2022; Nugent et al., 2022; O’Connor et al., 2022). Photometric and
spectroscopic follow-up observations of the afterglow and the host galaxy can
then unveil crucial aspects of the GRB-host connection that encode information
about the progenitors of BNS systems themselves, such as the stellar mass and
star formation rate (SFR) of the galaxy, its star formation history (SFH), and
the galactocentric offset of the GRB (Fong et al., 2013; Fong & Berger, 2013).

Here, we consider direct implications of GRB–galaxy connections on compact
object binary formation and the evolution of their progenitor stars,
leveraging an updated sample of GRB host galaxies and their properties
presented in Fong et al. (2022) and Nugent et al. (2022) to place novel
constraints on the DTD of short GRBs. This expansive catalog of short GRB
hosts is likely to outnumber the hosts of multimessenger GW events for years
to come, and depending on the highly uncertain rate of discovery of EM
counterparts to GW events, may not be met with multimessenger events until the
onset of third-generation GW detectors. Besides probing a much larger
cosmological volume than a multimessenger sample, future comparisons of DTD
constraints between short GRBs and GWs with EM counterparts could help to
determine whether the BNS merger paradigm for short GRBs is universal.
Furthermore, DTD constraints from such an extra-galactic population are
complementary to those from compact object binaries observed in the Milky Way
and with GWs, allowing for comparisons of these distinct probes to help
unravel selection effects that impinge upon the detection of each population
individually.

The format of this Letter is as follows: In Section 2, we summarize the
catalog of GRB host galaxies used in our analysis and assumptions regarding
their SFHs. Our main results are in Section 3, where we describe our inference
methodology and present our constraints on the DTD of short GRBs. We show
implications of our constraints on the predicted compact object binary merger
rate in Section 4. In Section 5 we discuss caveats of our analysis, and
highlight other aspects of binary stellar evolution and compact object binary
formation that will be addressed in future work. We summarize our results in
Section 6. Throughout this work we employ a standard WMAP cosmology of $H_{0}$
= 69.6 km s-1 Mpc-1, $\Omega_{M}$ = 0.286, $\Omega_{\Lambda}$ = 0.714 (Hinshaw
et al., 2003; Bennett et al., 2014), consistent with Nugent et al. (2022).

##  2 Catalog of Short GRB Hosts

Recently, Fong et al. (2022) presented an updated catalog of $90$ short GRBs
with broadband photometry and well-studied explosion environments. Combined
with the literature, the catalog comprises $335$ photometric data points, $40$
spectra, and $83$ offset measurements, and is comprehensive for Swift short
GRBs and neutron star GW merger events discovered in 2005–2021.111The sample
in Fong et al. (2022) also includes long-duration GRB060614 and GRB211211A,
which are thought to be the result of neutron star mergers due to the lack of
a supernova to deep optical limits and the observation of a kilonova,
respectively (Gehrels et al., 2006; Rastinejad et al., 2022). This catalog, as
well as other recent short GRB catalogs (e.g., O’Connor et al., 2022), have
more than doubled the number of short GRBs with confident host associations,
although a fraction are found to be “host-less” with no coincident host to
deep optical limits, depending on the threshold defining host galaxy
association (Berger, 2010).

Of the $90$ short GRBs with broadband photometry presented in Fong et al.
(2022), $67$ were determined to have sufficient observational data (detected
in $\geq 3$ photometric bands) to model the galactic spectral energy
distribution and infer aspects of the hosts; the detailed properties of this
population of short GRB host galaxies are analyzed in Nugent et al. (2022).
Approximately $75\%$ of hosts have spectroscopic redshifts, with $60\%$ of the
total having detailed spectra. One of the $67$ GRBs was determined to not have
a confident enough host association (see Section 2.2) and is not used in this
analysis. Additionally, we include in our sample NGC4993, the host galaxy of
multimessenger event GW170817/GRB170817.222We do not include the long-duration
GRB060614 and GRB211211A in our main analysis despite that they may be the
result of compact object binary mergers, though we perform additional fits to
the DTD with the inclusion of the host galaxies of these two events and
comment on the impact in Section 5. This totals in $67$ host galaxies used in
our analysis, the properties of which are modeled consistently as described
below.

###  2.1 Host Galaxy Properties

Properties of the host galaxies, such as host redshift, stellar mass, stellar
population ages, metallicity, and SFR, are inferred using the Prospector
package (Leja et al., 2017; Johnson et al., 2021), as described in Nugent et
al. (2022). Prospector uses the dynesty nested sampler (Speagle, 2020) to
constrain aspects of galaxies based on available photometric and spectroscopic
data. Redshifts are fixed in this inference for galaxies that have
spectroscopically-determined redshifts, as the redshift measurement
uncertainties for these hosts are negligible.

The SFH of each host, $\psi(t)$, is modeled assuming a delayed-$\tau$
functional form,333 A number of other parametric functional forms (Carnall et
al., 2019) or non-parametric approaches (Leja et al., 2019) can be used when
reconstructing the SFH, sometimes yielding discrepant results. However,
parametric fits are used in this work to better establish uniformity in our
dataset, which has an inconsistent amount and quality of data across host
observations. Parametric SFHs are also more commonly used in host galaxy
literature, making direct comparisons more seamless. We comment on this
systematic uncertainty in Section 5.

$\psi(t)=M_{\mathrm{F}}\times\frac{te^{-t/\tau}}{\int_{0}^{t_{\mathrm{SF}}}t^{\prime}e^{-t^{\prime}/\tau}dt^{\prime}},$
(1)

where $M_{\mathrm{F}}$ is the total stellar mass formed in the galaxy,
$t_{\mathrm{SF}}$ is the time at which star formation commenced relative to
the observation time, and $\tau$ is the e-folding factor. Posterior samples
for the parameters $M_{\mathrm{F}}$, $t_{\mathrm{SF}}$, and $\tau$, as well as
the redshift $z$ for hosts that only have photometric redshifts, fully
construct the SFH for each Prospector posterior sample. We convert Eq. 1 into
a redshift-dependent SFH, $\psi(z)$, in our DTD parameter hyperlikelihood
described in Section 3.1.

###  2.2 Criteria for Host Associations

The probability of chance coincidence ($\mathit{P_{\mathrm{cc}}}$) metric is
used for determining the confidence of a particular host association (Bloom et
al., 2002). This metric decreases (i.e., results in a higher likelihood that
the GRB is correctly associated with a particular host) with decreasing
optical magnitude of the potential host and decreasing angular offset of the
GRB. We use the $\mathit{P_{\mathrm{cc}}}$ criteria of Fong et al. (2022) to
broadly categorize GRBs in the catalog based on host association confidence: a
“Gold Sample” with $\mathit{P_{\mathrm{cc}}}\leq 0.02$, a “Silver Sample” with
$0.02<\mathit{P_{\mathrm{cc}}}\leq 0.10$, and a “Bronze Samples” with
$0.10<\mathit{P_{\mathrm{cc}}}\leq 0.20$.

In this work, our full sample considers the $67$ hosts observed in $\geq 3$
photometric bands with $\mathit{P_{\mathrm{cc}}}\leq 0.20$. We also present
results using three subsets of this sample based on their redshift uncertainty
and/or $\mathit{P_{\mathrm{cc}}}$: a subset of hosts with spectroscopic
redshifts (and therefore negligible redshift uncertainty, $47$ hosts), the
subset of hosts with $P_{cc}\leq 0.02$ (Gold Sample, $38$ hosts), and the
subset of hosts in the Gold Sample that also have spectroscopic redshifts
($29$ hosts). Note that in our main analysis we include GRB170817, which falls
under the Gold Sample and has a spectroscopic redshift, and do not include
long-duration GRB060614 and GRB211211A, hence the slight differences in the
number of hosts in our subsamples compared to Nugent et al. (2022) (see their
Figure 1). We comment on completeness and potential biases from excluding
certain GRB hosts from our sample in Section 5.

##  3 Delay-Time Distributions

The ‘delay-time’ of compact object binaries is defined as the time between the
formation of the progenitor stars and the merger of the two compact objects,
such that
$\mathit{t_{\mathrm{d}}}=\mathit{t_{\star}}+\mathit{t_{\mathrm{insp}}}$ where
$\mathit{t_{\star}}$ is the time between stellar birth and compact object
binary formation and $\mathit{t_{\mathrm{insp}}}$ is the compact object binary
inspiral time. The DTD for isolated compact binaries is typically
parameterized as a power-law distribution444Another commonly-used functional
form for the DTD of short GRBs is a log-normal distribution (Nakar et al.,
2006; Berger et al., 2007; Wanderman & Piran, 2015), though this functional
form is in tension with the growing number of high-redshift short GRBs and is
not consistent with the expected DTD from compact object binary mergers. with
a minimum delay time $\mathit{t_{\mathrm{min}}}$, which encodes the minimum
stellar evolutionary timescales required to form a compact object binary
($\mathcal{O}(10~{}\mathrm{Myr})$ for BNS progenitors) and the minimum orbital
separation of a compact object binary at formation, and can include a
(somewhat arbitrary) maximum delay time $\mathit{t_{\mathrm{max}}}$ that can
be much larger than the Hubble time,

$p(\mathit{t_{\mathrm{d}}}|\alpha,\mathit{t_{\mathrm{min}}},\mathit{t_{\mathrm{max}}})=\begin{cases}\mathcal{N}\mathit{t_{\mathrm{d}}}^{\alpha},&\mathit{t_{\mathrm{min}}}\leq\mathit{t_{\mathrm{d}}}\leq\mathit{t_{\mathrm{max}}}\\\
0,&\text{otherwise}\\\ \end{cases}$ (2)

where

$\mathcal{N}=\begin{cases}\left[\mathrm{log}_{10}\left(\frac{\mathit{t_{\mathrm{max}}}}{\mathit{t_{\mathrm{min}}}}\right)\right]^{-1},&\alpha=-1\\\
(1+\alpha)(\mathit{t_{\mathrm{max}}}^{1+\alpha}-\mathit{t_{\mathrm{min}}}^{1+\alpha})^{-1},&\text{otherwise}\\\
\end{cases}$ (3)

is the normalization constant for the probability distribution function.

The slope of the DTD for compact object binary mergers is typically assumed to
be $\alpha\approx-1$ (e.g., Piran, 1992). This assumption stems from the fact
that in general $\mathit{t_{\mathrm{insp}}}\gg\mathit{t_{\star}}$, and from
coupling the equations that govern GW inspiral with an assumption that the
orbital separation distribution of compact object binaries follows the
distribution of massive O/B stars, which is approximately flat in log with
$dN/da\propto a^{-1}$ where $a$ is the semimajor axis (for example, the best-
fit orbital separation distribution from a survey of massive-star binary
initial properties in Sana et al. 2012 finds $dN/da\propto a^{-0.83}$). Since
the inspiral time scales as $t_{\mathrm{insp}}\propto a^{4}$, it follows that
$da/dt_{\mathrm{insp}}\propto t_{\mathrm{insp}}^{-3/4}$ and, assuming
$dN/da\propto a^{-1}$, $dN/da\propto t_{\mathrm{insp}}^{-1/4}$. The
distribution of inspiral times is thus $dN/dt_{\mathrm{insp}}\propto
t_{\mathrm{insp}}^{-1}$. However, hardening phases during the coevolution of
massive-star binaries, such as common envelopes or stable mass transfer phases
from a more massive donor star to a less massive compact object accretor, can
lead to steeper orbital separation distributions and therefore steeper slopes
in the DTD (e.g., Belczynski et al., 2018). The DTD inferred from short GRBs
is thus useful in constraining a variety of uncertain aspects of massive-star
binary evolution.

Multiple observational probes have been explored in constraining aspects of
the short GRB/BNS DTD. The redshift distribution of short GRBs has been able
to place broad constraints on the DTD (Nakar et al., 2006; Berger et al.,
2007; Jeong & Lee, 2010; Hao & Yuan, 2013; Wanderman & Piran, 2015; Anand et
al., 2018), with high-redshift GRBs in particular allowing for stronger
constraints on the minimum delay time and power-law slope of the DTD (Paterson
et al., 2020; O’Connor et al., 2022; Nugent et al., 2022). The binary–host
connection interpreted via host galaxy demographics and galaxy scaling
relations has also shown promise in constraining BNS kick velocities, delay
times, and the properties of galaxies that host BNS mergers and short GRBs
(Zheng & Ramirez-Ruiz, 2007; Kelley et al., 2010; Fong et al., 2013; Behroozi
et al., 2014; Adhikari et al., 2020; Santoliquido et al., 2022). The Milky Way
offers a limited sample of BNS systems with well-characterized orbital
properties and inspiral times, many of which are much longer than the Hubble
time, that are useful for examining the BNS DTD (Beniamini & Piran, 2016;
Tauris et al., 2017; Vigna-Gómez et al., 2018; Andrews & Zezas, 2019;
Beniamini & Piran, 2019). However, uncertain selection effects inherent to
this sample may affect inference of DTD parameters (e.g., Tauris et al.,
2017). All these observational probes are complemented by population modeling
of compact object binaries that explore variations in binary evolution
physical assumptions and therefore the predicted DTD of compact object
binaries (Belczynski et al., 2018; Chruslinska et al., 2018; Broekgaarden et
al., 2021; Santoliquido et al., 2022); these predicted compact object
populations can be compared and constrained with the population properties and
merger rates observed by GW detectors (Abbott et al., 2021a). Though GW
observations of compact object binary mergers also hold promise in
constraining the DTD of compact object mergers, both through the use of galaxy
scaling relations (Safarzadeh & Berger, 2019; Safarzadeh et al., 2019b) and
multimessenger events whose host galaxies have been identified and
characterized (Safarzadeh et al., 2019a), it may require the accumulation of
hundreds of host galaxy identifications to achieve strong constraints on the
DTD.

###  3.1 Inferring the DTD Parameters

We follow Safarzadeh et al. (2019a) in constructing our likelihood function
for the DTD, which is the key component of our population inference. From the
Prospector modeling of each host galaxy, we have a set of posterior samples
that define the SFH and, for hosts without spectroscopically-measured
redshifts, the galaxy redshift; a fixed redshift is used for the galaxies that
have spectroscopic redshifts. We draw $N_{\mathrm{samp}}=100$ Prospector
samples from each host, where we weight our draws with the inverse of the
prior to get draws from the likelihood rather than the posterior.555All
relevant Prospector parameters used in this analysis were sampled with uniform
priors, making this prior reweighting irrelevant. For a given host galaxy $i$
and likelihood sample $j$, the expected merger rate at the measured redshift
$z_{i}^{j}$ is

$\dot{n}_{i}^{j}=\int_{z^{\prime}=\infty}^{z^{\prime}=z_{i}^{j}}p(t^{\prime}_{\mathrm{lb}}-\mathit{t_{\mathrm{d}}}|\alpha,\mathit{t_{\mathrm{min}}},\mathit{t_{\mathrm{max}}})\lambda\psi_{i}^{j}(z^{\prime})\frac{dt}{dz}(z^{\prime})dz^{\prime},$
(4)

where $t^{\prime}_{\mathrm{lb}}$ is the lookback time at redshift
$z^{\prime}$, $\lambda$ is the BNS formation efficiency which we assume is
$10^{-5}\mathit{M_{\odot}}^{-1}$ and does not evolve with redshift (see e.g.,
Broekgaarden et al., 2021), $\psi(z)$ is the SFR as a function of redshift
defined in Section 2.1, $dt/dz=[H_{0}E(z)(1+z)]^{-1}$, and
$E(z)=\sqrt{\Omega_{M}(1+z)^{3}+\Omega_{\Lambda}(z)}$.

The probability of observing a single short GRB in a host galaxy over a set
interval of time $\Delta t$ follows a Poisson distribution based on the merger
rate of Equation 4. Approximating the marginalization over likelihood samples
as a discrete sum, the hyperlikelihood for observing a single short GRB in a
particular host galaxy given the set of DTD parameters becomes

$\mathcal{L}(\mathrm{obs}_{i}|\alpha,\mathit{t_{\mathrm{min}}},\mathit{t_{\mathrm{max}}})\approx\frac{1}{N_{\mathrm{samp}}}\sum_{j=1}^{N_{\mathrm{samp}}}(\dot{n}_{i}^{j}\Delta
t)e^{-\dot{n}_{i}^{j}\Delta t}$ (5)

where we assume that the fiducial observation time is $\Delta
t=10~{}\mathrm{yr}$, though this choice does not impact results. We note that
this formulation of the hyperlikelihood assumes that the DTD parameters are
not prone to selection effects, that is, the probability of a detection does
not depend on $\alpha$, $\mathit{t_{\mathrm{min}}}$, and
$\mathit{t_{\mathrm{max}}}$. Selection effects may impact the detection of
more highly offset (i.e., “host-less”) systems with less luminous counterparts
and/or weaker host associations, as well as systems at high redshifts with
dimmer hosts; we discuss potential implications of selection effects further
in Section 5.

Assuming the observed short GRB observations are independent, the
hyperposterior for the DTD parameters is

$P(\alpha,\mathit{t_{\mathrm{min}}},\mathit{t_{\mathrm{max}}}|\mathrm{obs})\propto\prod_{i=1}^{N_{\mathrm{GRB}}}\mathcal{L}(\mathrm{obs}_{i}|\alpha,\mathit{t_{\mathrm{min}}},\mathit{t_{\mathrm{max}}})\\\
\times\pi(\alpha,\mathit{t_{\mathrm{min}}},\mathit{t_{\mathrm{max}}}),$ (6)

where $N_{\mathrm{GRB}}$ is the number of GRB hosts used in our sample and
$\pi(\alpha,\mathit{t_{\mathrm{min}}},\mathit{t_{\mathrm{max}}})$ is the prior
distribution on the DTD parameters. We assume that the priors on $\alpha$,
$\mathit{t_{\mathrm{min}}}$, and $\mathit{t_{\mathrm{max}}}$ are independent
and uniform on the ranges of $[0,3]$, $[5~{}\mathrm{Myr},500~{}\mathrm{Myr}]$,
and $[3~{}\mathrm{Gyr},100~{}\mathrm{Gyr}]$, respectively.

Sampling the likelihood of Equation 5 is computationally prohibitive, as it
involves calculating an integral at each step in the parameter space of the
DTD distribution. To avoid this, we precompute a regular grid of likelihoods
for each host galaxy and create an interpolant of these grids when evaluating
the likelihood in Equation 5. This grid contains $20$ points for each DTD
parameter uniformly spaced across their prior range, resulting in $8,000$
likelihood interpolants for each host galaxy. We find that this approach
provides consistent results compared to directly evaluating the likelihood, as
the likelihood surface varies smoothly across the DTD parameter space. For
generating posterior distributions of the DTD parameters, we use the dynesty
nested sampler (Speagle, 2020) as implemented in Bilby (Ashton et al., 2019).

Figure 1: Constraints on DTD parameters when including all host associations
with $\mathit{P_{\mathrm{cc}}}\leq 0.20$ in the inference (filled blue), as
well as three subsets that only include samples that have
$\mathit{P_{\mathrm{cc}}}\leq 0.02$ (yellow; Gold Association), samples that
have spectroscopic redshifts (green), and samples that have spectroscopic
redshifts and $\mathit{P_{\mathrm{cc}}}\leq 0.02$ (pink). Contours show the
50% and 90% credible regions, with dotted lines marking the median for each
marginalized distribution. Gray dashed lines in the panels plotting the
marginal distributions show the prior distribution for each parameter, which
is uniform for all parameters of the DTD. Both kernel density estimates and
binned histograms are shown for the marginal distributions.

We note that this approach only considers host galaxies in which a short GRB
has been observed, and not galaxies in which no short GRB occurred during the
fiducial observing period. Locally, this has been shown to lead to only a mild
bias in the inferred DTD parameters (Safarzadeh et al., 2019a).

###  3.2 DTD Constraints

Our main results are in Figure 1, which shows the constraints on the DTD
parameters from our hierarchical inference of Section 3.1. We show posterior
distributions for the DTD parameters using the full population, as well as the
three subsets of the population described in Section 2.2.

Considering our full sample of $67$ GRBs in the inference, we find that the
posterior distribution of DTD parameters significantly deviates from the
prior, with a power-law slope of ${\alpha=-1.82^{+0.40}_{-0.41}}$ and a
minimum delay time of
${\mathit{t_{\mathrm{min}}}=171^{+74}_{-91}~{}\mathrm{Myr}}$, where we quote
the median and 90% symmetric credible interval. These parameters, as expected,
exhibit a strong correlation since the bulk of the short GRB sample is
consistent with relatively short delay times of
$\mathcal{O}(100~{}\mathrm{Myr})$; higher values of $\alpha$ (i.e., a
shallower power-law slope) correspondingly decrease the inferred values for
$\mathit{t_{\mathrm{min}}}$ so that the DTD probability distribution function
still has enough support at the low-end of the distribution to adequately
explain population as a whole.

Constraints on the DTD parameters are fairly consistent when using the full
sample of GRB hosts as opposed to the subsamples described in Section 2.2.
When only using samples with $\mathit{P_{\mathrm{cc}}}\leq 0.02$ (i.e., Gold
Associations) in the inference, the constraints on DTD parameters show little
variation compared to the full sample where hosts only satisfy
$\mathit{P_{\mathrm{cc}}}\leq 0.20$. This suggests that the DTD constraints
are not strongly dependent on the specific criteria used for host
associations, namely the optical magnitude of the host and the offset of the
GRB with respect to the host. On the other hand, when considering only
associations for which the host galaxy has a spectroscopically-measured
redshift, we find mild changes in our inferred DTD parameters. In particular,
the minimum delay time pushes to slightly larger values and becomes more
tightly constrained with a median value that is $33~{}\mathrm{Myr}$ larger
than when all hosts are included, though still within the range of
uncertainty. Furthermore, the subsample with spectroscopic redshifts strongly
excludes extremely short delay times, with
$\mathit{t_{\mathrm{min}}}>80~{}\mathrm{Myr}$ at 99% credibility. The power-
law index shifts to slightly shallower slopes as opposed to the full sample,
with $\alpha=-1.67^{+0.39}_{-0.39}$. As the spectroscopic sample has precise
redshift measurements, the modeled SFHs are much less uncertain compared to
the photometric redshift sample, which is the main driver of these tighter
constraints. We also note that most of the host galaxies in the spectroscopic
redshift sample have redshifts $z<1$ (Fong et al., 2022), which may lead to
mild differences in our DTD constraints relative to the full sample as
discussed in Section 5.

Neither $\alpha$ or $\mathit{t_{\mathrm{min}}}$ are strongly correlated with
$\mathit{t_{\mathrm{max}}}$, though for $\mathit{t_{\mathrm{max}}}\lesssim
30~{}\mathrm{Gyr}$, decreasing the value of $\mathit{t_{\mathrm{max}}}$ leads
to slightly more support for smaller values of $\mathit{t_{\mathrm{min}}}$ and
larger values of $\alpha$ (i.e., shallower power-law slopes). The posterior
distribution for $\mathit{t_{\mathrm{max}}}$ peaks at $\sim
13~{}\mathrm{Gyr}$, and the upper limit for $\mathit{t_{\mathrm{max}}}$ rails
against the upper bound of our prior. Thus, we cannot make meaningful
statements for the upper limit of $\mathit{t_{\mathrm{max}}}$, which
translates to a lack of constraints on the maximum orbital separation of
compact object binary systems containing a neutron star at formation.
Posterior support for values below $\mathit{t_{\mathrm{max}}}\sim
10~{}\mathrm{Gyr}$ drops precipitously; using all host galaxies, we constrain
$\mathit{t_{\mathrm{max}}}>7.73~{}\mathrm{Gyr}$ at 99% credibility. This
constraint is driven by the systems in our population associated with
quiescent host galaxies and old stellar populations.

Figure 2: DTDs constructed from random draws of our posterior distribution.
Colored lines match the set of samples used in the inference as in Figure 1.
Black solid, dashed, and dotted lines show fiducial DTDs with power-law slopes
$\alpha$ of $-1$, $-1.5$, and $-2$, respectively, for comparison. We show this
set of fiducial distributions for
$\mathit{t_{\mathrm{min}}}=10~{}\mathrm{Myr}$ and
$\mathit{t_{\mathrm{min}}}=100~{}\mathrm{Myr}$, with
$\mathit{t_{\mathrm{max}}}$ fixed at $100~{}\mathrm{Gyr}$.

In Figure 2, we show DTDs constructed from our posterior distributions, as
well as DTDs using fiducial values for comparison. As described above, we find
a preference for minimum delay times of $\mathit{t_{\mathrm{min}}}\sim
150~{}\mathrm{Myr}$ and power-law slopes of $\alpha\sim-1.5$. These
constraints can help inform binary evolution modeling of BNS progenitors, as
slopes with $\alpha<-1$ may be indicative of hardening phases in the evolution
of the progenitor such as mass transfer or common envelope phases (Belczynski
et al., 2018), and larger minimum delay times hint at proposed fast-merging
BNS channels (e.g., the Case BB mass transfer scenario where the progenitor of
the second-born neutron star proceeds through a second mass transfer episode
as an evolved naked helium star, see Dewi et al. 2002; Ivanova et al. 2003)
may operate inefficiently.

Figure 3: Expected merger rate evolution of BNSs based on our DTD constraints.
The black line (right axis) shows the SFR density evolution from Madau &
Fragos (2017). Colored lines (left axis) match the set of samples used in the
inference as in Figure 1 and show the predicted BNS comoving merger rate
density for random draws from our DTD parameter inference, assuming their
progenitor stars are born according to the SFR density of Madau & Fragos
(2017). The shaded gray region shows the local BNS merger rate density
constraints from Abbott et al. (2021b).

##  4 Implications for gravitational-wave observations

Given our constraints on the parameters of the short GRB DTD, we can construct
the expected merger rate evolution of BNS mergers under the assumption that
all short GRBs are the result of this class of compact object binary merger.
We show the expected merger rate density evolution in Figure 3, assuming that
BNSs follow the underlying cosmic SFH of Madau & Fragos (2017). The BNS merger
rate density peaks at a redshift of $z\sim 1.6\mbox{--}1.7$. The overall
normalization of the merger rate evolution is directly related to the assumed
BNS production efficiency, which we take to be fixed at
$10^{-5}~{}\mathit{M_{\odot}}^{-1}$. However, this assumption does not affect
the shape of the merger rate evolution or the location of its peak. Given this
fixed value for the BNS production efficiency, the local merger rate density
is predicted to be $\sim 150~{}\mathrm{Gpc}^{-3}~{}\mathrm{yr}^{-1}$,
consistent with the current constraints from the LIGO Scientific and Virgo
Collaboration (Abbott et al., 2021b).

The inferred BNS merger rate evolution does not vary significantly for the
different subsets of samples described in Section 2.2 since the inferred DTD
is relatively robust. Though the redshift of the peak merger rate density is
well beyond the BNS horizon of current-generation GW detectors operating at
design sensitivity (Abbott et al., 2018), it should be within the range of
planned third-generation GW detectors (Evans et al., 2021). If inconsistencies
are found between the peak of the BNS merger rate density inferred from short
GRBs and the peak of the BNS merger rate observed by future GW detectors, it
may indicate an additional mechanism for producing short GRBs is at play or
that other transients (e.g., long GRBs) are also caused by compact object
binary mergers.

##  5 Discussion

Although the DTD constraints of short GRBs in this work are not dependent on
the specific type of formation mechanism, these constraints can help inform
uncertain aspects of massive-star binary evolution within the BNS paradigm for
short GRBs. Assuming that the merger time is purely driven by GW inspiral
(e.g., in the case of isolated binary evolution), the DTD encodes information
about the orbital properties of compact object binary systems at formation.
Given the short stellar evolutionary timescales of BNS progenitors prior to
the second supernova that forms the compact object binary ($\sim
10$–$50~{}\mathrm{Myr}$, Andrews & Zezas 2019), the delay time is
approximately the inspiral time, except for the potential low-end portion of
the DTD where a short inspiral time becomes comparable to the stellar
evolution timescale.

The DTD can thus be used to place constraints on the minimum separation of
BNSs at formation, as well as evolutionary phases in the progenitors that can
steepen the compact object binary orbital separation distribution and
therefore the DTD (Belczynski et al., 2018; Broekgaarden et al., 2021). For
example, given the $1^{\mathrm{st}}$ percentile of our recovered
$\mathit{t_{\mathrm{min}}}$ distribution when using all host galaxies
($t_{\mathrm{min}}^{1\%}=46~{}\mathrm{Myr}$) and assuming a BNS formation
timescale of $t_{\star}=30~{}\mathrm{Myr}$, the minimum orbital separation at
BNS formation is $\approx 0.9~{}R_{\odot}$. The preference of our results for
longer minimum delay times $\gtrsim 100~{}\mathrm{Myr}$ may indicate that
late-stage hardening phases, such as the Case BB mass transfer scenario (e.g.,
Dewi et al., 2002; Ivanova et al., 2003), may operate inefficiently or be
nonexistent. This would have implications for the ability of BNS systems to
provide the $r$-process enrichment observed in low escape velocity
environments such as ultra-faint dwarf galaxies (Safarzadeh et al., 2018) or
environments with rapid star formation episodes such as globular clusters
(Zevin et al., 2019), since enrichment in these environments requires short
delay times of $\mathcal{O}(10~{}\mathrm{Myr})$. In addition to informing the
distribution of orbital properties at compact object binary formation, delay
times can be used to place joint constraints on the properties of a binary
prior to the second supernova and the strength of the natal kick that formed
the second neutron star (e.g., Andrews et al., 2015; Andrews & Mandel, 2019;
Tauris et al., 2017). Combining this information with the physical offset of
the GRB with respect to its host galaxy can further improve such constraints
(e.g., Abbott et al., 2017b; Zevin et al., 2020). We will explore the
constraints that this catalog places on BNS progenitor pre-supernova orbital
properties and natal kicks based on delay times and offsets in future work.

A significant amount of work has explored constraints on BNS progenitor
properties based on the $\sim 20$ BNS systems observed in the Milky Way
(Andrews et al., 2015; Tauris et al., 2017; Andrews & Zezas, 2019; Andrews &
Mandel, 2019). The Galactic BNS population has uncertain selection effects
that can be complemented with short GRB host galaxy associations. In
particular, the shortest delay time BNS systems in the Milky Way may be missed
from surveys due to Doppler shifting from the orbit making their detection
difficult, as well as their short lifetimes making their existence at any
given point in time relatively rare (Tauris et al., 2017). On the other hand,
long delay-time BNS systems in the Milky Way may be hidden from surveys due to
their long orbital periods, making their proper motion on the sky minuscule.
We find a minimum delay time that is consistent with the shortest delay time
systems in the Milky Way. For example, the BNS system J1757-1854 is estimated
to have one of the shortest delay times of the population observed in the
Milky Way; it will merge in $\sim 70~{}\mathrm{Myr}$ and has a characteristic
age of $\sim 130~{}\mathrm{Myr}$ (Andrews & Mandel, 2019) based on the spin-
down of the pulsar in the system, giving it a delay time of $\sim
210\mbox{--}250~{}\mathrm{Myr}$ with the inclusion of a $\sim
10\mbox{--}50~{}\mathrm{Myr}$ timescale between stellar birth and BNS
formation. This may be an indication that selection effects do not strongly
impinge the detection of short delay time systems in the Milky Way, or that
certain selection effects against observing short delay-time systems also
affect the short GRB population. The complementarity of these local and
cosmological samples may also help us understand selection effects impacting
the detection of counterparts to short GRBs, such as the possibility that
highly offset short GRBs have fainter afterglows due to the tenuous nature of
gas in the outskirts of galaxies (e.g., Berger, 2010; Tunnicliffe et al.,
2014).

The largest uncertainty in our analysis likely lies in the reconstruction of
the SFH. Though the parametric delayed-$\tau$ SFH has been shown to be an
adequate representation of quiescent galaxies, other parametric or non-
parametric approaches may be more physically realistic for SFHs of host
galaxies that still have substantial star formation or bursty histories
(Carnall et al., 2019; Leja et al., 2019). However, due to the robustness of
our results with differing subsets of GRB hosts, we expect our main
conclusions to hold despite this systematic uncertainty. In future work, we
will explore how the choice of SFH reconstruction affects our results by
considering alternative parametric and non-parametric modeling techniques.

The host galaxy of the multimessenger event GW170817/GRB170817, NGC4993, is
quiescent with stellar age estimates of $\gtrsim 3\mbox{--}10~{}\mathrm{Gyr}$
(Blanchard et al., 2017; Im et al., 2017; Levan et al., 2017; Pan et al.,
2017; Nugent et al., 2022). Though more quiescent galaxies host short GRBs at
low redshifts ($z<0.5$) than at higher redshifts (Nugent et al., 2022),
properties of NGC4993 such as its mass-weighted age and specific SFR make it a
mild outlier relative to the other host galaxies used in this work (see Nugent
et al. 2022). Despite this, the exclusion of GRB170817 in our analysis has a
very minor impact on our DTD constraints, only decreasing (increasing) our
median recovered value of $\alpha$ ($\mathit{t_{\mathrm{min}}}$) by $5\%$,
well within the bounds of the uncertainty in these parameters. The DTD
constraints in this work are consistent with a long tail of delay times with
$\mathit{t_{\mathrm{max}}}>7.73~{}\mathrm{Gyr}$ at high credibility, and thus
can adequately account for the relatively old stellar population of NGC4993
and the long inferred delay time of GW170817 (Blanchard et al., 2017; Adhikari
et al., 2020).

Extremely short delay times of $\lesssim 10~{}\mathrm{Myr}$ have been argued
for the progenitor of GRB060505, which is spatially associated with an active
star-forming region (Ofek et al., 2007). This GRB had a $T_{90}$ of $\sim
4~{}\mathrm{s}$, longer than the typical delineation between short and long
GRBs of $2~{}\mathrm{s}$. However, deep imaging ruled out the presence of a
supernova (typically associated with long GRBs) to deep optical limits. Since
$10~{}\mathrm{Myr}$ is approximately the shortest possible evolutionary
timescale that a massive-star binary can form a BNS, for a delay time of $\sim
10~{}\mathrm{Myr}$ the inspiral timescale would need to be extremely short,
requiring a BNS birth semimajor axis of $a_{\mathrm{BNS}}\ll
1~{}\mathit{R_{\odot}}$ unless the BNS was born with a high eccentricity. Such
a scenario is inconsistent with our constraints on the minimum delay time
$\mathit{t_{\mathrm{min}}}$; we find
$\mathit{t_{\mathrm{min}}}>46~{}\mathrm{Myr}$ at 99% credibility when using
all host galaxies in our sample and
$\mathit{t_{\mathrm{min}}}>80~{}\mathrm{Myr}$ at 99% credibility when only
considering host galaxies with spectroscopically-measured redshifts. Given our
DTD constraints, this spatial coincidence with a star-forming region may
instead be a lucky coincidence following the post-supernova migration of the
BNS, or this particular event may be a separate class of long GRBs without a
supernova (e.g., Fryer et al., 2006).

However, a number of other long GRB detections are lending credence to the
possibility that some of these events are caused by compact object binary
mergers. The discovery of a potential kilonova associated with the minute-long
GRB211211A (Rastinejad et al., 2022) hints at BNS mergers causing some
fraction of the long GRB population. Similarly, GRB060614 had a long duration
(although is also classified as a possible short GRB with extended emission;
Lien et al. 2014) and no coincident supernova was detected to deep optical
limits (Gehrels et al., 2006). As the properties of the host galaxies of
GRB211211A and GRB060614 were also modeled in Nugent et al. (2022), we
performed additional analyses with the inclusion of these two hosts in our
general short GRB population. These provided nearly identical results. Even
though our DTD results are insensitive to the inclusion of these two long
GRBs, the number of compact object binary mergers consistent with this picture
of extended gamma-ray emission is uncertain and could potentially add more
support for short delay times. However, it is unlikely that the majority of
long GRBs result from a compact object merger paradigm, as many long GRBs in
the local ($z<0.5$) universe have been followed up extensively by
electromagnetic observatories and GRB211211A is the only one observed thus far
to have a plausible kilonova counterpart (Rastinejad et al., 2022).

In addition to some fraction of BNS mergers masquerading as long GRBs, our
sample used to constrain the DTD may suffer from other issues of
incompleteness. Since we rely on modeling of host galaxies when constraining
the DTD, we do not consider short GRBs that do not have a confident host
association. Though the properties of short GRB hosts do not seem to deviate
strongly as a function of host association confidence (see e.g. Figure 4 of
Nugent et al. 2022), neglecting these events may have a potential impact on
both the low and high ends of our inferred DTD. For example, GRBs that are
highly offset from their hosts may have afterglows with much lower
luminosities, making precise localization (and therefore host identification)
difficult (Perna et al., 2022). Such systems may have migrated over longer
timescales to reach the highly offset locations of the burst, and therefore
may have longer delay times than the general population. Furthermore, the
$\mathit{P_{\mathrm{cc}}}$ method for host identification may incorrectly
associate a GRB with a faint underlying host rather than a bright host at a
larger offset, though Fong et al. (2022) predicted this to be an effect only
at the $\lesssim 7\%$ level. On the other hand, if such poorly-associated GRBs
are instead truly associated with faint galaxies that are below detection
limits, we may be excluding additional systems with short delay times since
these faint, low-mass galaxies are typically star-forming. Furthermore, though
Swift can detect GRBs out to $z\sim 3$, there is likely some fraction of short
GRBs that occur beyond this horizon, when the universe was $\lesssim
2~{}\mathrm{Gyr}$ old. Short GRBs that occur at these early stages in the
history of the universe must have short delay times, and this selection effect
may bias the general population in our analysis to longer delay times. This
would lead to a larger inferred $\mathit{t_{\mathrm{min}}}$, and, due to the
correlation between $\mathit{t_{\mathrm{min}}}$ and $\alpha$, more negative
values of $\alpha$. However, this population of high-redshift short GRBs is
likely small; assuming the SFH from Madau & Fragos (2017), $<10\%$ of stars
are born beyond $z=3$, and the fraction of compact object binary mergers
beyond this redshift will be even smaller due to the delay time between
formation and merger.

##  6 Conclusions

We have placed constraints on the delay-time distribution (DTD) of short GRBs
using the largest catalog of short GRB host associations to date with the
inclusion of inference on host galaxy properties. Assuming these transient
events result from the merger of BNS systems, the DTD for short GRBs can be
directly translated to the DTD of BNS mergers. This allows for predictions of
the expected merger rate evolution, and can constrain aspects of massive-star
binary evolution physics and compact object binary formation. Our main results
are as follows:

  1. 1.

Based on a catalog of $67$ short GRB host galaxies, we constrain the DTD of
short GRBs to have a power-law slope steeper than flat-in-log, with a power-
law index of $\alpha=-1.82^{+0.40}_{-0.41}$, a minimum delay time of
$\mathit{t_{\mathrm{min}}}=171^{+74}_{-91}~{}\mathrm{Myr}$, and a maximum
delay time that is $>7.73~{}\mathrm{Gyr}$.

  2. 2.

Using different subsets of the full dataset that make cuts based on host
association confidence or whether the host galaxy has a spectroscopically-
measured redshift, we find our results to be robust. However, when using the
spectroscopic redshift sample (which is dominated by $z<1$ bursts), the DTD
pushes to slightly larger values of the minimum delay time and slightly
shallower power-law slopes.

  3. 3.

Assuming short GRBs are the result of BNS mergers, we construct the expected
merger rate evolution of BNS mergers, which we predict to peak at a redshift
of $z\sim 1.6\mbox{--}1.7$.

Constraints on the DTD is the tip of the iceberg of what can be accomplished
given a large sample of compact object binary host galaxies. This binary-host
sample can enable novel constrains on properties of the compact object binary
at birth, and pairing such inference with the physical offsets of short GRBs
with respect to their hosts can help to unveil supernova mechanisms and the
strength of kicks that neutron stars receive at formation. Pairing this sample
with multimessenger GW events and the already-existing population of BNSs in
the Milky Way will help to unravel uncertain selection effects and determine
the ubiquity of BNS mergers as the cause of short GRBs, all of which we aim to
explore in future work. Even when multimessenger BNS observations surpass the
number of host-identified short GRBs, they will probe a lower redshift regime,
and pairing these samples will be paramount for understanding the binary-host
connection and its evolution over cosmic time.

We thank Jeff Andrews, Gourav Khullar, Joel Leja, and Vicky Kalogera for
useful discussions. Support for this work and for M.Z. was provided by NASA
through the NASA Hubble Fellowship grant HST-HF2-51474.001-A awarded by the
Space Telescope Science Institute, which is operated by the Association of
Universities for Research in Astronomy, Incorporated, under NASA contract
NAS5-26555. A.E.N. acknowledges support from the Henry Luce Foundation through
a Graduate Fellowship in Physics and Astronomy. W.F. gratefully acknowledges
support by the David and Lucile Packard Foundation, the Alfred P. Sloan
Foundation, and the Research Corporation for Science Advancement through
Cottrell Scholar Award #28284. D.E.H is supported by NSF grants PHY-2006645
and PHY-2110507, as well as by the Kavli Institute for Cosmological Physics
through an endowment from the Kavli Foundation and its founder Fred Kavli.
This work was performed in part at the Aspen Center for Physics, which is
supported by NSF grant PHY-1607611.

## References

  * Abbott et al. (2017a) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017a, Physical Review Letters, 119, 161101, doi: 10.1103/PhysRevLett.119.161101
  * Abbott et al. (2017b) —. 2017b, The Astrophysical Journal Letters, 850, L40, doi: 10.3847/2041-8213/aa93fc
  * Abbott et al. (2018) —. 2018, Living Reviews in Relativity, 21, 3, doi: 10.1007/s41114-018-0012-9
  * Abbott et al. (2021a) Abbott, R., Abbott, T. D., Abraham, S., et al. 2021a, arXiv e-prints, doi: 10.1103/PhysRevX.11.021053
  * Abbott et al. (2021b) —. 2021b, arXiv e-prints.  https://arxiv.org/abs/2111.03634
  * Adhikari et al. (2020) Adhikari, S., Fishbach, M., Holz, D. E., Wechsler, R. H., & Fang, Z. 2020, The Astrophysical Journal, 905, 21, doi: 10.3847/1538-4357/abbfb7
  * Allam et al. (2017) Allam, S., Annis, J., Berger, E., et al. 2017, GCN 21530
  * Anand et al. (2018) Anand, N., Shahid, M., & Resmi, L. 2018, Monthly Notices of the Royal Astronomical Society, 481, 4332, doi: 10.1093/mnras/sty2530
  * Andrews et al. (2015) Andrews, J. J., Farr, W. M., Kalogera, V., & Willems, B. 2015, The Astrophysical Journal, 801, 32
  * Andrews & Mandel (2019) Andrews, J. J., & Mandel, I. 2019, The Astrophysical Journal, 880, L8, doi: 10.3847/2041-8213/ab2ed1
  * Andrews & Zezas (2019) Andrews, J. J., & Zezas, A. 2019, Monthly Notices of the Royal Astronomical Society, 486, 3213, doi: 10.1093/mnras/stz1066
  * Arcavi et al. (2017) Arcavi, I., Howell, D. A., McCully, C., et al. 2017, GCN 21538
  * Ashton et al. (2019) Ashton, G., Hübner, M., Lasky, P. D., et al. 2019, The Astrophysical Journal Supplement Series, 241, 27, doi: 10.3847/1538-4365/ab06fc
  * Behroozi et al. (2014) Behroozi, P. S., Ramirez-Ruiz, E., & Fryer, C. L. 2014, The Astrophysical Journal, 792, 19, doi: 10.1088/0004-637X/792/2/123
  * Belczynski et al. (2018) Belczynski, K., Bulik, T., Olejak, A., et al. 2018, arXiv e-prints.  https://arxiv.org/abs/1812.10065
  * Beniamini & Piran (2016) Beniamini, P., & Piran, T. 2016, Monthly Notices of the Royal Astronomical Society, 456, 4089, doi: 10.1093/mnras/stv2903
  * Beniamini & Piran (2019) —. 2019, arXiv e-prints.  https://arxiv.org/abs/1903.11614
  * Bennett et al. (2014) Bennett, C. L., Larson, D., Weiland, J. L., & Hinshaw, G. 2014, Astrophysical Journal, 794, 2, doi: 10.1088/0004-637X/794/2/135
  * Berger (2010) Berger, E. 2010, The Astrophysical Journal, 722, 1946, doi: 10.1088/0004-637X/722/2/1946
  * Berger (2014) —. 2014, Annual Review of Astronomy and Astrophysics, 52, 43, doi: 10.1146/annurev-astro-081913-035926
  * Berger et al. (2007) Berger, E., Fox, D. B., Price, P. A., et al. 2007, The Astrophysical Journal, 664, 1000, doi: 10.1086/518762
  * Blanchard et al. (2017) Blanchard, P. K., Berger, E., Fong, W., et al. 2017, The Astrophysical Journal, 848, L22, doi: 10.3847/2041-8213/aa9055
  * Bloom et al. (2002) Bloom, J. S., Kulkarni, S. R., & Djorgovski, S. G. 2002, The Astronomical Journal, 123, 1111, doi: 10.1086/338893
  * Broekgaarden et al. (2021) Broekgaarden, F. S., Berger, E., Stevenson, S., et al. 2021, arXiv e-prints.  https://arxiv.org/abs/2112.05763
  * Carnall et al. (2019) Carnall, A. C., Leja, J., Johnson, B. D., et al. 2019, The Astrophysical Journal, 873, 44, doi: 10.3847/1538-4357/ab04a2
  * Chruslinska et al. (2018) Chruslinska, M., Belczynski, K., Klencki, J., & Benacquista, M. 2018, Monthly Notices of the Royal Astronomical Society, 474, 2937, doi: 10.1093/mnras/stx2923
  * Colombo et al. (2022) Colombo, A., Salafia, O. S., Gabrielli, F., et al. 2022, arXiv e-prints.  https://arxiv.org/abs/2204.07592
  * Coughlin et al. (2018) Coughlin, M. W., Tao, D., Chan, M. L., et al. 2018, Monthly Notices of the Royal Astronomical Society, 478, 692, doi: 10.1093/mnras/sty1066
  * Coulter et al. (2017a) Coulter, D. A., Kilpatrick, C. D., Siebert, M. R., et al. 2017a, GCN 21529
  * Coulter et al. (2017b) Coulter, D. A., Foley, R. J., Kilpatrick, C. D., et al. 2017b, Science, 358, 1556, doi: 10.1126/science.aap9811
  * Dewi et al. (2002) Dewi, J. D., Pols, O. R., Savonije, G. J., & Van Den Heuvel, E. P. 2002, Monthly Notices of the Royal Astronomical Society, 331, 1027, doi: 10.1046/j.1365-8711.2002.05257.x
  * Dichiara et al. (2020) Dichiara, S., Troja, E., O’Connor, B., et al. 2020, Monthly Notices of the Royal Astronomical Society, 492, 5011, doi: 10.1093/MNRAS/STAA124
  * Evans et al. (2021) Evans, M., Adhikari, R. X., Afle, C., et al. 2021, arXiv e-prints.  https://arxiv.org/abs/2109.09882
  * Fong & Berger (2013) Fong, W., & Berger, E. 2013, The Astrophysical Journal, 776, 18, doi: 10.1088/0004-637X/776/1/18
  * Fong et al. (2022) Fong, W., Nugent, A., Dong, Y., et al. 2022, arXiv e-prints
  * Fong et al. (2013) Fong, W., Berger, E., Chornock, R., et al. 2013, The Astrophysical Journal, 769, 18, doi: 10.1088/0004-637X/769/1/56
  * Fryer et al. (2006) Fryer, C. L., Young, P. A., & Hungerford, A. L. 2006, The Astrophysical Journal, 650, 1028, doi: 10.1086/506250
  * Gehrels et al. (2006) Gehrels, N., Norris, J. P., Barthelmy, S. D., et al. 2006, Nature, 444, 1044, doi: 10.1038/nature05376
  * Hao & Yuan (2013) Hao, J. M., & Yuan, Y. F. 2013, Astronomy and Astrophysics, 558, 4, doi: 10.1051/0004-6361/201321471
  * Hinshaw et al. (2003) Hinshaw, G., Spergel, D. N., Verde, L., et al. 2003, The Astrophysical Journal Supplement Series, 148, 135, doi: 10.1086/377225
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 99, doi: 10.1109/MCSE.2007.55
  * Im et al. (2017) Im, M., Yoon, Y., Lee, S.-K., et al. 2017, arXiv e-prints.  https://arxiv.org/abs/1710.05861
  * Ivanova et al. (2003) Ivanova, N., Belczynski, K., Kalogera, V., Rasio, F. A., & Taam, R. E. 2003, The Astrophysical Journal, 592, 475, doi: 10.1086/375578
  * Jeong & Lee (2010) Jeong, S., & Lee, C. H. 2010, Journal of the Korean Physical Society, 56, 1619, doi: 10.3938/jkps.56.1619
  * Johnson et al. (2021) Johnson, B. D., Leja, J., Conroy, C., & Speagle, J. S. 2021, The Astrophysical Journal Supplement Series, 254, 22, doi: 10.3847/1538-4365/abef67
  * Kelley et al. (2010) Kelley, L. Z., Ramirez-Ruiz, E., Zemp, M., Diemand, J., & Mandel, I. 2010, The Astrophysical Journal Letters, 725, 91, doi: 10.1088/2041-8205/725/1/L91
  * Leja et al. (2019) Leja, J., Carnall, A. C., Johnson, B. D., Conroy, C., & Speagle, J. S. 2019, The Astrophysical Journal, 876, 3, doi: 10.3847/1538-4357/ab133c
  * Leja et al. (2017) Leja, J., Johnson, B. D., Conroy, C., van Dokkum, P. G., & Byler, N. 2017, The Astrophysical Journal, 837, 170, doi: 10.3847/1538-4357/aa5ffe
  * Levan et al. (2017) Levan, A. J., Lyman, J. D., Tanvir, N. R., et al. 2017, The Astrophysical Journal, 848, L28, doi: 10.3847/2041-8213/aa905f
  * Lien et al. (2014) Lien, A., Sakamoto, T., Barthelmy, S. D., et al. 2014, Proceedings of Science, 829, 1, doi: 10.3847/0004-637x/829/1/7
  * Lipunov et al. (2017) Lipunov, V., Gorbovskoy, E., Kornilov, V., et al. 2017, GCN 21570
  * Madau & Fragos (2017) Madau, P., & Fragos, T. 2017, The Astrophysical Journal, 840, 39, doi: 10.3847/1538-4357/aa6af9
  * Margutti & Chornock (2021) Margutti, R., & Chornock, R. 2021, Annual Review of Astronomy and Astrophysics, 59, 155, doi: 10.1146/annurev-astro-112420-030742
  * McKinney (2010) McKinney, W. 2010, in Proceedings of the 9th Python in Science Conference, ed. S. van der Walt & J. Millman, 51–56, doi: 10.25080/Majora-92bf1922-00a
  * Nakar et al. (2006) Nakar, E., Gal‐Yam, A., & Fox, D. B. 2006, The Astrophysical Journal, 650, 281, doi: 10.1086/505855
  * Nugent et al. (2022) Nugent, A., Fong, W., Dong, Y., et al. 2022, arXiv e-prints
  * O’Connor et al. (2022) O’Connor, B., Troja, E., Dichiara, S., et al. 2022, arXiv e-prints.  https://arxiv.org/abs/2204.09059
  * Ofek et al. (2007) Ofek, E. O., Cenko, S. B., Gal‐Yam, A., et al. 2007, The Astrophysical Journal, 662, 1129, doi: 10.1086/518082
  * Oliphant (2006) Oliphant, T. E. 2006, A guide to NumPy (USA: Trelgol Publishing).  https://web.mit.edu/dvp/Public/numpybook.pdf
  * Pan et al. (2017) Pan, Y. C., Kilpatrick, C. D., Simon, J. D., et al. 2017, The Astrophysical Journal Letters, 848, L30, doi: 10.3847/2041-8213/aa9116
  * Paterson et al. (2020) Paterson, K., Fong, W., Nugent, A., et al. 2020, The Astrophysical Journal, 898, L32, doi: 10.3847/2041-8213/aba4b0
  * Pérez & Granger (2007) Pérez, F., & Granger, B. E. 2007, IEEE Journals & Magazines, 9, 21, doi: 10.1109/MCSE.2007.53
  * Perna et al. (2022) Perna, R., Artale, M. C., Wang, Y.-h., et al. 2022, Monthly Notices of the Royal Astronomical Society, 512, 2654, doi: 10.1093/mnras/stac685
  * Piran (1992) Piran, T. 1992, The Astrophysical Journal Letters, 389, 45, doi: 10.20595/jjbf.19.0_3
  * Price-Whelan et al. (2018) Price-Whelan, A. M., Sipocz, B. M., Günther, H. M., et al. 2018, The Astronomical Journal, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Rastinejad et al. (2022) Rastinejad, J. C., Gompertz, B. P., Levan, A. J., et al. 2022, arXiv e-prints.  https://arxiv.org/abs/2204.10864
  * Robitaille et al. (2013) Robitaille, T. P., Tollerud, E. J., Greenfield, P., et al. 2013, Astronomy & Astrophysics, 558, A33, doi: 10.1051/0004-6361/201322068
  * Safarzadeh & Berger (2019) Safarzadeh, M., & Berger, E. 2019, The Astrophysical Journal, 878, L12, doi: 10.3847/2041-8213/ab24df
  * Safarzadeh et al. (2019a) Safarzadeh, M., Berger, E., Leja, J., & Speagle, J. S. 2019a, The Astrophysical Journal, 878, L14, doi: 10.3847/2041-8213/ab24e3
  * Safarzadeh et al. (2019b) Safarzadeh, M., Berger, E., Ng, K. K. Y., et al. 2019b, The Astrophysical Journal, 878, L13, doi: 10.3847/2041-8213/ab22be
  * Safarzadeh et al. (2018) Safarzadeh, M., Ramirez-Ruiz, E., Andrews, J. J., et al. 2018, The Astrophysical Journal, 872, 105, doi: 10.3847/1538-4357/aafe0e
  * Sana et al. (2012) Sana, H., De Mink, S. E., De Koter, A., et al. 2012, Science, 337, 444, doi: 10.1126/science.1223344
  * Santoliquido et al. (2022) Santoliquido, F., Mapelli, M., Artale, M. C., & Boco, L. 2022, arXiv e-prints.  https://arxiv.org/abs/2205.05099
  * Soares-Santos et al. (2017) Soares-Santos, M., Holz, D. E., Annis, J., et al. 2017, The Astrophysical Journal, 848, L16, doi: 10.3847/2041-8213/aa9059
  * Speagle (2020) Speagle, J. S. 2020, Monthly Notices of the Royal Astronomical Society, 493, 3132, doi: 10.1093/mnras/staa278
  * Tanvir et al. (2017) Tanvir, N. R., Collaboration, A. J. L., & Vinrouge Collaboration. 2017, GCN 21544
  * Tauris et al. (2017) Tauris, T. M., Kramer, M., Freire, P. C. C., et al. 2017, The Astrophysical Journal, 846, 170, doi: 10.3847/1538-4357/aa7e89
  * Tunnicliffe et al. (2014) Tunnicliffe, R. L., Levan, A. J., Tanvir, N. R., et al. 2014, Monthly Notices of the Royal Astronomical Society, 437, 1495, doi: 10.1093/mnras/stt1975
  * Van Der Walt et al. (2011) Van Der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science and Engineering, 13, 22, doi: 10.1109/MCSE.2011.37
  * Vigna-Gómez et al. (2018) Vigna-Gómez, A., Neijssel, C. J., Stevenson, S., et al. 2018, Monthly Notices of the Royal Astronomical Society, 481, 4009, doi: 10.1093/mnras/sty2463
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Wanderman & Piran (2015) Wanderman, D., & Piran, T. 2015, Monthly Notices of the Royal Astronomical Society, 448, 3026, doi: 10.1093/mnras/stv123
  * Yang et al. (2017) Yang, S., Valenti, S., Sand, D., et al. 2017, GCN 21531
  * Zevin et al. (2020) Zevin, M., Kelley, L. Z., Nugent, A., et al. 2020, The Astronomical Journal, 904, 190, doi: 10.3847/1538-4357/abc266
  * Zevin et al. (2019) Zevin, M., Kremer, K., Siegel, D. M., et al. 2019, The Astrophysical Journal, 886, 4, doi: 10.3847/1538-4357/ab498b
  * Zheng & Ramirez-Ruiz (2007) Zheng, Z., & Ramirez-Ruiz, E. 2007, The Astronomical Journal, 655, 1220

  *[EM]: electromagnetic
  *[BNS]: binary neutron star
  *[GRB]: gamma-ray burst
  *[GRBs]: gamma-ray burst
  *[DTD]: delay-time distribution
  *[GW]: gravitational wave
  *[SFR]: star formation rate
  *[SFH]: star formation history
  *[GWs]: gravitational wave
  *[SFHs]: star formation history
  *[DTDs]: delay-time distribution
  *[BNSs]: binary neutron star


structures include the $(k,\delta,\epsilon)$ spiked covariance models and
eigendecay rates given by [Bar+20] such as
$\lambda_{i}=i^{-\alpha}\ln^{-\beta}(i+1)$ for $\alpha=1,\beta>1$. Target data
is generated from random Gaussians with covariances that lead to beneficial
and malignant shifts based on our theories and modifications of the
aforementioned source covariance structures.

All ground truth models are sampled uniformly on the $p$-dimensional
hypersphere, as ${\theta_{\mathsf{s}}^{*}}\sim\mathcal{S}^{p-1}$. Label noise
is sampled as $\varepsilon\sim\mathcal{N}(0,1)$, unless otherwise specified.
For a data matrix $X\in\mathbb{R}^{n\times p}$, training labels are obtained
as $y=X{\theta_{\mathsf{s}}^{*}}+\varepsilon$. Excess risk is computed for
unseen testing data from source and target distributions of interest using
clean labels.

In Figure 4 we take the source to be the $(k,\delta,\epsilon)$ spiked model
with parameters given by $k=70$, $\delta=1$, and $\epsilon=0.005$. The
beneficial shift scales the first $k$ eigenvalues by $\alpha=1.125$ and the
last $p-k$ eigenvalues by $\beta=0.65$. For the malignant shift we use
$\alpha=0.875$ and $\beta=1.35$. The minimum-norm linear interpolator is fit
to 500 data points sampled from a centered multivariate Gaussian with unit
variance and dimension $p=4900$. The model vector is sampled from a centered
Gaussian and scaled to unit norm. The x-axis represents the amount of additive
label noise in training. All evaluation is done on clean data. Each point is
the average of 40 runs.

In Figure 5, we take the source to be the $(k,\delta,\epsilon)$ spiked model
with source parameters as $k=10,\delta=1.0,\epsilon=1e^{-6}$ and target
parameters $\tilde{k}=10,\tilde{\delta}=1.35,\tilde{\epsilon}=6.5e^{-7}$. We
use $n=50$ training data points, $10k$ held-out testing data points in each
OOD test set, and vary $p$ from 75 to 1000 dimensions. We solve OLS using the
closed-form MNI solution on the source data. Each experiment is averaged over
100 independent runs.

In Figure 2 we train fully-connected neural networks with ReLU activation
functions. Data is sampled as above from the covariance structures given by
$\lambda_{i}=i^{-\alpha}\ln^{-\beta}(i+1)$ with varying $\beta$ to obtain
beneficial and malignant shifts. The network architecture is $3$ hidden
layers, with hidden widths $512$ and $2048$. Networks are trained with
stochastic gradient descent with momentum $0.9$ until the training MSE has
reached $<5e^{-6}$. We start with a learning rate of 0.01 and decay by a
stepped cosine schedule for 1,500 epochs. We take batch size of 64 and train
without weight decay. Each experiment is averaged over 20 independent runs. We
train in PyTorch with a single A100 NVIDIA GPU. In these experiments we take
$n=200$ and compare $p=20$ with $p=2000$. Label noise is sampled as
$\mathcal{N}(0,\sigma^{2})$ and we vary $\sigma^{2}$ to show the behaviors at
varying train label noise.

In Figure 8 we train full-connected neural networks with ReLU activation
functions. Source data is sampled from a mean-centered Gaussian with diagonal
covariance matrix with eigenvalues $\lambda_{i}=i^{-1}\ln^{-1.5}(i+1)$. Target
covariate shifts are implemented in the style of Theorem 3.4 where the top $k$
source eigenvalues are multiplied by $\alpha$ and the bottom $p-k$ source
eigenvalues are multiplied by $\beta$. In this experiment, we take
$k=10,\alpha=2,\beta=0.1$ and experiment with $n=400$ source data samples for
$p=200$ and $p=4,000$. The network architecture is 3 hidden layers with hidden
width $2,048$. Our training setup is the same as given above for prior MLP
experiments.

###  I.2 CIFAR-10 and CIFAR-10C experiments

In Figures 3 and 9 we use a binary variant of CIFAR-10 and CIFAR-10C. For
details on the CIFAR-10C dataset, see [HD19]. The binary problem is
constructed by selecting only the dog and truck classes. To stay
overparameterized, we subsample $n=500,1000,2000$ points in a class-balanced
manner. Images are flattened into $p=3072$ dimensional vectors. We fit our
model using the OLS solution for the MNI against $\\{0,1\\}$ class labels. We
test on the same two classes from CIFAR-10 and CIFAR-10C Gaussian blur and
Gaussian noise corruptions. Recall that these two corruptions were selected
for their eigenspectra’s similarity to beneficial and malignant shifts,
respectively. Label noise is injected by flipping class labels with a given
probability.

In Figure 10, we train ResNet18 models on the entire CIFAR-10 dataset and
evaluate on the CIFAR-10C test sets for the Gaussian blur and Gaussian noise
corruptions. The setting is not high-dimensional because we train on 50000
images with 3072 dimensions. However, the ResNet18 architecture has around
11.7 million parameters, so the level of overparameterization is very high.
The training procedure is similar to that used for our MLP experiments.
Networks are trained with stochastic gradient descent with a learning rate of
0.1 and stepped cosine decay schedule for 60 epochs. Each point in the plot is
an average over 30 independent runs. As before, we train in PyTorch with a
single A100 NVIDIA GPU. Label noise takes the form of random label flips with
probabilities 0.1 to 0.9.

##  Appendix J Additional experiments

We present a number of additional supporting experiments that show: (1) more
cases of the behavior of the MNI and MLPs under covariate shift on synthetic
datasets; (2) underparameterized and overparameterized regimes for linear
regression under covariate shift for more realistic eigendecay rates outside
of $(k,\delta,\epsilon)$ spiked covariance models; (3) cases in which the MNI
is overfit in a tempered or catastrophic manner and evaluated on OOD datasets
constructed based on our results in Theorem 3.4, indicating that our insights
hold up for the MNI even when benign source conditions are not satisfied; (4)
the value of overparameterization for the MNI trained on CIFAR-10 and
evaluated on CIFAR-10C; (5) experiments training ResNet-18 models to
interpolation on the full CIFAR-10 dataset and evaluated on CIFAR-10C blur and
noise corruptions.

###  J.1 MNI on Synthetic Data

Figure 4: We fit interpolating linear models to random Gaussian data sampled
from spiked covariance models with parameters $k,\delta,\epsilon$. In this
setting, $k=70,n=500$, $p=4900$, $\delta=1$ and $\epsilon=0.005.$ To
illustrate a beneficial shift, we scale the first $k$ eigenvalues by
$\alpha=1.125$ and the last $p-k$ eigenvalues by $\beta=0.65$. Similarly, for
the malignant shift we use $\alpha=0.875$ and $\beta=1.35$. All experiments
are averaged over 25 independent runs with standard error bars displayed. Note
that the bias is consistently below $10^{-16}$.

In Figure 4, we experiment with interpolating linear models where
$\Sigma_{\mathsf{s}},\Sigma_{\mathsf{t}}$ are given by
$(k,\delta,\epsilon)$-spike covariances with $k=70$, $n=500$, and $p=4900$. We
design problem parameters to show settings in which
$\mathrm{tr}(\Sigma_{\mathsf{t}})>\mathrm{tr}(\Sigma_{\mathsf{s}})$ and we get
a beneficial shift, and
$\mathrm{tr}(\Sigma_{\mathsf{t}})<\mathrm{tr}(\Sigma_{\mathsf{s}})$ and we get
a malignant shift. To do this, the source covariance matrix is constructed
using $\delta=1$ and $\epsilon=0.005.$ To illustrate a beneficial shift, we
scale the first $k$ eigenvalues by $\alpha=1.125$ and the last $p-k$
eigenvalues by $\beta=0.65$. Similarly, for the malignant shift we use
$\alpha=0.875$ and $\beta=1.35$. The resulting plots are significant because
they highlight the distinct effects that the first $k$ and last $p-k$
components have on the excess risk.

As illustrated by our main theorems, increasing the energy of an eigenvalue
has a negative impact on the risk. Nonetheless, these plots show that where
the increase happens plays an important role on how the shift affects
generalization. We are able to improve performance by decreasing the energy on
the tail and increasing the energy on the head in such a way that the total
energy is increased. In short, this setting is a direct connection to our
theory and shows clearly that our constructions for beneficial and malignant
shifts, when mildly overparameterized, hold up in low and high train label
noise regimes, with higher noise exacerbating the effects of the shifts. In
addition, Figure 4 demonstrates that the variance generally contributes much
more significantly to the overall risk.

Figure 5: We experiment with the $(k,\delta,\epsilon)$ spiked covariance
models and examine conditions for beneficial and malignant shifts as given in
Theorem 3.4. We take
$n=50,k=10,\delta=1.0,\epsilon=1e^{-6},\tilde{\delta}=1.5,\tilde{\epsilon}=5e^{-7}$,
and vary $p$. In all cases,
$\mathrm{tr}(\Sigma_{\mathsf{t}})>\mathrm{tr}(\Sigma_{\mathsf{s}})$, showing
that beneficial shifts of this form can occur. As we increase $p$ while
keeping other problem parameters fixed we observe the transition from mild to
severe overparameterization and see the cross-over point between the shift
going from beneficial to malignant. For both ID and OOD excess risk, we
observe that excess risk is a decreasing function of input dimension. Curves
are averaged over 100 independent runs.

Figure 5 shows another example of the transition from mild
overparameterization to severe overparameterization in the case of
$(k,\delta,\epsilon)$ spiked covariance models. In this example we take
$k=10,n=50,\delta=1.0,\epsilon=1e^{-6}$. Using our shifts defined in Theorem
3.4 we set $\alpha=1.5$ and $\beta=0.5$. We plot excess MSE on both ID and OOD
test sets vs. the input data dimension, while holding all other problem
parameters fixed and clearly observe the transition from beneficial to
malignant shifts in keeping with our theorem.

Next, we experimentally show that while our theory is built for benign source
covariance structures it holds for non-benign covariances. In particular, we
examine eigendecay rates that are known to lead to tempered overfitting and
catastrophic overfitting [Mal+22]. [Bar+20] identify the covariance structure
given by $\lambda_{i}=i^{-1}\textrm{ln}^{-2}(i+1)$ as sufficient for benign
overfitting. The rate of $i^{-\alpha}$ for $\alpha>1$ is akin to a ridgeless
Laplace kernel and corresponds to tempered overfitting. Finally, the rate of
$i^{-\ln(i)}$ is akin to a ridgeless Gaussian kernel and corresponds to
catastrophic overfitting. This relative ordering is determined by how high-
dimensional the tail eigenvalues are, in decreasing order.

(a) Ridgeless, $\lambda_{i}=i^{-1}\ln^{-2}(i+1)$ (b) Ridgeless,
$\lambda_{i}=i^{-2}$ (c) Ridgeless, $\lambda_{i}=i^{-\ln(i)}$

Figure 6: Comparing covariate shift in underparameterized vs.
overparameterized linear regression for three different eigendecay rates. In
the $p>n$ setting: (a) leads to benign overfitting, (b) leads to tempered
overfitting, and (c) leads to catastrophic overfitting. We implement simple
multiplicative shifts with $\alpha,\beta$ as defined in Section 3.1 where we
take $n=50,k=10$. Ground truth models are sampled uniformly from
$\mathcal{S}^{p-1}$ and training label noise is sampled from
$\mathcal{N}(0,2)$. Every curve is averaged over 50 independent runs.

It is clear that even though Theorem 3.4 is for the case in which
$\Sigma_{\mathsf{s}}$ satisfies benign source conditions, the style of
beneficial and malignant shift we identify holds for the MNI even when overfit
in a non-benign manner. That is, when $\Sigma_{\mathsf{s}}$ has eigendecay
rates that are tempered or catastrophic we can still obtain non-trivial
beneficial and malignant shifts by changing the energy on the signal and noise
components in a heterogeneous way.

We also notice in Figure 6 that even when varying the dimension up to $p=2000$
at $n=50,k=10$ we do not quite observe the cross-over from beneficial to
malignant shifts in the overparameterized regime. However, we observe that in
Figure 6(a) that the two OOD curves begin to cross-over. Given compute budget,
we run a variant of Figure 6 where we extend up to $p=5000$ and take smaller
$n$, e.g. $n=20,30,40$, in order to closer examine the different regimes of
overfitting. In addition, we experimentally show results for $p=5,10$ which we
liken to the classical linear regression regime in which $k=p<n$, meaning all
of the signal is captured in the $p$ components. In this setting, $\alpha$
shifts are all that influence the distribution shift behavior. We show these
behaviors in Figure 7.

(a) $n=20$ (b) $n=30$ (c) $n=40$

Figure 7: Comparing covariate shift in underparameterized vs.
overparameterized linear regression for when
$\lambda_{i}=i^{-1}\ln^{-2}(i+1)$. We implement simple multiplicative shifts
with $\alpha,\beta$ as defined in Section 3.1 where we take $k=10$ and vary
$n$. Ground truth models are sampled uniformly from $\mathcal{S}^{p-1}$ and
training label noise is sampled from $\mathcal{N}(0,2)$. Every curve is
averaged over 100 independent runs.

###  J.2 MLP on Synthetic Data

We now show additional results for MLPs trained to interpolation on synthetic
datasets. This experiment is analogous to that of Figure 2 except that we
implement shifts in the style of Theorem 3.4.

(a) $n=400,p=200,h=2048$ (b) $n=400,p=4k,h=2048$

Figure 8: We implement multiplicative shifts for interpolating 3-layer ReLU
MLPs in the style of Theorem 3.4. Source data, $X$, is sampled from a mean-
centered Gaussian with diagonal covariance given by
$\lambda_{i}=i^{-1}\ln^{-1.5}(i+1)$. Ground truth models are sampled as
${\theta_{\mathsf{s}}^{*}}\sim\mathcal{S}^{p-1}$ and training label noise is
samples as $\varepsilon=\mathcal{N}(0,\sigma_{x}^{2})$. Noisy training labels
are obtained as $y=X{\theta_{\mathsf{s}}^{*}}+\varepsilon$. The target
covariances are obtained by multiplying the top $k=10$ source eigenvalues by
$\alpha$ and the bottom $p-k$ source eigenvalues by $\beta$. From our theory,
we expect that $\alpha=2,\beta=0.1$ leads to beneficial shifts while
$\alpha=0.1,\beta=2$ leads to malignant shifts. We see this holds up when
$p>n$, and that $h>n$ does not change this relationship. All curves are
averaged over 20 independent runs and each training run reaches MSE loss $\leq
5e^{-6}$.

In Figure 2 we sampled the ID dataset from a mean-centered Gaussian with
diagonal covariance that has eigenvalues $\lambda_{i}=i^{-1}\ln^{-3}(i+1)$ and
we examined the behavior for OOD datasets under covariate shift where the
eigenvalues of the OOD covariance are given by
$\lambda_{i}=i^{-1}\ln^{-2}(i+1)$ and $\lambda_{i}=i^{-1}\ln^{-4}(i+1)$. In
Figure 8 we take the ID data to be sampled from a mean-centered Gaussian with
diagonal covariance that has eigenvalues $\lambda_{i}=i^{-1}\ln^{-1.5}(i+1)$.
For the covariance of the OOD datasets, we shift the top $k=10$ eigenvalues by
a factor of $\alpha$ and the bottom $p-k$ eigenvalues by a factor of $\beta$,
as in the setting of Theorem 3.4. We experiment here with $\alpha=2,\beta=0.1$
and $\alpha=0.1,\beta=2$. Each model achieves training MSE $\leq 5e^{-6}$. We
see the same trends as in Figure 2 with respect to $p>n$ versus $h>n$. In the
$p<n$ case, even though $h>n$ we do not clearly observe a beneficial shift as
predicted by our high-dimensional linear theory. However, when $p>n$ we do
observe beneficial shifts for $\alpha=2,\beta=0.1$, as suggested by our
theorem for the mildly overparameterized case.

###  J.3 MNI on CIFAR-10C Experiments

(a) $p>n$, In-distribution (b) $p>n$, Blur

(c) $p>n$, Noise

(d) $p<n$, In-distribution (e) $p<n$, Blur

(f) $p<n$, Noise

Figure 9: We fit the ridgeless OLS solution to binary CIFAR-10 (dog vs. truck)
and test on binary CIFAR-10C under Gaussian blur and noise corruptions. In the
top row, we vary the level of overparameterization as $n=500,1k,2k$ and
average each curve over 50 independent runs. In this $p>n$ setting the
ridgelss OLS solution results in the MNI. In the bottom row we obtain a non-
interpolating, ridgeless linear solution. Evaluations are done on severity 3
of CIFAR-10C, however the results hold up across all severities. We plot
excess classification error vs training label noise, which is class label flip
probability. We see that overparameterization improves robustness of the MNI
at all noise levels.

In Figure 9 (a-c) we show that overparameterization improves OOD excess
classification error for the MNI fit to binary CIFAR-10 and evaluated on
binary CIFAR-10C under Gaussian blur and noise corruptions. The details of
these datasets and setups are given in Appendix I. We note that all of the
curves in the top row of this figure are in the overparameterized regime,
meaning they are on the right side of the double descent curve. Flattened
CIFAR images have $p=3,072$ and so we vary the number of training subsample
sizes over $n=500,1000,2000$ in order to remain in an overparameterized
setting. We find that when we are overparameterized, as we reduce $n$ we
obtain improved performance. We average over 50 independent runs in each
setting and provide standard error bars to show that this observation is not
due to specific random samples. We also see that at higher levels of
overparameterization, the relative difference in excess classification error
between ID, blur, and noise test sets lessens. For example, at $0.0$ label
noise and $n=2k$ the average excess error varies from $0.3028$ on the blur set
to $0.4668$ on the noise set for an absolute difference of $0.164$, whereas at
$n=500$ the average excess error only varies from $0.252$ on the blur set to
$0.3131$ on the noise set for an absolute difference of $0.0611$.

For completeness, in Figure 9 (d - e) we show the above setting in the
underparameterized regime where we obtain the linear solution via the
ridgeless OLS solution. As these plots are on the left side of the double
descent peak, we see that adding more data improves OOD excess classification
error. While these models are not interpolating, we observe that noise
corruptions lead to nearly catastrophic performance, meaning random guessing,
on the OOD test sets, whereas blur corruptions lead to more benign
performance. Finally the ID performance appears to be tempered, in showing a
nearly linear relationship between train label noise and test excess
classification error.

###  J.4 ResNet on CIFAR-10C Experiments

(a) Interpolating ResNet tested on Blur (b) Interpolating ResNet tested on
Noise

Figure 10: We train ResNet18 on clean CIFAR-10 and evaluate on test sets that
has been corrupted by Gaussian blur and Gaussian noise, which correspond to
beneficial and malignant shifts, respectively. Labels are flipped with
probability 0.1 through 0.9, seen on the x-axis. The setting is not high-
dimensional because the training data contains 50000 images, each of which are
3072-dimensional. ResNet18 contains around 11.7 million parameters, so the
setting is very overparameterized. We observe that while both shifts
negatively affect generalization, the beneficial shift isn’t as bad as the
malignant shift. This result is similar to those seen in subfigures (a) and
(b) in Figure 2, where the data is not high-dimensional but the MLP is
overparameterized.

Figure 10 shows the behavior of interpolating ResNets trained on the full
CIFAR-10 dataset and evaluated on CIFAR-10C blur and noise corruptions. While
these numbers are suboptimal with respect to CNNs on CIFAR-10 we note that
they are justified in our setting as our goal is to study interpolating
models. At 90% label noise it takes a lot of compute to interpolate the entire
CIFAR-10 dataset, especially if using data augmentations, weight decay, or
other regularizations. As such, we turn off weight decay and data
augmentations for these models to be able to tractably interpolate CIFAR-10 at
high noise levels.