# Predicting Short-duration GRB Rates in the Advanced LIGO Volume

Tzvetelina A. Dimitrova Earth and Space Exploration, Arizona State University,
P.O. Box 871404, Tempe, AZ 85287-1404, USA  Nathaniel R. Butler Earth and
Space Exploration, Arizona State University, P.O. Box 871404, Tempe, AZ
85287-1404, USA  Srihari Ravi  Earth and Space Exploration, Arizona State
University, P.O. Box 871404, Tempe, AZ 85287-1404, USA

###### Abstract

Starting with models for the compact object merger event rate, the short-
duration Gamma-ray Burst (sGRB) luminosity function, and the Swift/BAT
detector, we calculate the observed Swift sGRB rate and its uncertainty. Our
probabilistic sGRB world model reproduces the observed number distributions in
redshift and flux for 123 Swift/BAT detected sGRBs and can be used to predict
joint sGRB/LIGO detection rates. We discuss the dependence of the rate
predictions on the model parameters and explore how they vary with increasing
experimental sensitivity. In particular, the number of bursts in the LIGO
volume depends strongly on the parameters that govern sGRB beaming. Our
results suggest that nearby sGRBs should be observed to have broader jets on
average ($\theta_{\rm jet}\gtrsim 30$ degrees), as compared to the narrowly-
beamed appearance of cosmological sGRBs due to detection selection effect
driving observed jet angle. Assuming all sGRBs are due to compact object
mergers, within a $D<200$ Mpc aLIGO volume, we predict $0.18^{+0.19}_{-0.08}$
sGRB/GW associations all-sky per year for on-axis events at Swift
sensitivities, increasing to $1.2^{+1.9}_{-0.6}$ with the inclusion of off-
axis events. We explore the consistency of our model with GW170817/GRB 170817A
in the context of structured jets. Predictions for future experiments are
made.

##  1 Introduction

Knowledge of Gamma-Ray Bursts (GRBs) has greatly benefited from the dramatic
success of the Swift satellite (Gehrels et al., 2005), capable of rapid
localization and highly-sensitive multi-wavelength observations that have
enabled unprecedented ground-based followup and a high-rate of GRB redshift
determinations. However, there remains relatively few short-duration GRB
(sGRB) events with redshift measurements, and relatively little is understood
about the properties of these populations (for an overview, see; Berger,
2014).

The probable origin of sGRBs from mergers between binary compact objects such
as binary neutron stars (BNS) and neutron star black hole (NSBH) binaries make
them promising gravitational wave (GW) electromagnetic (EM) counterparts
capable of detection by advanced interferometers. Following the success of the
Laser Interferometer Gravitational-Wave Observatory (LIGO) in detecting
GW150914, interest spiked in the sGRB population. The advanced aLIGO-Virgo
discovery of GW170817 – near-simultaneous with GRB 170817A detected by the
Fermi GRB Monitor (GBM) and International Gamma-Ray Astrophysics Laboratory
(Abbott et al., 2017a) – provides direct evidence for BNS mergers as sGRB
progenitors.

Here we seek to determine whether and how the cosmological sGRB population is
consistent with the source progenitor merger population constrained by LIGO.
We employ a sample of 123 sGRBs detected by the Burst Alert Telescope (BAT)
onboard Swift, and we carefully treat the Swift/BAT detection threshold to
account for the majority of cosmological sGRBs which are undetected. We assume
a universal model for the merger event rate, along with parameterized
luminosity functions and beaming distributions, in order to better understand
the sGRB population. The derived probabilistic sGRB world model successfully
reproduces the observed Swift/BAT number distributions in $z$ and flux, and it
allows for predictions of the joint GW/sGRB detection rates for both on- and
off-axis sGRBs. While the GW signal does not depend on beaming, the inclusion
of distributions to describe sGRB beaming – and the resulting number of events
pointed toward the observer – are essential to establishing a rate
correspondence.

To set the stage, we carry-out an initial, approximate calculation. Swift
detects about 2 sGRBs per year within 1 Gpc, based on 7 sGRBs with measured
redshift over 15.7 years and a 20% rate of redshift determination. This
corresponds to a mean, all-sky rate density for cosmological sGRBs of $20\pm
6$ Gpc-3 yr-1. Two aLIGO BNS events in two years of observation (within 200
Mpc; The LIGO Scientific Collaboration et al., 2021) implies a rate of 125
Gpc-3 yr-1, assuming constant density. Roughly then, we expect an all-sky,
joint detection rate of $0.16\pm 0.05$ per year.

Separating sGRBs from their progenitor path is difficult (e.g., Sarin et al.,
2022), and we will consider a rate model including both BNS and NSBH mergers.
Due to a large uncertainty in the underlying mass ratios, The LIGO Scientific
Collaboration et al. (2021) estimates a possible rate of BNS events as large
as $\approx 2000$ Gpc-3 yr-1, which dominates among the progenitor channels.
This provides some flexibility for the source sGRB rate to be $\gtrsim 100$
times the observed rate. Assuming that all sGRBs are due to compact object
mergers, that factor will limit the extent to which the sGRB population – not
all of which are above detection threshold – can be beamed. We estimate a
detection fraction $\lesssim 25$% within 1 Gpc, constraining our fits to yield
a mean $\theta_{\rm jet}>17$ degrees (Section 4.3, below). We find that the
observed jet angle is driven by the detection selection effect, where bias
favoring bright events at increasing redshift makes the cosmological
population appear quite narrowly-beamed, whereas events within the aLIGO’s
lower distance volume are predicted to allow a wider mean $\theta_{\rm
jet}\gtrsim 30$ degrees to be seen.

###  1.1 Previous Estimates

Prior estimates of the GW/sGRB joint detection rate span a large range.
Partially, the variation is due to pre-LIGO estimates without the benefit of
LIGO constraints. It is also due in large part to differing assumptions
regarding the extent of sGRB beaming.

Guetta & Piran (2005) consider sGRB/BNS associations and the effect of time-
delay until merger on the redshift distribution and luminosity function and
find a local sGRB rate of $\sim 0.8$ Gpc-3 yr-1. Liu & Yu (2019) find a local
sGRB rate of $\sim 3-4$ Gpc-3 yr-1. Guetta & Piran (2006) find a much higher
rate of $\sim 8-30$ Gpc-3yr-1 from redshift and luminosity distributions
constrainted by the Swift/HETE-II sample. Guetta & Stella (2008) predict that
a large fraction of detectable GW events will coincide with sGRBs. Nakar et
al. (2006) find a higher local rate as well of $>10$ Gpc-3 yr-1. Coward et al.
(2012) account for dominant detection biases to find an sGRB rate density of
$\sim 8^{+5}_{-3}-1100^{+700}_{-470}$ Gpc-3 yr-1 out to $z\approx 0.5$,
assuming isotropic emission and correcting for beaming. Assuming that all BNS
mergers produce an sGRB, Siellez et al. (2013) estimate a simultaneous BNS-GW
observation rate of $0.11-4.2$ yr-1 with Swift and Fermi. Under the assumption
that a notable fraction of BNS mergers result in an sGRB, Coward et al. (2012)
find a lower and upper detection rate limit for a aLIGO-Virgo search of
$1-180$ yr-1.

Other recent studies attempt to also constrain sGRB beaming and present rates
using these assumptions. Ghirlanda et al. (2016) estimate an sGRB rate within
the aLIGO volume of $0.007-0.03$ yr-1 with jet opening angles between 3 and 6
degrees. Liu & Yu (2019) estimate a similar rate of Swift detectable sGRBs in
the aLIGO horizon of $0.032$ yr-1. Sarin et al. (2022) find a BNS merger rate
of $384^{+431}_{-213}$ Gpc-3 yr-1, with an average $\theta_{\rm jet}\sim 15$
degrees of BNS produced sGRBs. They report that an estimated $40\%$ of BNS
mergers produce jets, and measure a fraction $0.02^{+0.02}_{-0.01}$ of BNS
events and $0.01\pm 0.01$ NSBH mergers to result in observable sGRBs. Jin et
al. (2018) use $\theta_{\rm jet}\sim 0.1$ rad and a local BNS rate density of
$\sim 1109^{1432}_{-657}$ Gpc-3yr-1, which decreases to $\sim 162^{140}_{-83}$
Gpc-3yr-1 when excluding a narrowly beamed sGRB. They find that off-axis sGRB
events enhance the GW/sGRB rate to a $10\%$ association probability.

Figure 1: The separation of the GRB sample into short ($T_{\rm 90}<2$s) and
long-duration class, including events with and without measured $z$.

##  2 Data Retrieval and Sample Selection

We use publicly-available GRB data from the Swift BAT Gamma-Ray Burst
Catalog111https://swift.gsfc.nasa.gov/results/batgrbcat (Lien et al., 2016)
for 1389 GRBs detected between December 17, 2004 and August 29, 2020. We
retrieve burst duration intervals $T_{\rm 90}$, signal-to-noise ratios $SNR$,
partial coding fractions $f_{\rm p}$, redshifts $z$ when available, and energy
fluxes $F$ (15-150 keV band). A detailed description of these quantities and
relevant discussion may be found in the BAT catalog or in Butler et al.
(2007). We apply a cutoff $T_{\rm 90}<2$ s to define the short-duration burst
population (e.g., Kouveliotou et al., 1993, ; also, Figure 1). To determine
rest-frame luminosities, we assume a cosmology with h $=0.71$, $\Omega_{\rm
m}=0.3$, and $\Omega_{\rm\lambda}=0.7$. All error regions below correspond to
the 90% confidence intervals, unless otherwise noted.

Figure 2: Left: Observed $SNR$ values (blue points) – after correction to a
fixed integration time and detector area (orange points) – correlate tightly
with the observed 15-150 keV flux. The $SNR$ histogram is shown with a fitted
curve that yields a detection threshold in $\log(SNR)$ of $0.9\pm 0.1$ (black;
blue dashed). Right: The observed number of sGRBs over a range of Fluence $S$
for Swift ($15-150$ keV band; blue) and Fermi ($10-10^{3}$ keV band; green).
The best-fit GRB world models for Swift (orange) and for Fermi (orange dashed)
derived below are also plotted, including an $S^{0.5}$ bandpass correction
factor for Swift.

A plot of the number of events versus a certain quantity can provide important
lessons. The number of events with $SNR$ above a given $SNR$, is quite steep,
with a slope $<-2$ (Figure 2, Left). This may suggest a steep luminosity
function, but it certainly indicates – as is well-known – that sGRBs tend to
be weakly detected. The number of events with fluence (or flux) above a
certain level, the so-called logN-logS relation, can be used to to infer the
nature of the source population. A logN-logS slope of $-3/2$ for a relatively
nearby population (such as sGRBs) is only weekly-dependent upon the luminosity
function and is indicative of a non-evolving source population (e.g., Longair,
1966). Such a slope appears to be present at the bright end of the logN-logS
relation (Figure 2, Right). The faint-end slope for sGRBs may also be near
$-3/2$ after accounting for threshold effects (e.g., Petrosian & Lee, 1996).
However, to disentangle the physical effects that drive the number counts – to
separately measure the luminosity function and the number density versus
redshift – redshift measurements and a careful treatment of the detection
limit are required. Parametric modelling is beneficial given the typically low
$SNR$ values.

In order to use the measured data to constrain intrinsic GRB properties, we
need a prescription for mapping rest frame luminosities to observer frame
fluxes and for deciding whether the observed fluxes are above the BAT
detection threshold. We find that $SNR$ for Swift can be predicted from the
flux $F$ with 0.1 dex scatter (i.e., $1\sigma$ spread; see, Figure 2, Left),
as:

$SNR=(F/F_{\circ})\sqrt{T_{\rm 90}f_{p}}/\sqrt{2950+F/F_{\circ}}.$ (1)

Here, the zero point $F_{\circ}=1.2\times 10^{-10}$ erg/cm2/s, and 2950
represents a typical background count rate. We note that the observed $SNR$
distribution is very steep ($dN/dSNR\propto SNR^{-3}$) – indicating relatively
little dynamic range in the sample – and is well-fit by a powerlaw truncated
to zero at threshold and convolved with a Gaussian. This yields a stochastic
$SNR$ threshold of $10^{0.9\pm 0.1}$.

Figure 3: Left: Binary merger rates as a function of $z$ from Safarzadeh et
al. (2019) for different minimum time delays $\tau_{\rm min}$ between
formation – which is assumed to track the SFR and Universal volume – and
merger. Right: The $E_{\rm iso}$ luminosity function – found by integrating
Equation 3 for the $L_{\gamma}$ luminosity function over $f_{b}$ – includes
the detector threshold which increases with increasing $z$. The asymptotic
power law behavior is shown with dotted lines. The integral of the model (see
“norm” in legend) provides the fraction of sGRBs detected.

We assume flux $F$ and fluence $S$ in the Swift/BAT bandpass are related as
$S=FT_{\rm 90}$. Due to the low $SNR$ values and the soft Swift bandpass,
spectral parameters are typically poorly-constrained and the conversion to
bolometric $S$ is challenging. To simplify matters, we assume $S_{\rm
bol}=xS$, where $x$ is a bandpass conversion factor applied to all bursts. We
expect a large uncertainty of about a factor of 2 (0.3 dex) for $x$. We assume
a mean correction $x=3.77$ which is suitable for a typical, soft Swift sGRB
spectrum ($E_{\rm peak}\approx 250$ keV). A change to this value would alter
the energetics scale (e.g., shift locations on some plots below).

This simple approach to a bolometric correction appears to break down for
bright sGRBs (see, Figure 2, Right). In comparison with sGRB fluences from
Fermi/GBM (Gruber et al., 2014; von Kienlin et al., 2014; Bhat et al., 2016;
von Kienlin et al., 2020), we find that a larger correction is required for
bright bursts. For bursts with $S>10^{-7}$ erg/cm2 (15-150 keV), we assume
that the bolometric correction factor increases as $S^{0.5}$. This observer-
frame effect likely arises due to the increasing hardness of bright bursts and
the resulting lack of flux in the soft 15-150 keV band as compared to the
harder, and approximately-bolometric 10-1000 keV Fermi/GBM band. Although the
effect is confirmed to correlate with hardness in the Swift 50-150 keV to
15-50 keV bands for 12 bursts measured in common between Swift and Fermi, we
model it indirectly (on $S$ alone) because a correction using the Swift
hardness would introduce large errors and require spectral modelling not
possible (with sufficient accuracy) given only the soft Swift bandpass. While
a bandpass correction based on $S$ only would not be an acceptable approach
for individual burst modelling, it does capture the ensemble bandpass effects
well. From Figure 2 (Right), we see that both Swift and Fermi have similar
sGRB sensitivities (see also, Band, 2003), and both have logN-logS
distributions scaling roughly as $S^{-3/2}$ at the bright end.

Finally, we note that there are 3 events with anomalously high luminosities
($>10^{52}$ erg/s) which we exclude from the analysis, but return to in
Section 3.3 below. GRBs 090426 and 111117A have spectroscopic $z>2$ and may
lie on the short-duration end of the long-duration GRB population (see,
Antonelli, L. A. et al., 2009; Levesque et al., 2010; Selsing, J. et al.,
2018). GRB 120804A has an uncertain photometric $z=1.3^{+0.3}_{-0.2}$ (1
sigma) based on a study of the host galaxy (Berger et al., 2013). Our final
sample includes 123 sGRBs between and including GRB 050502 and GRB 20022A. Of
these, 27 (22%) have measured redshifts in the catalog.

##  3 sGRB World Model Summary

Following Butler et al. (2010), to model the observed rates (number versus
luminosity, number versus redshift $z$, number versus beaming fraction, etc.),
we define functional forms for the probability distributions and, for a given
set of the parameters, use the intrinsic model to generate sGRBs reaching the
detector. Not all of these will be detected, and evaluation of the extent to
which the predicted rates from a given set of parameters match the observed
rates requires the detector model above (e.g., Figure 2).

We assume a sGRB rate density model $\dot{\rho}(z)$ due to compact object
mergers from Safarzadeh et al. (2019), which is a convolution of the star
formation rate (SFR) and a merger delay time $\tau_{\rm d}$ distribution. To
convert from rate density to number $N$ per redshift interval:

$\vspace{-0.05in}P_{z}(z)={dN\over dz}=\dot{\rho}(z){dV\over dz}{1\over 1+z}.$
(2)

Here, $V$ is the Universal volume at $z$, and the factor of $1+z$ converts
rates at $z$ to now. We assume a $\tau_{\rm d}$ distribution
$\propto\tau_{d}^{-1}$, where the maximum delay time is 10 Gyr and the minimum
delay time ($\tau_{\rm min}$) is a free-parameter. In principle, the $-1$
slope can also be varied; however this provides no additional model
flexibility over the $z$ range under study here. Resulting $dN/dz$
distributions are shown in Figure 3 (Left). The rate density tracks the SFR
for $\tau_{\rm min}\lesssim 100$ Myr and peaks at lower $z$ for larger
$\tau_{\rm min}$.

The sGRB luminosity function is generally modelled as a single or broken power
law, similar to studies of long-duration GRBs (see, e.g., Guetta & Piran,
2006; Salvaterra et al., 2008; Virgili et al., 2011; D’Avanzo et al., 2014;
Ghirlanda et al., 2016). To get the number of GRBs over a range of
$\gamma$-ray Luminosity $L_{\gamma}$, we assume $P_{L}=dN/dL_{\gamma}$, and

$\begin{split}E_{\circ}{dN\over
dL_{\gamma}}=A(L_{\gamma}/E_{\circ})^{\alpha},\;\;\;\;E_{\rm
min}<L_{\gamma}<E_{\circ}\\\
=A(L_{\gamma}/E_{\circ})^{\beta}.\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;L_{\gamma}\geq
E_{\circ}\end{split}$ (3)

The luminosity function is defined above a minimum $E_{\rm min}$, allowing for
it to be normalized (fixing $A$) assuming $\alpha\geq\beta$ and $\beta<-1$.

Choices for physical quantities to associate with $L_{\gamma}$ include the
total $\gamma$-ray energy reservoir, otherwise known as the beaming-corrected
energy release $E_{\gamma}$, or the energy per solid angle as represented by
the isotropic-equivalent energy release $E_{\rm iso}$. Here,
$E_{\gamma}=f_{b}E_{\rm iso}$, for beaming fraction $f_{b}=1-\cos(\theta_{\rm
jet})$. Also, we are assuming the simple top-hat GRB jet model (e.g., Rhoads,
1999; Frail et al., 2001), where the GRB is detected only if it launches a jet
pointed in the direction of the observer. The rate of such events is
proportional to the fraction of sky illuminated (i.e., $f_{b}$), establishing
a direct relationship between outflow collimation and rates (see, also,
Berger, 2014).

To allow some flexibility between these possibilities, we write
$L_{\gamma}=E_{\rm iso}f_{b}^{c}$. The index $c$ represents a potential
correlation between $f_{b}$ and $E_{\gamma}$. We allow $c$ to vary between 0
and 1, allowing the luminosity argument to vary between $E_{\rm iso}$ ($c=0$)
and $E_{\gamma}$ ($c=1$). Viewed differently, the luminosity function is
allowed to regulate GRBs in terms of total energetics ($c=1$) or local
properties ($c=0$). We can motivate allowing this flexibility by noting that,
during a GRB, the outside of the jet is not in causal contact with the inside
of the jet; the processes that govern the observed GRB brightness distribution
could depend weakly on the global energetics.

Since detection depends on $E_{\rm iso}$, the observed $f_{b}$ distribution
will change with $z$ for $c>0$. With $c=0$, the luminosity function depends on
$E_{\rm iso}$, and the observed $f_{b}$ will take on an average value for all
$z$. We assume a powerlaw $P_{f_{b}}=dN/df_{b}$ distribution (as in, e.g.,
Frail et al., 2001), with

$dN/df_{b}=Bf_{b}^{n}$ (4)

The normalization $B$ is set by integrating between a minimum $f_{b,{\rm
min}}=10^{-4}$ (corresponding to $\theta_{\rm jet,min}=0.8$ degrees) and a
maximum $f_{b}=1$. We apply a weak constraint on the beaming distribution:

$n<-1.5-c\;(1+\alpha).$ (5)

This prevents a flat, observed $\theta_{\rm jet}$ distribution and corresponds
to the statement “short-duration GRBs are beamed.”

We assume a log-normal distribution of $\log(T_{\rm 90})=-0.5\pm 0.4$ and an
exponential distribution for $\log(f_{p})$ with mean -0.23 as in Butler et al.
(2010).

Combining the above distributions, the predicted, observed event rate density
$r$ is:

$\begin{split}r=r_{\circ}f_{b}P_{f_{b}}(f_{b}|n)P_{z}(z|\tau_{\rm
min})P_{L}(L_{\gamma}|\alpha,\beta,c,E_{\rm min},E_{\circ})\\\
\times\;H(L_{\gamma}-L_{\rm th})P_{th}(L_{\rm t}|L_{\gamma,{\rm
min}}).\end{split}$ (6)

Here, $L_{\gamma}=4\pi D_{L}^{2}/(1+z)xSf_{b}^{c}$, and the minimum detectable
$L_{\gamma,{\rm min}}=4\pi D_{L}^{2}/(1+z)xS_{\rm min}(T_{\rm
90},f_{p})f_{b}^{c}$, and $H$ is the Heaviside step function. The threshold
$L_{\rm th}$ distribution $P_{th}(L_{\rm th})$ is assumed to be log-normal
with a mean of $L_{\gamma,{\rm min}}$ and a standard deviation of 0.35 dex.
This accounts for the uncertain conversions from limiting $SNR$ to limiting
Swift 15-350 keV $S$ (0.14–0.18 dex, depending on whether we fix $T_{\rm 90}$
and $f_{b}$ or average over the observed values) and then from Swift $S$ to
bolometric $S$ (0.3 dex). The contribution to the scatter from averaging over
$T_{\rm 90}$ and $f_{p}$ is negligible, and we can effectively marginalize
over those distributions when calculating the distribution normalization by
assuming their mean values.

###  3.1 Model Normalization

We assume all sGRBs are due to compact object mergers. The normalization
$r_{\circ}$ in Equation 6 includes 3 factors, 1 fixed and 2 variable which
become model parameters requiring priors. First is the fixed $\Omega/(4\pi)$
fraction of sky covered by Swift, for which we assume $\Omega$ = 1.3 str,
multiplied by a Swift observation period of 15.7 years. Second is the rate of
merger events in the $D<200$ Mpc aLIGO volume. This normalization factor is
represented by a gamma function prior corresponding to 2 BNS events (GW170817
and GW190425, Abbott et al., 2017b, 2020) detected in 2 years (O1 through O3).
Third is a factor that scales the $D<200$ Mpc yearly rate out to a $D=1$ Gpc
yearly rate appropriate for sGRBs. The third normalization factor is highly-
uncertain based on LIGO team estimates (The LIGO Scientific Collaboration et
al., 2021) and is assumed to be distributed uniformly between 1 and 10. We
choose a distance of 200 Mpc as a representative LIGO sensitivity volume for
BNS events; LIGO is expected to detect BNS events out to 190 Mpc in the O4
observing run (2023–2024) and out to 325 Mpc in O5
(2026–2028)222https://observing.docs.ligo.org/plan.

Figure 4: The observed (blue) and predicted (orange) rates of sGRBs with
respect to $z$ (Left), $E_{\rm iso}$ (Middle), and $S$ (Right). The right
panel includes all sGRBs in the blue histogram and also displays the sGRBs
with measured $z$ as a superimposed, unshaded orange histogram. GRB 150101B is
discussed in Section 4.1.

###  3.2 The Distribution of $E_{\rm iso}$, given $z$

To fit the data, we require $P_{E}$: the distribution of $E_{\rm iso}$, given
$z$. We construct $P_{E}$ by multiplying the luminosity function by the
beaming distribution, and the detection distributions by $f_{b}$, evaluating
at $L_{\gamma}=f_{b}^{c}E_{\rm iso}$, and integrating over $f_{b}$. The
integration over $f_{b}$ results in a smoothly-broken powerlaw with low-energy
slope $\alpha$ and high-energy slope $\beta$ (see, Figure 3; Right). The
transition between slopes occurs over an extended range in $E_{\rm iso}$ from
$E_{\circ}$ until $E_{\circ}/f^{c}_{b,{\rm min}}$. To be explicit, starting
from Equation 6:

$r=r_{\circ}P_{z}(z|\tau_{\rm min})P_{E}(E_{\rm iso}|z,n,\alpha,\beta,c,E_{\rm
min},E_{\circ}),$ (7)

with

$P_{E}=\int f_{b}P_{f_{b}}P_{L}(E_{\rm iso}f_{b}^{c})H(E_{\rm
iso}f_{b}^{c}-L_{\rm thresh})P_{t}df_{b}.$ (8)

We first calculate $P_{E}$ analytically at $z=0$ (i.e., all events detected),
and ensure that it is normalized. The detection distributions $HP_{t}$ are
then applied as a $z$-dependent cutoff at $L_{\rm thresh}=E_{\rm
iso}f_{b}^{c}$, followed by a Gaussian convolution to get the observed
distribution in $E_{\rm iso}$ at a given $z$.

From $P_{E}$ and $P_{z}$, it is then possible to calculate the joint
distribution of $S$ and $z$ and the marginal distributions for $E_{\rm iso}$,
$S$, or $z$. The integrations over $E_{\rm iso}$, $S$, and $z$ are carried out
numerically over fine grids, with the $z$ grid extending from $z=7.1\times
10^{-4}$ (i.e., $D=3$ Mpc) to $z=10$. Figure 3 (Right) displays the $E_{\rm
iso}$ luminosity function at different $z$ values.

###  3.3 Model Fitting

In fitting the above models by matching the model and observed rates, we treat
the measured values of $z$, $T_{\rm 90}$, $f_{p}$, and $S$ (and $E_{\rm iso}$
fixed by $z$ and $S$) as point estimates only (i.e., having no error). The
uncertainty in these quantities is placed in the models, which can be
considered smoothed to avoid over-fitting.

The parameters $\vec{\theta}$ are $n$, $c$, $E_{\rm min}$, $\alpha$, $\beta$,
$E_{\circ}$, and $\tau_{\rm min}$. The observables $\vec{D}$ are $z$, $S$,
$T_{\rm 90}$, and $f_{p}$. The rate evaluated for a given event $i$ is written
$r_{i}$. In the case of GRBs with no measured $z$, we integrate $r_{i}$ over
$z$. To find the parameters that best fit the model, we maximize the Poisson
likelihood:

$\mathcal{L}(\vec{D}|\vec{\theta})=\exp\left[-\int
rd\vec{D}\right]\prod_{i=1}^{N}r_{i}(\vec{\theta})$ (9)

The integral in the exponential in Equation 9, with $d\vec{D}=dz\,dS\,dT_{\rm
90}\,df_{p}$, yields the total number of predicted events. Prior to crawling
the parameter space $\vec{\theta}$ via Markov Chain Monte Carlo (MCMC)
maximization of $\log(\mathcal{L})$, we first marginalize analytically over
$r_{\circ}$ (see, Equation 7).

The MCMC analysis is carried out with PyMC3 (Salvatier et al., 2016), using
multiple chains at a variety of starting parameter values to establish
uniqueness of the derived solution. We assume uniform prior distributions on
all parameters except for $E_{\rm min}$ and $E_{\circ}$. Also included are the
luminosity function slope constraints (Section 3) and constraints on the
$f_{b}$ distribution (Equation 5). For $E_{\rm min}$ and $E_{\circ}$ we assume
log-normal priors, $\log(E_{\rm min})=47.8\pm 0.5$ and
$\log(E_{\circ})=50.5\pm 0.3$. The $\log(E_{\rm min})$ prior is weak and non-
constraining. The tighter $E_{\circ}$ prior is designed to test for the
possibility of a break in the luminosity function.

Figure 5: The all-sky predicted number of sGRBs per year as a function of distance, detectable by Swift. On-axis events are displayed in green, and on- and off-axis bursts are displayed in black. Increased detector area (blue and orange curves) strongly-affect the predicted rate at $D\gtrsim 1$ Gpc. Structured jets (dashed green curve) are discussed in Section 5. Table 1: MCMC Parameters 90% Confidence Intervals Parameter | Value (all GRBs)
---|---
$f_{b}$ index $n$ | $-0.7\pm 0.3$
correlation index $c$ | $>0.6$
$\log(E_{\rm min})$ | $47.3\pm 0.4$
$\alpha$ | $-2.1\pm 0.2$
$\beta$ | $-2.5^{+0.3}_{-2.1}$
$\log(E_{\circ})$ | $50.5^{+0.6}_{-0.3}$
$\tau_{min}$ |  $<700$ Myr

The posterior parameter distributions are estimated by drawing $10^{4}$
parameter $\vec{\theta}$ samples from the MCMC, after thinning by a factor of
4 to limit draw-to-draw correlations. These are summarized in Table 1 and
discussed more in Section 4.3 below. We generate posterior distributions for
the observables $\vec{D}$ by averaging their distributions from Section 3.2
over the sampled parameters.

##  4 Discussion

###  4.1 Predicted Observable Distributions

Predicted distributions and uncertainties for $z$, $E_{\rm iso}$, and $S$ are
overplotted on the observed data in Figure 4. The models fit the observed data
well. The best-fit fluence model (Figure 4; Right), for example, has a reduced
$\chi^{2}_{\nu}=1.1$ for $\nu=4$ degrees of freedom, using the Gehrels (1986)
approximate for the standard deviation of Poisson counts.

While the observed redshift distribution peaks at $z\approx 0.2$, this is
largely due to the decline in detection efficiency at high-$z$ (Figure 3;
Right). The input redshift distribution – Equation 2 averaged over draws of
the $\tau_{\rm min}$ parameter – peaks at $z=1.35$, as compared to $z=1.45$
for the SFR. We estimate a modest, 73% confidence that the rate model is
delayed with respect to the SFR by counting the fraction of draws with a peak
below $z=1.4$.

The predicted distributions for $E_{\rm iso}$ and fluence $S$ are also closely
consistent with the data. However, we note that in the middle panel of Figure
4, there is a possible outlier. GRB 150101B was an un-triggered event also
detected by Fermi/GBM with a significantly higher reported flux (leading to
$E_{\rm iso}$ = $1.1\times 10^{49}$ erg instead of $1.6\times 10^{47}$ erg
(Burns et al., 2018). GRB 150101B may have also been an off-axis event (Troja
et al., 2018).

In the fitting above, we have assumed that events with $z$ should be fit
alongside events without redshift. This can be justified by the similarity of
the observed $S$ distributions in Figure 4 (Right). From a two-sample
Kolmogorov-Smirnov test, we find that the distributions cannot be
significantly differentiated ($\tau=0.168$; p-value 0.58).

Figure 6: The number of on-axis sGRBs expected within the aLIGO volume,
displayed versus the different parameters. Quadratic curves are plotted
through the parameter samples to highlight possible trends. Overall
distributions are displayed as histograms on both axes (see also, Table 1).

The evidence for a break in the luminosity function is formally strong (99%
confidence), given the confidence intervals for $\alpha$ and $\beta$. However,
the measurement of $E_{\circ}$ without the restrictive prior is highly
uncertain given the smooth shape of the luminosity function that results from
the uncertainty in our $E_{\rm iso}$ measurements and the additional smoothing
caused by integration over $f_{b}$. Moreover, it is important to note that the
possible measurement of a break depends on the sample selection. If we include
the 3 high-luminosity events in our fitting (GRBs 090426, 111117A, and
120804A; Section 2), the slopes become entirely consistent
($\alpha=-2.1^{+0.2}_{-0.1}$, $\beta=-2.2\pm 0.1$) with a single power-law
luminosity function. Also, including these events lowers the upper limit on
$\tau_{\rm min}$ (30 versus 700 Myr). The other model parameters in Table 1
are not significantly altered.

###  4.2 GRB Rate Predictions

The integral of Equation 8 over $E_{\rm iso}$ at a given $z$, averaged over
MCMC parameters (Table 1), provides the total, predicted event rate for on-
axis sGRBs as a function of redshift. If we omit the first $f_{b}$ factor from
the integral, we can determine the predicted event rate for all sGRBs on- or
off-axis as a function of redshift. These curves and their uncertainties given
the spread in model parameters are displayed in Figure 5. The curves are
generated assuming the Swift sensitivity but all-sky coverage.

At large distances ($D\gtrsim 1$ Gpc), the satellite sensitivity becomes
important as many events are lost below threshold. Curves for hypothetical
satellites with sensitivities two times and ten times greater than Swift are
displayed. For distances below 1 Gpc and above 200 Mpc, the dominant rate
uncertainty is due to the uncertainty in the total number of BNS mergers as
constrained by LIGO (The LIGO Scientific Collaboration et al., 2021). This
uncertainty is displayed using a dot-dashed curve at the top of Figure 5.

Below 200 Mpc, essentially all sGRBs are detected by Swift according to our
modelling. The number of events detectable all-sky per year at Swift
sensitivities within the aLIGO volume is $0.18^{+0.19}_{-0.08}$ for on-axis
events – in agreement with our rough calculation in Section 1 – and increases
to $1.2^{+1.9}_{-0.6}$ with the inclusion of off-axis events.

###  4.3 Rate Dependence on Beaming, Other Parameters

In Figure 6, we display the predicted rate determined for each parameter draw,
with a separate sub-panel dedicated to each parameter. A quadratic curve is
fit in each sub-panel to highlight trends. Significant trends are observed for
3 parameters: $n$, $c$, and $\alpha$.

When the low-energy powerlaw index for the luminosity function $\alpha$ is
steep, there is naturally an increased rate of faint, nearby events. This
effect can be observed weakly in the $E_{\rm min}$ parameter as well. However,
the impact of that parameter is limited by the fact that $E_{\rm min}$ values
below the observed range of measured $E_{\rm iso}$ values would lead to an
over-large model normalization and be penalized by the exponential term in
Equation 9.

Several of the parameters are correlated with each other as a result of the
flux limit and the beaming constraint Equation 5. The rate-density parameter,
however, is uncorrelated with the other parameters, except for a weak
correlation with $\alpha$ for long delay times ($\tau_{\rm min}>500$ Myr).

We note that our derived luminosity function slopes $\alpha\approx-2$ are
considerably steeper than $\alpha\approx-0.5$ found in Ghirlanda et al.
(2016). While we do not claim strong evidence for a break in the luminosity
function (Section 4.1), and therefore a clean separation between the slopes
$\alpha$ and $\beta$, we do find that the overall luminosity function is
steep.

Figure 7: Left: Mean inverse beaming fraction $\langle f_{b}^{-1}\rangle$ –
proportional to the number of missed events – and confidence intervals (dotted
lines, blue shaded regions) of sGRBs as a function of $z$. Right: The mean
distribution of sGRB jet angles dNobs/d$\theta_{\rm jet}$. We show the
predicted observed (blue, $z>1$ green) distributions as well as those expected
for an ideal instrument (i.e., $z=0$; red dotted). Structured jets (dashed
curves) are discussed in Section 5.

The steeper luminosity function likely arises from our approach to the
modelling of beaming. A large, negative value for $n$ corresponds to a highly-
beamed population where the bulk of sGRBs are pointed away from the observer.
It is not possible to generate the observed sGRB population from compact
object mergers if $n$ is too negative. However, with $f_{b}$ present in the
luminosity function, an interesting effect arises where a steep $\alpha<-1$
permits a larger $n$ (i.e., more GRBs pointed at the observer; see, Equation
5).

The effect of a correlation index $c$ near unity is to make the luminosity
function depend on $E_{\gamma}$. The data prefer a luminosity function
depending on $E_{\gamma}$ rather than $E_{\rm iso}$ ($c>0.5$ at 98%
confidence). In this case, there arises a significant detection bias in favor
of narrowly-beamed events with increasing redshift (see, Figure 7). The
observed cosmological sGRB population (i.e., at high-$z$) appears narrowly
beamed (satisfying the constraint Equation 5), while the low-$z$ population
tends to be more broadly beamed, leading to a higher predicted rate.

Direct constraints on sGRB jet opening angles from afterglow jet-break
measurements are rare. Fong (2014) present a complete tabulation of sGRB
beaming angles measurements and limits available prior to 2014. From 4
measurements of jets with $\theta_{\rm jet}\approx 5$ degrees and 9 upper
limits (as large as 25 degrees), Berger (2014) estimate a mean $\theta_{\rm
jet}>10$ degrees. Our predictions are consistent with these constraints (see,
Figure 7; Right), allowing for large beaming angles but implying that the
cosmological population should appear significantly more tightly beamed than
the general population. We find a most-likely $\theta_{\rm jet}=2.6$ degrees
for $z>1$ as compared to $\theta_{\rm jet}=9$ degrees for all $z$.

An additional important model parameter not displayed in Figure 6 is the
normalization. To reproduce all cosmological sGRBs we must exploit the large
(factor of 10) uncertainty in the BNS rate within 1 Gpc (Sections 1, 3.1). We
find that this factor must be $>5$ (90% confidence), corresponding to a BNS
rate $>10^{3}$ Gpc-3.

##  5 Conclusions

We have demonstrated that a sample of 123 sGRBs from Swift – including 27
events with measured redshifts – can be generated from the source aLIGO merger
population, in turn deriving probabilistic estimates for the sGRB rate
density, the luminosity function, and the beaming distribution. The rate-
density appears to allow for a range of merger delay times ($\tau_{\rm
d,min}<700$ Myr) and suggests contributions from BNS and NSBH progenitor
channels (e.g., Safarzadeh et al., 2019; Sarin et al., 2022), although we have
not attempted to disentangle these contributions. The luminosity function is a
relatively steep ($\alpha\approx-2$), possibly-broken powerlaw operating on
$E_{\rm iso}$ and the beaming fraction $f_{b}$ and leading to a correlation
between the two. As a result, we predict relatively wide sGRB jet angles,
particularly for nearby events.

We predict $0.18^{+0.19}_{-0.08}$ GW/sGRB associations per year, all-sky, for
on-axis events at Swift sensitivities within the $D<200$ Mpc aLIGO volume.
This increases to $1.2^{+1.9}_{-0.6}$ yr-1 with the inclusion of off-axis
events. These very nearby events have a broad, jet-opening angle distribution,
with mean $\theta_{\rm jet}\gtrsim 30$ degrees, while cosmological sGRBs above
the detection limit appear to have narrower jet opening angles as a result of
a correlation between $E_{\rm iso}$ and the beaming factor (e.g., Figure 7).

Figure 8: The joint $E_{\rm iso},f_{b}$ distribution at $z=0$. The correlation
between $E_{\rm iso}$ and $f_{b}$ appears as a broad ridge around $E_{\rm
iso}f_{b}^{c}=10^{47.5}$ (solid blue curve) erg as a consequence of the steep
luminosity function with a minimum energy (Section 4.1). The similar curve
corresponding to a “standard” sGRB $E_{\gamma}$ ($E_{\rm
iso}\propto\theta_{\rm jet}^{2}$) using the small-angle formula is shown for
reference (dashed red curve).

In the context of our sGRB world model (Section 3), it is interesting to
explore the uniqueness of GRB 170817A, associated with GW170817 (Abbott et
al., 2017). The Fermi/GBM fluence – measured over the $T_{\rm 90}$ duration
interval – for GRB 170817A was $2.8\times 10^{-7}$ erg cm-2(10-1000 keV;
Goldstein et al., 2017). The corresponding $E_{\rm iso}=5.35\times 10^{46}$
erg. In our model (see, Figure 8), such a low $E_{\rm iso}$ would be
associated with a spherical, or nearly spherical explosion. However, GRB
170817A was an exceptionally off-axis event, with an X-ray afterglow that
peaked 160 days after the sGRB (Troja et al., 2020). Although the on-axis
energetics remain unconstrained, it is likely that the GRB was viewed at an
angle $\theta_{\rm view}\sim 5\theta_{\rm jet}$ (see, Nakar & Piran, 2021, and
references therein). In our model, a broadly-jetted $\theta_{\rm jet}=15$
degree and $E_{\rm iso}=10^{49}$ erg event is $10^{2}$–$10^{3}$ times more
likely than a $\theta_{\rm jet}=5$ degree event with energies comparable to
those of the brightest, cosmological sGRBs ($E_{\rm iso}\approx
10^{51}$–$10^{52}$ erg; Figure 4).

Although the off-axis nature of GRB 1701817A is clearly an important clue as
to the GW association, it does not suggest that we should increase our on-axis
rate predictions significantly toward our (larger) on- and off-axis rate
predictions (Figure 5). In the flat-top jet model (e.g., Ioka & Nakamura,
2017) or structured jet models (e.g., Nakar & Piran, 2021), the on-axis event
could be very bright ($E_{\rm iso}\gtrsim 10^{52}$ erg), hence highly-unusual
at such a nearby distance of 40 Mpc. Similarly, if $E_{\rm iso}=5.35\times
10^{46}$ erg is instead taken to be typical of the nearly-spherical emission
that may accompany sGRBs, it is only detectable at distances $D\lesssim 60$
Mpc by Swift or Fermi, leading to a predicted rate increase of only 16% in a
200 Mpc aLIGO volume.

A Gaussian structured jet model can increase the observable solid angle
significantly for nearby events (by a factor of 1.4 for $D<200$ Mpc, and 2.2
for $D<40$ Mpc; Figure 5; see, also Howell et al., 2019). Here, instead of
$f_{b}=1-\cos(\theta_{\rm jet})$, we write:

$f_{b}=\int\exp(-(\theta/\theta_{\circ})^{2})\sin(\theta)d\theta.\vspace{0.05in}$
(10)

To capture the same asymptotic behavior for small $\theta_{\rm jet}$ and
$\theta_{\rm jet}=90$ degrees, we set $\theta_{\circ}=\tan(\theta_{\rm jet})$.
Equation 10 does not change the above analysis for cosmological GRBs which
sample only the core of the jet; however the observable angles – for nearby
and narrowly-jetted events where the minimum-detectable $E_{\rm
iso,min}<E_{\rm min}/f_{b}^{c}$ – can be significantly wider than $\theta_{\rm
jet}$. This is because the first factor of $f_{b}$ in Equation 6 is replaced
with $1-\cos(\theta_{\rm view})$, where $\theta_{\rm
view}^{2}=\theta_{\circ}^{2}\ln(L_{\gamma}/L_{\gamma,\rm min})$. This also
tends to favor narrowly-jetted events. For $D<40$ Mpc, the most likely
$\theta_{\rm jet}=29$ degrees (see, Figure 7) and the mean viewing angle is
$\theta_{\rm view}\approx 60$ degrees.

The resulting sGRB flux is near the detection threshold. It is possible that
such events could be identified, likely in a non-decisive manner, after an
optical counterpart is found and the number of search trials is reduced.

Predicted on-axis sGRBs within the aLIGO volume, however, are well-above
(factor of 26 on average) the Swift/BAT threshold. We estimate that the
majority would be detected even using a down-scaled instrument with 80 times
smaller effective area. More important than sensitivity, Swift itself can be
expected to detect very few associations (0.018 yr-1 per 1.3 sr) due to it’s
narrow field of view. Experiments with large fields of view, like Fermi/GBM
(0.09 yr-1 per $2\pi$ sr), are far more likely to be successful. To the extent
that GRB 170817A is in any way typical, and off-axis emission is generally-
detectable, a sensitive and wide-field experiment like Fermi/GBM is optimal
(approaching 0.6 yr-1 per $2\pi$ sr).

Nonetheless, our estimates are fundamentally-limited by a scarcity of
constraints on sGRB beaming. Additional sGRB detections and followup –
including deep followup with Rubin Observatory (Ivezić et al., 2019) to
measure or limit the orphan afterglow (Rhoads, 1999) rate – is needed to
better constrain sGRB jetting and to further probe the relation between sGRBs
and GW events.

Similarly, future GW detections in O4 and O5 with increasing sensitivity
(e.g., Williams et al., 2018), can further constrain the BNS rate beyond 200
Mpc and place more challenging constraints on sGRB beaming. Under the
uncertainty in the fraction of GW events that result in jetting (The LIGO
Scientific Collaboration et al., 2021), reproducing the GRB rate from the GW
rate will be a greater challenge if not all GW events can produce a jet. If
future GW observations constrain the BNS rate to $\lesssim 1$ Gpc-1 (see,
Section 4.3), our starting assumption that all – or even most – sGRBs are due
to compact object mergers is likely incorrect.

## References

  * Abbott et al. (2017a) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017a, The Astrophysical Journal, 848, L13, doi: 10.3847/2041-8213/aa920c
  * Abbott et al. (2017b) —. 2017b, Phys. Rev. Lett., 119, 161101, doi: 10.1103/PhysRevLett.119.161101
  * Abbott et al. (2017) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, ApJ, 848, L12, doi: 10.3847/2041-8213/aa91c9
  * Abbott et al. (2020) —. 2020, ApJ, 892, L3, doi: 10.3847/2041-8213/ab75f5
  * Antonelli, L. A. et al. (2009) Antonelli, L. A., D’Avanzo, P., Perna, R., et al. 2009, A&A, 507, L45, doi: 10.1051/0004-6361/200913062
  * Band (2003) Band, D. L. 2003, ApJ, 588, 945, doi: 10.1086/374242
  * Berger (2014) Berger, E. 2014, Annual Review of Astronomy and Astrophysics, 52, 43, doi: 10.1146/annurev-astro-081913-035926
  * Berger et al. (2013) Berger, E., Zauderer, B. A., Levan, A., et al. 2013, The Astrophysical Journal, 765, 121, doi: 10.1088/0004-637x/765/2/121
  * Bhat et al. (2016) Bhat, P. N., Meegan, C. A., von Kienlin, A., et al. 2016, The Astrophysical Journal Supplement Series, 223, 28, doi: 10.3847/0067-0049/223/2/28
  * Burns et al. (2018) Burns, E., Veres, P., Connaughton, V., et al. 2018, The Astrophysical Journal, 863, L34, doi: 10.3847/2041-8213/aad813
  * Butler et al. (2010) Butler, N. R., Bloom, J. S., & Poznansk, D. 2010, The Astrophysical Journal, 711, 495, doi: 10.1088/0004-637X/711/1/495
  * Butler et al. (2007) Butler, N. R., Kocevski, D., Bloom, J. S., & Curtis, J. L. 2007, The Astrophysical Journal, 671, 656, doi: 10.1086/522492
  * Coward et al. (2012) Coward, D. M., Howell, E. J., Piran, T., et al. 2012, Monthly Notices of the Royal Astronomical Society, 425, 2668, doi: 10.1111/j.1365-2966.2012.21604.x
  * D’Avanzo et al. (2014) D’Avanzo, P., Salvaterra, R., Bernardini, M. G., et al. 2014, MNRAS, 442, 2342, doi: 10.1093/mnras/stu994
  * Fong (2014) Fong, W.-f. 2014, PhD thesis, Harvard Smithsonian Center for Astrophysics
  * Frail et al. (2001) Frail, D. A., Kulkarni, S. R., Sari, R., et al. 2001, ApJ, 562, L55, doi: 10.1086/338119
  * Gehrels (1986) Gehrels, N. 1986, ApJ, 303, 336, doi: 10.1086/164079
  * Gehrels et al. (2005) Gehrels, N., Chincarini, G., Giommi, P., et al. 2005, The Astrophysical Journal, 621, 558, doi: 10.1086/427409
  * Ghirlanda et al. (2016) Ghirlanda, G., Salafia, O. S., Pescalli, A., et al. 2016, A&A, 594, A84, doi: 10.1051/0004-6361/201628993
  * Goldstein et al. (2017) Goldstein, A., Veres, P., Burns, E., et al. 2017, ApJ, 848, L14, doi: 10.3847/2041-8213/aa8f41
  * Gruber et al. (2014) Gruber, D., Goldstein, A., von Ahlefeld, V. W., et al. 2014, The Astrophysical Journal Supplement Series, 211, 12, doi: 10.1088/0067-0049/211/1/12
  * Guetta & Piran (2005) Guetta, D., & Piran, T. 2005, A&A, 435, 421, doi: 10.1051/0004-6361:20041702
  * Guetta & Piran (2006) —. 2006, A&A, 453, 823, doi: 10.1051/0004-6361:20054498
  * Guetta & Stella (2008) Guetta, D., & Stella, L. 2008, Astronomy & Astrophysics, 498, 329, doi: 10.1051/0004-6361:200810493
  * Howell et al. (2019) Howell, E. J., Ackley, K., Rowlinson, A., & Coward, D. 2019, Monthly Notices of the Royal Astronomical Society, 485, 1435, doi: 10.1093/mnras/stz455
  * Ioka & Nakamura (2017) Ioka, K., & Nakamura, T. 2017, Progress of Theoretical and Experimental Physics, 2018, doi: 10.1093/ptep/pty036
  * Ivezić et al. (2019) Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111, doi: 10.3847/1538-4357/ab042c
  * Jin et al. (2018) Jin, Z.-P., Li, X., Wang, H., et al. 2018, ApJ, 857, 128, doi: 10.3847/1538-4357/aab76d
  * Kouveliotou et al. (1993) Kouveliotou, C., Meegan, C. A., Fishman, G. J., et al. 1993, ApJ, 413, L101, doi: 10.1086/186969
  * Levesque et al. (2010) Levesque, E. M., Bloom, J. S., Butler, N. R., et al. 2010, Monthly Notices of the Royal Astronomical Society, 401, 963, doi: 10.1111/j.1365-2966.2009.15733.x
  * Lien et al. (2016) Lien, A., Sakamoto, T., Barthelmy, S. D., et al. 2016, The Astrophysical Journal, 829, 7, doi: 10.3847/0004-637X/829/1/7
  * Liu & Yu (2019) Liu, H.-Y., & Yu, Y.-W. 2019, Research in Astronomy and Astrophysics, 19, 118, doi: 10.1088/1674-4527/19/8/118
  * Longair (1966) Longair, M. S. 1966, MNRAS, 133, 421, doi: 10.1093/mnras/133.4.421
  * Nakar et al. (2006) Nakar, E., Gal-Yam, A., & Fox, D. B. 2006, ApJ, 650, 281, doi: 10.1086/505855
  * Nakar & Piran (2021) Nakar, E., & Piran, T. 2021, ApJ, 909, 114, doi: 10.3847/1538-4357/abd6cd
  * Petrosian & Lee (1996) Petrosian, V., & Lee, T. T. 1996, in American Institute of Physics Conference Series, Vol. 366, High Velocity Neutron Stars, ed. R. E. Rothschild & R. E. Lingenfelter, 170–174
  * Rhoads (1999) Rhoads, J. E. 1999, ApJ, 525, 737, doi: 10.1086/307907
  * Safarzadeh et al. (2019) Safarzadeh, M., Berger, E., Ng, K. K. Y., et al. 2019, The Astrophysical Journal, 878, L13, doi: 10.3847/2041-8213/ab22be
  * Salvaterra et al. (2008) Salvaterra, R., Cerutti, A., Chincarini, G., et al. 2008, MNRAS, 388, L6, doi: 10.1111/j.1745-3933.2008.00488.x
  * Salvatier et al. (2016) Salvatier, J., Wiecki, T. V., & Fonnesbeck, C. 2016, PeerJ Computer Science, 2, e55
  * Sarin et al. (2022) Sarin, N., Lasky, P. D., Vivanco, F. H., et al. 2022, Phys. Rev. D, 105, 083004, doi: 10.1103/PhysRevD.105.083004
  * Selsing, J. et al. (2018) Selsing, J., Krühler, T., Malesani, D., et al. 2018, A&A, 616, A48, doi: 10.1051/0004-6361/201731475
  * Siellez et al. (2013) Siellez, K., Boër, M., & Gendre, B. 2013, Monthly Notices of the Royal Astronomical Society, 437, 649, doi: 10.1093/mnras/stt1915
  * The LIGO Scientific Collaboration et al. (2021) The LIGO Scientific Collaboration, the Virgo Collaboration, the KAGRA Collaboration, et al. 2021, arXiv e-prints, arXiv:2111.03634.  https://arxiv.org/abs/2111.03634
  * Troja et al. (2018) Troja, E., Ryan, G., Piro, L., et al. 2018, Nature Communications, 9, 4089, doi: 10.1038/s41467-018-06558-7
  * Troja et al. (2020) Troja, E., van Eerten, H., Zhang, B., et al. 2020, Monthly Notices of the Royal Astronomical Society, 498, 5643, doi: 10.1093/mnras/staa2626
  * Virgili et al. (2011) Virgili, F. J., Zhang, B., O'Brien, P., & Troja, E. 2011, The Astrophysical Journal, 727, 109, doi: 10.1088/0004-637x/727/2/109
  * von Kienlin et al. (2014) von Kienlin, A., Meegan, C. A., Paciesas, W. S., et al. 2014, The Astrophysical Journal Supplement Series, 211, 13, doi: 10.1088/0067-0049/211/1/13
  * von Kienlin et al. (2020) —. 2020, The Astrophysical Journal, 893, 46, doi: 10.3847/1538-4357/ab7a18
  * Williams et al. (2018) Williams, D., Clark, J. A., Williamson, A. R., & Heng, I. S. 2018, The Astrophysical Journal, 858, 79, doi: 10.3847/1538-4357/aab847


#  A data processing system for balloon-borne telescopes

Valentina Scotti , Giuseppe Osteria, and Francesco Perfetto Manuscript
received October 30, 2020. This work was partially supported by the Italian
Space Agency through the ASI INFN agreement n. 2017-8-H.0.V. Scotti is with
the Department of Physics of University of Naples Federico II and INFN Naples,
Italy (email: scottiv@na.infn.it)G. Osteria is with INFN Naples, ItalyF.
Perfetto is with INFN Naples, Italy

###### Abstract

The JEM-EUSO Collaboration aims at studying Ultra High Energy Cosmic Rays
(UHECR) from space. To reach this goal, a series of pathfinder missions has
been developed to prove the observation principle and to raise the
technological readiness level of the instrument. Among these, the EUSO-SPB2
(Extreme Universe Space Observatory on a Super Pressure Balloon, mission two)
foresees the launch of two telescopes on an ultra-long duration balloon. One
is a fluorescence telescope designed to detect UHECR via the UV fluorescence
emission of the showers in the atmosphere. The other one measures direct
Cherenkov light emission from lower energy cosmic rays and other optical
backgrounds for cosmogenic tau neutrino detection.

In this paper, we describe the data processing system which has been designed
to perform data management and instrument control for the two telescopes. It
is a complex which controls front-end electronics, tags events with arrival
time and payload position through a GPS system, provides signals for time
synchronization of the event and measures live and dead time of the telescope.
In addition, the data processing system manages mass memory for data storage,
performs housekeeping monitor, and controls power on and power off sequences.

The target flight duration for the NASA super pressure program is 100 days,
consequently, the requirements on the electronics and the data handling are
quite severe. The system operates at high altitude in unpressurised
environment, which introduces a technological challenge for heat dissipation.

###### Index Terms:

Balloon experiments, Front-end electronics for detector readout, On-board data
handling, On-board space electronics, Trigger concepts and systems (hardware
and software)

##  I Introduction

The EUSO-SPB2 is a precursor mission for future space observatories for multi-
messenger astrophysics.

The innovative approach of EUSO-SPB2 is to use the Earth as a giant particle
detector, studying fluorescence and Cherenkov light from the edge of the
space.

The main goal of this mission is to detect fluorescence light generated in the
Earth’s atmosphere by Ultra High Energy Cosmic Rays (UHECR) from above, to
confirm the expectations from ground observations [1, 2].

In addition, this mission will explore new areas such as detecting Cherenkov
light emission from air showers and measuring the background of up-going tau
decays from cosmogenic neutrinos [3].

This mission is part of the JEM-EUSO science program. JEM-EUSO is a
collaboration of 350 researchers from 16 countries, whose scientific case has
been positively evaluated by several space agencies with funding ongoing in
all the participating countries.

The aim of JEM-EUSO is to observe from space UV light in the terrestrial
atmosphere to study Extreme Energy Cosmic Rays (EECR): above $5\times
10^{19}$eV. The main challenge for the study of EECR is the extremely low flux
of cosmic rays at the highest energies.

The scientific program of JEM-EUSO includes a series of pathfinder experiments
using fluorescence detectors to make a proof-of-principle of the EECR
observation from space and to raise the technological level of the
instrumentation to be employed in a space mission.

Among these, the EUSO-SPB2 is build upon the EUSO-SPB1 [4] experience to pave
the way towards the Probe of Extreme Energy Multi-Messenger Astrophysics
(POEMMA) space mission [3]. POEMMA is designed to discover the nature and
origin of EECR and to discover astrophysical neutrino emission.

The launch of the EUSO-SPB2 for an ultra-long duration balloon flight is
foreseen in 2022 from Wanaka in New Zealand. The target flight duration for
the NASA super pressure balloon program is 100 days at a floating altitude of
33 km.

A description of the mission is given in [5], in this paper, we present the
concept of the Data Processor system designed to manage the main instruments
and all the devices on board while respecting the severe constraints imposed
by a balloon-borne experiment.

##  II The EUSO-SPB2 telescopes

EUSO-SPB2 will host onboard two independent telescopes, a fluorescence
detector (FT) and a Cherenkov detector (CT), in order to be a scientific and
technical sub-orbital altitude precursor for the POEMMA. A complete
description of the two instruments can be found respectively in [6] and [7],
here we will briefly report only on the details relevant to the functioning of
the Data Processing system. Fig. 1 shows a drawing of the gondola structure.

Figure 1: A pictorial view of the gondola hosting the two telescopes and all
the ancillary devices of the EUSO-SPB2 mission.

In table I the main characteristics of the two telescopes are summarized.

TABLE I: Main parameters of the EUSO-SPB2 telescopes. The two telescopes share the same mechanical structure, mirror, and a 1-meter wide PMMA corrector plate. Fluorescence |  | Cherenkov
---|---|---
| Energy threshold (eV) |
| Sensor type |
300$\div$420 nm, BG3 filter | Wavelenght sensitivity | 400$\div$800 nm, no filter
| Time unit |
| Pointing (zenith angle) |
| FoV (instrumented) |
| Number of pixels |
| Optics (modified Schmidt) |
| Payload mass (kg) |

The fluorescence detector is build on the experience of EUSO-SPB1 to detect
ultra high energy cosmic rays via the UV fluorescence emission of the
extensive air showers generated by UHECR in the atmosphere. The main
difference with respect to EUSO-SPB1 is that in EUSO-SPB2 the light is focused
on the camera through a Schmidt optical system, instead of Fresnel lenses. In
addition, the focal surface of the FT will be composed of 3 Photo Detector
Modules (PDM), each composed of 36 MAPMT Hamamatsu R11265, for a total of 6912
pixels, compared to 2304 of EUSO-SPB1. The Field of View of each PDM is
$12^{\circ}\times 12^{\circ}$. The signal from MaPMT are with are readout by a
specific ASIC followed by multiple trigger levels to filter out noise and
identify events of interest. This system has single photoelectron counting
capability and a time unit of $1\mu s$. The FT points to nadir and the
expected trigger rate is around 0.2/hour for a trigger threshold of
$10^{18}$eV.

The CT points toward the limb to measure direct Cherenkov light emission from
lower energy cosmic rays and the optical backgrounds for cosmogenic tau
neutrino detection in the UV/VIS range. For the CT the mirror segments will be
aligned in a bifocal configuration which allows to distinguish between a
charged particle hitting camera which makes 1 spot, and a light pulse from
far-field outside the telescope makes 2 spots. The CT focal surface will be
composed by 32 Hamamatsu S14521 (4 × 4 pixel array) SiPMT, for a total of 512
pixels. The readout is based on a 100 MS/s ASIC for General Electronics for
TPC‘s (AGET) switch capacitor ring sampler.

Each telescope is completed by several ancillary devices:

  * •

Aperture shutter to prevent sun-melt and avoid dust entering the camera. The
shutter is light tight and makes the inside of the telescope a dark box.

  * •

Health LED Systems to check camera health and response.

  * •

Airglow MONitor (AMON) for an independent measure of airglow and night sky
radiation backgrounds.

  * •

Light sensor, called EMON, located inside the telescopes to verify darkness
condition during operation and testing.

  * •

A dual-band InfraRed camera to monitor the presence of clouds and to measure
cloud-top temperatures only for the fluorescence telescope.

##  III The data processing system

The focal plane detectors and the front-end electronics of the two telescopes
are different since the Cherenkov and fluorescence light have different
wavelengths and time development. However, from the point of view of the
control, configuration, monitoring, data management and working mode, the two
instruments are very similar. For this reason, the architecture of the Data
Processing system designed for the two telescopes is an evolution of the one
developed for the fluorescence telescope of the EUSO-SPB1 [8].

The block diagram of the DP of the fluorescence telescope and its connections
with the rest of the instrument is shown in Fig. 2.

Figure 2: Block diagram of the Data Processing system nd its connections with
the rest of the instrument for the fluorescence telescope of the EUSO-SPB2.
The one for the Cherenkov detector is quite similar.

The one for the Cherenkov detector is quite similar. In the following
paragraph we will describe the Data Processing system of the FT.

###  III-A Requirements

The EUSO-SPB2 telescopes will collect data during moonless nights to obtain
the best contrast for the signatures of scientific interest. The data
processing system (DP) links each telescope with the Gondola system. The DP is
a complex system which includes most of the digital electronics of the
instrument which allows to control, configure, monitor and operate each
telescope during the commissioning phase, the test campaigns and the flight.

Besides being the interface with the flight computer, the data processing
system performs the following tasks:

  * •

Main interface with Flight Computer (SIP) telemetry system

  * •

Data selection/compression and transmission to Flight Computer (SIP)

  * •

Power ON/OFF the whole instrument

  * •

Define Telescope operation mode (Day, Night, D-N-D transitions)

  * •

Configure the Front End electronics

  * •

Start/Stop of the data acquisition and calibration procedures

  * •

Tag events with GPS time and GPS position

  * •

Synchronization of the data acquisition

  * •

Manage trigger signals (L1, external, GPS, by CPU command etc)

  * •

Manages mass memory for data storage

  * •

Monitor/Control/DAQ of some Ancillary Devices

  * •

Monitor voltages, current and temperatures (LVPSs, boards, FPGAs)

The main challenges in designing and building the system are due to the
physical characteristics of the signal and to the harsh environment in which
the DP has to operate.

At the Peak energy sensitivity for the FT (around 6 EeV) the expected event
rate is around 1 event every 13 hours. For an 80 day flight with 15% duty
cycle the total number of events is around 22 events, but the trigger rate of
the FT is estimated to be between 5 and 10 Hz, so it’s necessary to
distinguish a few cosmic rays events among the background.

Furthermore, since the target flight duration for the NASA super pressure
program is 100 days, the requirements on the electronics and the data handling
are quite severe. Finally, The system operates at high altitude in
unpressurised environment, which introduces a technological challenge for heat
dissipation.

In addition to power and mass budget restrictions on the balloon payload,
there is also a limited telemetry budget of 1GB per day. For this reason,
there is the need to prioritize data for downloading on board, at the same
time with the want to only record and transfer high quality events. Everything
must be done on board with minimal intervention from the ground.

##  IV The subsystems

The block diagram (Fig. 2) shows that the data processing system is composed
of several subsystems. The whole system is capable to acquire 7000 channels
without exceeding mass and power budget. ———-

Different tasks are performed by different subsystems. All the subsystems,
together with the low voltage power supplies modules, are hosted in a DP box
(customized Eurocard chassis) equipped with cooling plate to dissipate heat
(Fig. 3).

Figure 3: The CAD of the DP box, designed to host all the Data Processing
subsystems and equipped with cooling plate.

  * •

CPU with Hot and Cold redundancy PCI/104 single board computer Core i7 3517UE
+ 5 SATA disks for data storage

  * •

HK board board based on STM32 µController and TIBO EM2000

  * •

Clock board [9] based on Xilinx XC7Z SoC: manages trigger and data
synchronization and also acts as an interface with the GPS receivers

  * •

GPS receivers Trimble BX992

  * •

Ethernet switches

  * •

Solid State Power Controller

We choose to use mostly COTS devices because reliability in a hostile
environment is one of the requirements. The electronic components are all
selected to operate in extended temperature range.

For the most fundamental component, the CPU, we have hot and cold redundancy,
hence the choice to use CAN and Ethernet protocol which guarantee us in case
one CPU fails.

##  V The software

The DP is the interface between the detector and the telemetry blocks and the
end-users.Instrument operations and monitoring will be shared between
collaborators in the US, Japan and Europe to provide, during local daytime,
control and monitoring of the telescopes 24/7, similarly to the on-ground
control that was done for EUSO-SPB1 [10].

The CPU software can be divided into two categories:

  * •

the control software

  * •

the data handling software

The control software has to be flexible to be used in every phase of the
commissioning and flight. It allows to control and monitor the detector, and
to store all the data (science, housekeeping and log). The design enables full
user control via telemetry. In fact, during the flight the status of the
apparatus must be continuously checked to distinguish condition to start or
stop the measurement.

The data handling software manages science data acquisition from the CLK and
Zynq boards. The acquisition software initializes and configures the
subsystems, monitors the connections with the CPU, and verifies the behavior
of the subsystems. It can perform several types of acquisitions, such as
recording of external triggered events or internal triggered events. Finally,
it manages the prioritization folder structure for down-link to cope with the
maximum amount of telemetry available each day (1 GB).

Science Stack data are available in ”quasi” real time on a WEB page of
Columbia Scientific Balloon Facilities and are used to monitor and to control
the status of the telescope. In particular, before starting the power-on or
power-off sequence, it is necessary to monitor:

  * •

Status of the power system

  * •

Temperatures of critical point of the telescope (Batteries, mirror, focal
surface, electronics, etc.)

  * •

Light level of the EMON (Day/Night transition, Night/Day transition)

  * •

Status of the aperture shutter (Open/Close)

##  VI Schedule

The launch of the EUSO-SPB2 is planned for 2022, in order to reach this goal,
the EUSO-SPB2 Collaboration need to hit several milestones.

CPUs, data storage, GPS receivers, Ethernet switches, Solid State Power
Controller are under test in different laboratories. The first prototype of
the Clock board has been manufactured and tested, the firmware is under
development. The design of the HK board has been completed and the first
prototype will be available before the end of 2020. The design of the DP
housing and cooling mechanics are in progress and will be soon finalized and
tested. The Collaboration is currently looking for facilities for thermal-
vacuum tests of the whole DP. Extensive laboratory testing is planned to
characterize and calibrated the telescopes prior to payload integration.

At the end of the assembly, and after thermal-vacuum chamber test, both the
telescopes will be transported to the Utah desert for a field test campaign.
The aim of this campaign is to use lasers and other light sources to have an
absolute calibration of the detector response. The field tests will include
observation time to record Fluorescence and Cherenkov light due to cosmic
rays. The results of the tests will be used also to fine-tune the trigger
algorithm. .

##  VII Conclusions

EUSO-SPB2 is a 2-telescope instrument to study UHECR, Cherenkov, and tau
neutrino backgrounds. The launch planned for 2022. The Data Processing system
of the telescopes have been designed and the first prototypes are under test.
The Collaboration is getting ready for thermal-vacuum test and instrument
integration.

## Acknowledgment

We would like to thank the staff of Electronics and Detectors Service (SER)
and Design and Mechanics Service of INFN Napoli for their strong and undivided
support all along this project.

## References

  * [1] C. Berat, for the Pierre Auger Collaboration2, The Pierre Auger Observatory status and latest results, EPJ Web of Conferences 206, 07001 (2019), https://doi.org/10.1051/epjconf/201920607001
  * [2] The IceCube Collaboration, The Pierre Auger Collaboration, The Telescope Array Collaboration, The ANTARES Collaboration, Search for correlations of high-energy neutrinos and ultra-high energy cosmic rays, 36th ICRC Madison Wisconsin, 2019.
  * [3] A. Olinto, et al., POEMMA (Probe of Extreme Multi-Messenger Astrophysics) design, arXiv:1907.06217 [astro-ph.HE]
  * [4] J. Eser, on behalf of the JEM-EUSO Collaboration, Results of the EUSO-SPB1 flight, POS ICRC2019, doi = 10.22323/1.358.0247
  * [5] V. Scotti, G. Osteria, The EUSO-SPB2 mission, Nuclear Instruments and Methods in Physics Research Section A, Volume 958, 2020, https://doi.org/10.1016/j.nima.2019.05.005.
  * [6] Abdellaoui, G., Abe, S., et al., Ultra-violet imaging of the night-time earth by EUSO-Balloon towards space-based ultra-high energy cosmic ray observations, Astroparticle Physics, 111, pp. 54-71. DOI: 10.1016/j.astropartphys.2018.10.008
  * [7] A. N. Otte, et al., Development of a Cherenkov Telescope for the Detection of Ultrahigh Energy Neutrinos with EUSO-SPB2 and POEMMA, arXiv:1907.08728 [astro-ph.IM]
  * [8] V. Scotti, J. Bayer, et al., The Data Processor system of EUSO-SPB1, Nuclear Instruments and Methods in Physics Research Section A, Volume 916, 2019, Pages 94-101, https://doi.org/10.1016/j.nima.2018.10.207.
  * [9] V. Scotti, G. Osteria, on behalf of the JEM-EUSO Collaboration, The JEM-EUSO time synchronization system, Nuclear Instruments and Methods in Physics Research Section A:, Volume 718, 2013, Pages 248-250, https://doi.org/10.1016/j.nima.2012.10.121.
  * [10] Fornaro, C, Cafagna, et al., The onboard software of the EUSO‐SPB pathfinder experiment. Softw: Pract Exper. 2019; 49: 524– 539. https://doi.org/10.1002/spe.2655


# Nonlinear Radiation Damping of Nuclear Spin Waves and Magnetoelastic Waves
in Antiferromagnets

Alexander V. Andrienko<EMAIL_ADDRESS>National Research Centre “Kurchatov
Institute”, 123182 Moscow, Russian Federation  Vladimir L. Safonov
<EMAIL_ADDRESS>Mag and Bio Dynamics, Inc., Arlington, TX 76001;
Tarrant County College, Fort Worth, TX 76119.

###### Abstract

Parallel pumping of nuclear spin waves in antiferromagnetic CsMnF3 at liquid
helium temperatures and magnetoelastic waves in antiferromagnetic FeBO3 at
liquid nitrogen temperature in a helical resonator was studied. It was found
that the absorbed microwave power is approximately equal to the irradiated
power from the sample and that the main restriction mechanism of absortption
in both cases is defined by the nonlinear radiation damping predicted about
two decades ago. We believe that the nonlinear radiation damping is a common
feature of parallel pumping technique of all normal magnetic excitations and
it can be detected by purposeful experiments.

††preprint: This line only printed with preprint option

76.50.+g, 76.90.+d

##  I introduction

Microwave parametric resonance of normal magnetic oscillations, such as
electronic spin waves, nuclear spin waves and magnetoelastic waves is a
powerful tool to study linear and nonlinear properties of magnetoordered
systems (ferromagnets, antiferromagnets and ferrites) suhl ; schloemann ; ZLS
; L ; GM ; VS  . Parallel pumping of magnetic excitations, in which the
microwave magnetic polarization is parallel to the external magnetic field, is
one of the most convenient and popular methods of parametric resonance. A
microwave magnetic field $\mathbf{h}(t)$ enhanced by a microwave resonator is
applied to the sample parallel to its equilibrium magnetization, which is
parallel to the steady external magnetic field $\mathbf{H}$. The alternating
magnetic field excites the parametric resonance of the form
$\omega_{p}=\omega_{\mathbf{k}}+\omega_{-\mathbf{k}}$, where $\omega_{p}$ is
the pumping field frequency and $\omega_{\mathbf{k}}=\omega_{-\mathbf{k}}$ are
the half-pump frequencies of excited in the sample parametric pair of waves
with oppositly oriented wave vectors ${\mathbf{k}}$ and $-\mathbf{k}$.

The excited waves amplitudes grow exponentially when the microwave field
amplitude $h$ exceeds the parametric resonance threshold $h_{c}$. According to
S theory ZLS ; L , this growth is restricted by the nonlinearities of the
magnetic system, which are exhibited by a) phase mismatching of the forced
magnetic oscillations with the microwave field and b) by the positive
nonlinear magnetic relaxation due to nonlinearities of magnetic system. In
this theoretical picture the resonator cavity is considered to be just as an
ancillary system that enhances microwave field amplitude on a sample. No other
effect associated with the microwave resonator cavity is assumed in this small
sample approximation approach.

Actually, the process of parallel pumping includes two steps. First, the
external microwave source excites the same frequency $\omega_{p}$ microwave
magnetic oscillation of the resonator cavity. Second, this magnetic
oscillation is absorbed by the parametic pair ($\omega_{\mathbf{k}}$ and
$\omega_{-\mathbf{k}}$) of magnetic excitations of the sample. In principle,
one can expect a backward radiation of the parametric pairs and a bunch of
associated effects in the system of two interacting in resonance oscillations.
However, in the simple picture of small sample approximation this backward
radiation is assumed to be negligibly small compared to absorption; only an
energy flow from the microwave pump to the sample occurs.

Notwithstading that examples of non-trivial role of microwave resonator to the
process of parallel pumping of magnetic excitations have already been
discussed in Refs.S jmmm ; SY ; ASirrad ; AShumps ; ASY , these facts did not
attract much attention and the small sample approximation approach is still in
wide use for the description of parallel pumping of magnetic oscillations. The
main focus of the present paper is to demonstrate that theoretically predicted
two decades ago S jmmm ; SY  nonlinear radiation damping, effect due to
backward irradiation of parametric pairs to the resonator, is a common and
dominant feature in the process of parallel pumping of magnetic excitations.

In this paper we studied parallel pumping in a helical resonator of a) nuclear
spin waves in an antiferromagnetic CsMnF3 at liquid helium temperatures and b)
magnetoelastic waves in an antiferromagnetic FeBO3 at liquid nitrogen
temperatures.

The concept of nuclear spin waves was introduced by de Gennes et al NSW . The
nuclear spin wave denotes the magnetic excitation of mixed electronic and
nuclear spin oscillations that is situated in the nuclear magnetic resonance
frequency range. The most remarkable property of these excitations is that, at
liquid helium temperatures, they exhibit the coupled oscillations of two
completely different in their magnetic properties subsystems. The electronic
spins are ordered while the state of the nuclear spins is paramagnetic; the
polarization is no more than several percent. As a result of mixing of these
two subsystems by hyperfine interaction, the frequency of the electronic spin
waves increases, and the nuclear magnetic resonance frequency $\omega_{n,0}$
decreases and becomes noticeably lower than the Larmor precession frequency
$\omega_{n}$ of nuclear spins. In other words, there arises the so-called
dynamic nuclear magnetic resonance pulling and the band of nuclear spin waves
$\omega_{n,k}$:

$\omega_{n,k}=\omega_{n}\left[1-\left(\frac{\gamma
H_{\Delta,hf}}{\omega_{e,k}}\right)^{2}\right]^{1/2},$ (1)

where $\omega_{e,k}=\gamma[H(H+H_{D})+H_{\Delta,hf}^{2}+(\alpha k)^{2}]^{1/2}$
is the frequency of electronic spin wave, $H_{D}$ is the Dzyaloshinskii field,
$H_{\Delta,hf}^{2}\propto 1/T$ is the gap due to hyperfine interaction,
$\alpha$ is the exchange constant and $\gamma$ is the gyromagnetic ratio. The
detailed review of nuclear spin wave properties in weakly anisotropic
antiferromagnets is given in Ref.AOSY .

Magnetoelastic waves describe normal modes of linearly coupled elastic waves
and electronic spin waves in magnetoordered crystals. So far as the
magnetoelastic waves contain both elastic and magnetic components, they can be
excited both by elastic vibrations and by alternating magnetic field. One of
the most interesting objects to study magnetoelastic waves is the high Néel
temperature antiferromagnet FeBO3 ($T_{N}=348$ K). Parallel pumping of
magnetoelastic waves in this crystal for the first time was observed in Ref.AP
. The spectrum of magnetoelastic waves in iron borate can be written as ASY :

$\omega_{me,k}=c_{e}k\left[1-\left(\frac{\gamma
H_{\Delta,ef}}{\omega_{ek}}\right)^{2}\right]^{1/2},$ (2)

where $c_{e}$ is the sound velocity, $H_{\Delta,ef}$ describes an efficiency
of linear interaction between spin and elastic subsystems,
$\omega_{e,k}=\gamma[H(H+H_{D})+H_{\Delta,me}^{2}+(\alpha k)^{2}]^{1/2}$ is
the frequency of electronic spin wave and $H_{\Delta,me}$ is the field which
corresponds to magnetoelastic gap.

We show that beyond the small sample approximation the resonator oscillation
dynamics plays an extremely important role in the process of parametric
resonance of nuclear spin waves and magnetoelastic waves and gives the
dominant mechanism of parallel pumping restriction in both cases by the
nonlinear radiation damping.

##  II experiment

The experimental absorbing cell is shown in Fig.1. The sample is placed in an
open helical resonator which is a half-wavelength dipole excited by the pulsed
microwave pumping field $\mathbf{h}(t)$. The inner diameter of the helix
equals 0.5 cm and the diameter of the copper wire is 0.5 mm. To a first
approximation the wire length needed to make the helix is $\simeq\lambda/2$
which is about 15 cm for 1 GHz. The effective volume of this resonator is
estimated as $\sim 200$ ${\rm mm}^{3}$. The effect of microwave absorption is
detected by the receiving antenna. This absorbing cell to study parallel
pumping of nuclear spin waves and magnetoelastic waves was used at different
temperature conditions.

Figure 1: Schematic diagramm of the experimental absorbing cell.

Parametric pairs of nuclear spin waves were excited by a pulsed ($300-2000$
$\mu{\rm s}$) parallel microwave pump with repeating frequency $10-100$ Hz in
the helical resonator with the quality factor $Q\sim 300-500$ over a wide
range of frequencies $\omega_{p}=600-1200$ MHz. The measurements were made on
single-crystal sample $v_{s}=3\times 3\times 5$ ${\rm mm}^{3}$ of the easy-
plane antiferromagnet CsMnF3 ($T_{N}=53.5$ K) at liquid helium temperatures
$T=1.9-4.2$ K and magnetic fields $H=500-2000$ Oe. The ratio of the sample
volume to the volume of resonator was $v_{s}/v_{R}\sim 0.2$. The relaxation
rate of parametically excited spin waves estimated by the threshold amplituede
was $\eta_{\mathbf{k}}/2\pi\sim 6-20$ kHz with the accuracy of 25 %.

A typical form of the microwave pump pulse passed thorough the resonator is
shown in the left side of Fig.2. There is a microwave absortpion by the
parametric pairs (upper part of the pulse) which is demonstrated by the
decrease of the pump pulse. At the end of the pulse one can see a general
phenomenon, a non-uniform time dependence with a peak of the microwave
radiation after the pump pulse, a typical oscillation of the microwave pulse
which appears above the threshold of parametric resonance. We could observe
this non-trivial radiation at $P/P_{c}-1\gg 1$. The peak demonstrates a
beating of magnetic oscillation of the resonator cavity mode with the
parametic pair ASirrad . In this case the lineshape of the cavity-sample
system becomes splitted into two humps AShumps  which is a direct indication
that the small sample approximation is not valid any more. Experimentally we
observed one peak if the pump frequency was equal to the frequency of the
resonator $\omega_{p}=\omega_{R}$ and up to three beating peaks if
$\omega_{p}\neq\omega_{R}$ . It should be noted that below the threshold of
parametric resonance the microwave radiation after the pump pulse demonstrates
just an exponential decrease (see, curve 2 in Fig.2), which corresponds to
unloaded resonator cavity irradiation.

Figure 2: LEFT: A typical form of the microwave pumping pulse passed through
the helical resonator. One can see a microwave absorption by the sample (upper
part) and a non-uniform radiation effect after the end of mirowave pumping.
RIGHT: Curve 1 demonstrates a non-monotonic radiation power signal from the
sample after the pump pulse was turned off. The pumping power $P\approx 2000$
$P_{c}$. Curve 2 demonstrates the case when $P<P_{c}$, when just an
exponentially decreasing radiation from the resonator cavity is observed. The
experimental parameters are: $T=2.08$ K, $\omega_{p}/2\pi=1094$ MHz and
$H=1840$ Oe.  Figure 3: Irradiation power of magnetoelastic waves versus time
after the end of the microwave pump pulse at two overcriticalities:
$P/P_{c}=13.2$ and$P/P_{c}=52.3$ at $T=77$ K, $\omega_{p}/2\pi=1109.7$ MHz and
$H=231$ Oe. Solid lines describe theoretical fit (see the text).

Parametric pairs of magnetoelastic waves were excited in the $v_{s}\simeq 20$
${\rm mm}^{3}$ sample of the “easy-plane” antiferromagnet FeBO3 by the pulsed
microwave field of the frequency $\omega_{p}/2\pi=900-1200$ MHz at magnetic
fields $H=30-500$ Oe at liquid nitrogen temperature $T=77$ K. The ratio of the
sample volume to the volume of resonator was $v_{s}/v_{R}\sim 0.1$. We
observed similar effects of the non-uniform radiation from the cavity-sample
system after the end of mirowave pump pulse as in the case of nuclear spin
waves. Typical experimental data of irradiation are shown in Fig.3.

We found very important feature of experiments with irradiation: the radiation
power is approximately equal to the absorption power. Thus, the stationary
state of parametric pairs is defined by the radiation from the sample through
the “sample-resonator” nonlinear interaction .

##  III discussion

Let us consider the monotonically decreasing time dependence of radiated power
behind the beating peak. The decrease of the parametric pairs number
$N_{k}(t)$ is described by the equation $dN_{k}=-2\eta(N_{k})dt$, where
$\eta(N_{k})=\eta_{k}+\eta_{nl}N_{k}$ is the relaxation rate, $\eta_{k}$ is
the linear and $\eta_{nl}N_{k}$ is the nonlinear parts, respectively.
Integrating of this equation, one obtains

$N_{k}(t)=\frac{\eta_{k}/\eta_{nl}}{u\exp[2\eta_{k}(t-t_{0})]-1},$ (3)

where $u=1+\eta_{k}/\eta_{nl}N_{k}(t_{0})$, $t_{0}$ is the starting time
($t\geq t_{0}$).

If we assume that the nonlinear part of damping is entirely defined by
nonlinear radiation damping, then the radiated power $P_{rad}(t)$ can be
expressed as

$P_{rad}(t)=-\hbar\omega_{p}\frac{dN_{k}}{dt}\frac{\eta_{nl}N_{k}(t)}{\eta_{k}+\eta_{nl}N_{k}(t)}=\hbar\omega_{p}\frac{2\eta^{2}_{k}/\eta_{nl}}{\\{u\exp[2\eta_{k}(t-t_{0})]-1\\}^{2}}.$
(4)

Figure 4: Radiation power (dots) from the parametrically pumped nuclear spin
waves versus time in CsMnF3 at $T=2.08$ K, $\omega_{p}/2\pi=1094$ MHz and
$H=1840$ Oe. Curve 1 schematically demonstrate the radiation power slope in
the case of linear damping. Curve 2 demonstrate the theoretical fit of formula
(4) (see, the text).

###  III.1 Nuclear Spin Waves

A typical time slope for radiation power is shown in Fig. 4. Mean-square fit
using formula (4) with $t_{0}=0.8$ $\mu\rm{s}$ gives: $\hbar\omega_{p}\cdot
2\eta^{2}_{k}/\eta_{nl}=1.8\cdot 10^{-4}$ W and
$\eta_{k}/\eta_{nl}N_{k}(0.8\,\mu\rm{s})=9.05\cdot 10^{-2}$. The linear
relaxation rate calculated from the threshold of parallel pumping is
$\eta_{k}=4.46\cdot 10^{4}$ s-1. Thus we obtain $\eta_{nl}=1.6\cdot 10^{-11}$
s-1 and $\eta_{nl}N_{k}(0.8\,\mu\rm{s})=4.93\cdot 10^{5}$ s-1 which is one
order greater than the linear relaxation rate $\eta_{k}$. The number of
parametric pairs at $t_{0}=0.8$ $\mu\rm{s}$ is equal to
$N_{k}(0.8\,\mu\rm{s})\simeq 3.1\cdot 10^{16}$. This esimate for the number of
parametric pairs is in agreement with the estimate obtained in Ref.AndriN
from the susceptibility in the overthreshold region.

Note that the obtained result is stable to the variation of $\eta_{k}$. For
example, if we take linear relaxation rate, say, 40% greater,
$\eta_{k}=6.24\cdot 10^{4}$ s-1, then from the fit one gets
$\eta_{nl}N_{k}(0.8\,\mu\rm{s})=4.64\cdot 10^{5}$ s-1, $\eta_{nl}=1.4\cdot
10^{-11}$ s-1 and $N_{k}(0.8\,\mu\rm{s})\simeq 3.3\cdot 10^{16}$. We see that
the accuracy of the threshold does not seriously affect the nonlinear damping
term due to relatively small value of linear damping.

Let us now compare experiment and theory. The theoretical formula for the
coefficient of nonlinear radiation damping can be expressed in the form:

$\eta_{nl}^{(theor)}\simeq\xi_{R}\cdot 2\pi\hbar Q\frac{V^{2}_{k}}{v_{R}},$
(5)

where $V_{k}$ is the coupling coefficient for the parametric pair with the
pump field in the resonator cavity, in other words, it is proportional to an
effective magnetic moment $\hbar{\partial\omega_{n,k}}/{\partial H}$ of
excited wave. For nuclear spin waves one has AOSY ; OzhYak :

$V_{k}=-\frac{1}{2}\frac{\partial\omega_{n,k}}{\partial
H}=\frac{\omega_{n}^{2}}{4\omega_{n,k}}\frac{\gamma^{4}(H_{\Delta,hf})^{2}(2H+H_{D})}{\omega_{e,k}^{4}}.$
(6)

The factor $\xi_{R}$ in Eq.(5) depends on the geometry of resonator cavity.
For a rectangular resonator cavity one has $\xi_{R}=1$. For a helical
resonator a compression of half wavelength $\lambda/2$ to the length of helix
$l$ occrus and it can result in $\xi_{R}\sim\lambda/2l$.

Let us estimate theoretical nonlinear radiation damping, Eq.(5) for the
experiment shown in Fig.4, using the following parameters:
$\omega_{n}=2\pi\cdot 666$ MHz, $H_{D}=0$, $H^{2}_{\Delta,hf}=6.4/T[K]$ kOe2,
$l\sim 1$ cm. One gets: $\eta_{nl}^{(theor)}\sim 0.6\cdot 10^{-11}$ which of
the order of magnitude is in a good agreement with the obtained experimental
result.

###  III.2 Magnetoelastic Waves

Let us consider the experimental results shown in Fig.3 for magnetoelastic
waves. The linear relaxation rate calculated from the threshold of parallel
pumping in this case is $\eta_{k}=3.2\cdot 10^{5}$ s-1. From the mean-square
fit using formula (4) with $t_{0}=0.4$ $\mu\rm{s}$ one gets: 1)
$\eta_{nl}N_{k}(0.4\,\mu\rm{s})=0.55\cdot 10^{6}$ s-1,
$N_{k}(0.4\,\mu\rm{s})\simeq 2.6\cdot 10^{16}$ for $P/P_{c}=13.2$ and 2)
$\eta_{nl}N_{k}(0.4\,\mu\rm{s})=0.94\cdot 10^{6}$ s-1,
$N_{k}(0.4\,\mu\rm{s})\simeq 4.4\cdot 10^{16}$ for $P/P_{c}=52.3$. For both
cases we obtain the same experimental coefficient of nonlinear radiation
damping $\eta_{nl}=2.1\cdot 10^{-11}$ s-1.

Figure 5: Magnetic field dependence for the nonlinear radiation damping
coefficient of magnetoelastic waves in FeBO3 at $T=77$ K and
$\omega_{p}/2\pi=1109.7$ MHz. Solid line is the theoretical fit.

In order to derive theoretical estimate, we find:

$V_{k}=-\frac{1}{2}\frac{\partial\omega_{me,k}}{\partial
H}=\frac{(c_{e}k)^{2}}{4\omega_{me,k}}\frac{\gamma^{4}(H_{\Delta,ef})^{2}(2H+H_{D})}{\omega_{e,k}^{4}}.$
(7)

Thus, using Eq.(5) and the following parameters for the iron borate:
$c_{e}\simeq 4.8\cdot 10^{5}$ cm/s, $H_{\Delta,ef}\simeq 2$ kOe,
$H_{\Delta,me}\simeq 2.2$ kOe, $H_{D}\simeq 100$ kOe, $\alpha\simeq 0.08$
Oe$\cdot$cm, one gets: $\eta_{nl}^{(theor)}\sim 2.9\cdot 10^{-11}$ s-1, which
of the order of magnitude is in a good agreement with the obtained
experimental result.

Dots in Fig. 5 show the magnetic field dependence of experimentally obtained
coefficient of nonlinear radiation damping. The solid line represents the
theoretical prediction of the field dependence. We see a perfect fit within
two orders of magnitude of experimental data for nonlinear radiation damping.

##  IV conclusion

In this work we have experimentally confirmed that the nonlinear radiation
damping is the main mechanism of parametic instability restriction during
parallel microwave pumping of two different types of normal magnetic
oscillations, nuclear spin waves and magnetoelastic waves in different
antiferromagnets. The obtained results are in a good agreement with the theory
by the field and overthreshold dependencies and are of the order of magnitude
of the theoretical prediction. We believe that the nonlinear radiation damping
is a common feature of parallel pumping technique and it can be detected by
the purposeful experiments with other types of normal magnetic oscillations in
magnetoordered systems. For example, a specific radiation after turning off
the pump of spin waves in YIG has already been observed in Ref.ZhitMel  and it
was not explained in the framework of small sample approximation.

## References

  * (1) H. Suhl, J. Phys. Chem. Solids 1, 209 (1957).
  * (2) E. Schoemann, Phys. Rev. 116, 827 (1959).
  * (3) V. E. Zakharov, V. S. L’vov, and S. S. Starobinets, Usp. Fiz. Nauk 114, 609 (1974) [Sov. Phys.- Usp. 17, 896 (1975)].
  * (4) V. S. L’vov, Wave Turbulence under Parametric Excitation (Springer-Verlag, Berlin, 1994).
  * (5) A. G. Gurevich and G. A. Melkov, Magnetization Oscillations and Waves (CRC Press, Boca Raton, 1996).
  * (6) V. L. Safonov, Nonequilibrium Magnons (Wiley-VCH, Weinheim, 2013).
  * (7) V. L. Safonov, J. Magn. Magn. Mater. 97, L1 (1991).
  * (8) V. L. Safonov and H. Yamazaki, J. Magn. Magn. Mater. 161, 275 (1996).
  * (9) A. V. Andrienko and V. L. Safonov, Pis’ma Zh. Eksp. Teor. Fiz. 60, 787 (1994) [JETP Lett. 60, 800 (1994)].
  * (10) A. V. Andrienko and V. L. Safonov, Pis’ma Zh. Eksp. Teor. Fiz. 62, 147 (1995) [JETP Lett. 62, 162 (1995)].
  * (11) A. V. Andrienko, V. L. Safonov, and H. Yamazaki, J. Phys. Soc. Jpn. 67, 2893 (1998).
  * (12) P. G. de Gennes, P. A. Pincus, F. Hartmann-Boutron, and J. M. Winter, Phys. Rev. 129, 1105 (1963).
  * (13) A. V. Andrienko, V. I. Ozhogin, V. L. Safonov, and A. Yu. Yakubovskii, Usp. Fiz. Nauk, 161, 1 (1991) [Sov. Phys.- Usp. 34, 843 (1991)].
  * (14) A. V. Andrienko and L. V. Podd’yakov, Zh. Eksp. Teor. Fiz. 95, 2117 (1989) [Sov. Phys.- JETP 68, 1224 (1989)].
  * (15) A. V. Andrienko, Zh. Eksp. Teor. Fiz. 101, 1644 (1992) [Sov. Phys.- JETP 74, 876 (1992)].
  * (16) V. I. Ozhogin and A. Yu. Yakubovskii, Zh. Eksp. Teor. Fiz. 67, 287 (1974) [Sov. Phys.- JETP 40, 144 (1975)].
  * (17) V. S. Zhitnyuk and G. A. Melkov, Zh. Eksp. Teor. Fiz. 75, 1755 (1978) [Sov. Phys.- JETP 48, 884 (1978)].