# Shock-driven synchrotron radio emission from the 2021 outburst of RS
Ophiuchi

Nayana A.J.,1,2 G.C. Anupama,2 Nirupam Roy3 Dipankar P. K. Banerjee4 Kulinder
Pal Singh5 Sonith L.S.,2 and U. S. Kamath2
1Department of Astronomy, University of California, Berkeley, CA 94720-3411,
USA.
2Indian Institute of Astrophysics, II Block, Koramangala, Bangalore 560034,
India.
3Department of Physics, Indian Institute of Science, Bangalore 560012, India.
4Astronomy & Astrophysics Division, Physical Research Laboratory, Ahmedabad
380009, India.
5 Department of Physical Sciences, IISER Mohali, Knowledge City, Sector 81,
Manauli PO, SAS Nagar, Punjab 140306, India E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We present low-frequency radio observations of the Galactic symbiotic
recurrent nova RS Ophiuchi during its 2021 outburst. The observations were
carried out with the upgraded Giant Metrewave Radio Telescope (uGMRT) spanning
a frequency range of 0.15$-$1.4 GHz during 23$-$287 days post the outburst.
The average value of the optically thin spectral index is $\alpha\sim$ $-$0.4
($F_{\nu}\propto\nu^{\alpha}$), indicating a non-thermal origin of the radio
emission at the observed frequencies. The radio light curves are best
represented by shock-driven synchrotron emission, initially absorbed by a
clumpy ionized circumbinary medium. We estimate the mass-loss rate of the red
giant companion star to be $\dot{M}\sim$ 7.5 $\times$ 10-8 $M_{\odot}$ yr-1
for an assumed stellar wind velocity of 20 km/s. The 0.15–1.4 GHz radio light
curves of the 2021 outburst are systematically brighter than those of the 2006
outburst. Considering similar shock properties between the two outbursts, this
is indicative of a relatively higher particle number density in the
synchrotron emitting plasma in the current outburst.

###### keywords:

Nova: general — Nova: RS Ophiuchi — radiation mechanisms: non-thermal — radio
continuum: general

††pubyear: 2023††pagerange: Shock-driven synchrotron radio emission from the
2021 outburst of RS Ophiuchi–Shock-driven synchrotron radio emission from the
2021 outburst of RS Ophiuchi

##  1 Introduction

RS Ophiuchi is a Galactic symbiotic recurrent nova system composed of a
primary white dwarf that accretes matter from its red giant companion. The
nova has undergone seven reported outbursts, at somewhat irregular intervals,
in the years 1898, 1933, 1958, 1967, 1985, 2006, and 2021 (Schaefer, 2010, and
references therein). The 2021 outburst was discovered on 2021 August 08.93
(UT) in the optical bands
111http://ooruri.kusastro.kyoto-u.ac.jp/mailarchive/vsnet-alert/26131. The
nova brightened to a visual magnitude of 4 mag, which is $\sim$ 1000 times
brighter than its quiescent magnitude. Subsequently, the nova showed luminous
emission over a wide range of frequencies from radio to the Very High Energy
(VHE) $\gamma$-rays.

$\gamma$-ray emission was detected with Fermi-LAT (Cheung et al., 2022; Zheng
et al., 2022), H.E.S.S (H. E. S. S. Collaboration et al., 2022), and MAGIC
(Acciari et al., 2022) telescopes in the GeV-TeV energy regime, marking the
first detection of VHE emission from a nova. The GeV to TeV emission has been
interpreted under either a hadronic (H. E. S. S. Collaboration et al., 2022;
Acciari et al., 2022; Zheng et al., 2022; Diesing et al., 2023) or a lepto-
hadronic (De Sarkar et al., 2023) particle acceleration scenario. The
evolution of the hard X-ray emission, arising from shock interaction between
the nova ejecta and the red giant wind, was found to be similar to the 2006
outburst (Page et al., 2022). On the other hand, the evolution in the soft
X-ray, during the super soft state (SSS), was found to be fainter compared to
the 2006 outburst (Page et al., 2022). High resolution X-ray spectra observed
during $t\sim 18-21$ days with Chandra and XMM Newton (Orio et al., 2022)
indicated the presence of multiple thermal components over a range of
temperatures $T=0.07-3.4$ keV.

Ness et al. (2023) analyzed the X-ray grating spectra of both the 2006 and
2021 outbursts and concluded the low SSS emission in 2021 to be due to higher
absorption in the line of sight compared to that of 2006. The early optical
spectrum of the 2021 outburst showed emission lines of full width at half
maximum $\sim$ 2900 km s-1 (Munari & Valisa, 2021b). The ejecta accelerated to
velocities $\sim$ 4700 km s-1 (Mikolajewska et al., 2021) in the first two
days followed by decreasing line velocities indicating deceleration (Pandey et
al., 2022). The optical photometric and spectroscopic evolution of the 2021
outburst presented by Munari & Valisa (2021a, 2022) indicate it to be very
similar to that of previous outbursts. The morpho-kinematic modeling of Hα
emission lines suggested a bipolar ejecta geometry (Pandey et al., 2022)
similar to that seen during the 2006 outburst (Rupen et al., 2008).

Radio emission was detected from RS Ophiuchi from $\sim$ 1 day post-outburst
at multiple frequencies (0.15 $-$ 35 GHz) (Williams et al., 2021; Sokolovsky
et al., 2021; Nayana et al., 2021; de Ruiter et al., 2023). Munari et al.
(2022) presented the European Very Large Baseline Interferometry Network (EVN)
observations of RS Oph at 5 GHz and reported the total extension of the radio-
emitting region to be $\sim 90$ mas at $t\sim 34$ days. Radio emission from
novae could be thermal free-free emission from hot ionized ejecta and/or non-
thermal synchrotron emission from shocks. Low-frequency radio emission are
synchrotron-dominated and are critical in understanding the properties of
circumbinary medium (CBM) and shocks.

In this paper, we present extensive low-frequency (0.15 $-$ 1.4 GHz) radio
observations and modelling of RS Ophiuchi during its 2021 outburst. The paper
is organized as follows: We present the radio observations and data reduction
in §2. Radio light curves and spectra are presented in §3 and the radio
emission model is explained in §5. We discuss our results in §6 and §7. We
summarize our findings in §8.

Throughout this paper, we use the time of nova outburst, 2021 Aug 08.5 UT
(Munari & Valisa, 2021a) as the reference time. All times $t$ mentioned are
with respect to this reference epoch.

##  2 UGMRT Observations and Data Reduction

uGMRT observations of RS Ophiuchi were carried out from 2021 August 31
($t\sim$ 23 days) to 2022 May 23 ($t\sim$ 287 days) under Director’s
Discretionary Time (DDTC204, DDTC205, DDTC207, DDTC211, and DDTC222) and
regular GTAC proposals (42-083). The observations were done in band-2
(125$-$250 MHz), band-3 (250$-$500 MHz), band-4 (550$-$850 MHz), and band-5
(1000$-$1460 MHz). The data were recorded with an integration time of 10
seconds in the standard continuum mode. We used a processing bandwidth of 400
MHz in bands 4 and 5, 200 MHz in band 3, and 100 MHz in band 2, split into
2048 channels. 3C286, 3C48, and 3C147 were used as flux density calibrators.
J1822$-$096 and J1743$-$038 were used as the phase calibrators.

The Astronomical Image Processing Software (aips; Greisen, 2003) was used for
data reduction, following standard radio continuum imaging analysis procedures
(Nayana et al., 2017). The data were initially inspected for non-working
antennae and radio frequency interference (RFI) affected channels and
corrupted data were flagged. The calibration was done using a central channel
and the solutions were applied to the entire band. The fully calibrated source
data were imaged using aips task IMAGR. We performed a few rounds of phase-
only self-calibration to improve the image quality.

We detected radio emission from RS Ophiuchi at all frequencies and the
emission was unresolved in all the maps. We estimated the flux density and
error of the target source by fitting a Gaussian at the source position using
aips task JMFIT. The size of the fitted Gaussian is consistent with that of a
point source in all the maps. Resolutions of the radio maps are typically 5
$\times$ 2, 7 $\times$ 5, 13 $\times$ 7, and 25 $\times$ 17 arcsec2 in
bands-5, 4, 3, and 2, respectively. In addition to the errors on flux
densities from task JMFIT, we add (in quadrature) a 10% error in bands-3, 4,
and 5 and a 15% error in band-2 to account for the calibration uncertainties.
The flux densities at multiple epochs are listed in Table 1.

##  3 Radio light curves and spectra

The uGMRT light curves at 1.36, 0.69, 0.44, and 0.15 GHz are presented in Fig.
1 (left panel). The 1.36 and 0.69 GHz light curves show a steady decline from
a possible maximum at $t\lesssim$ 20 days. The 0.44 GHz light curve peaks
slightly later (at $t\sim$ 35 days) than that of the 1.36 and 0.69 GHz light
curves. The first three epochs of 0.15 GHz flux measurements are possibly
indicating a rising phase, whereas the light curve is in the declining phase
at $t\sim 287$ days. However, we note that the error bars on these flux
measurements are large. de Ruiter et al. (2023) presented LOFAR (54 and 154
MHz) and MeerKAT (0.82 and 1.28 GHz) observations of RS Ophiuchi during
$t\sim$ 2$-$223 days post the outburst. The combined GMRT+LOFAR+MeerKAT
observations suggests the radio light curve evolution to be an initial steep
rise to peak followed by a steady decline where higher frequencies trace the
peak first (see Fig. 2). This kind of light curve evolution is similar to that
seen in radio supernovae (Weiler et al., 2002).

The 1.28 GHz light curve presented in de Ruiter et al. (2023) shows a plateau
between $t\sim 10$ to 40 days. Our band-5 light curve does not show such a
flattening towards the peak due to sparse temporal sampling, covering only two
measurements between 10 and 40 days.

The near-simultaneous spectral indices, $\alpha$ ($F\propto\nu^{\alpha}$), are
estimated at multiple epochs (see Fig. 1 and Table 2). The spectral index
between 1.36/0.69 GHz is $-$0.67 $\pm$ 0.22 at $t\sim$ 54 days and gradually
approaches a value $\alpha$ $\sim$ $-$ 0.03 by $t\sim$ 259 days. The average
value of spectral index between 0.69 and 0.44 GHz is $\sim-0.6$ during $t\sim
44-183$ days, and flattens to $\alpha\sim-0.2$ to $-0.3$ during $t\sim
204-232$ days. The 0.44/0.15 GHz spectral indices are $\sim$ 0.8 in the rising
phase of the light curve and then become $\alpha\sim-1$ at $t\sim 287$ days.
The negative spectral indices suggest that the radio emission is predominantly
non-thermal (synchrotron) at the observed frequencies. However, the
contribution from non-thermal processes seems to decrease at epochs $t\gtrsim
204$ days. This could possibly be due to the shock deceleration and the shock
expanding to a rarified CBM further away from the binary system.

##  4 Brightness Temperature

Synchrotron radio emission is characterized by high brightness temperature
($T_{\rm B}$ $\gtrsim$ 105 $-$ 106 K), significantly above the maximum
brightness temperature of a photo ionized thermal gas (104 K). We calculate
the brightness temperature of the observed radio flux densities using the
following equation (Chomiuk et al., 2021)

$\displaystyle\frac{T_{\rm B}}{(K)}=1200\left(\frac{F_{\nu}}{\rm
mJy}\right)\left(\frac{\nu}{\rm GHz}\right)^{-2}\left(\frac{\theta}{\rm
arcsec}\right)^{-2}$ (1)

Here, $F_{\nu}$ is the radio flux density at frequency $\nu$, and $\theta$ is
the angular size of the radio-emitting region. Estimation of the brightness
temperature requires knowing the velocity of the ejected material and hence
the geometry and kinematics of the ejecta. High-resolution HST and radio
imaging of both the 2006 and 2021 outbursts show that the nova ejecta is
bipolar (or equivalently hourglass) in shape (Munari et al., 2022; O’Brien et
al., 2006; Bode et al., 2007). Specifically, radio imaging of the 2021
eruption (Munari et al., 2022) supports a geometric model wherein the leading
edges of the bipolar lobes have a space velocity of 7550 km/s at $t\sim$ 34
days. In a bipolar morphology, since the velocity flow is expected to be
homologous (i.e., $v\propto r$), the leading edge of the bipolar lobes (having
the largest r) has the highest velocity, while the other material, especially
the constricted material at the waist, has lower velocities. As a rough
approximation, we assume the waist is constricted to the extent that the
material there has a space velocity of $\sim$ 3780 km/s i.e. half of that of
the leading edge which has a velocity of 7550 km/s. We thus adopt a
simplified, mean space velocity of 5660 km/s for the entire ejecta at $t\sim$
34 days. Mondal et al. (2018) reported a velocity deceleration index of the
nova ejecta $v\propto t^{-0.66}$ at $t\sim$ 5 to 70 days during the 2006
outburst of RS Ophiuchi. The velocity remains roughly constant up to $t\sim
247$ days after a slight increase around $t\sim$ 80 days. We use $\theta$ from
VLBI observations (Munari et al., 2022) at $t\sim$ 34 days and use the
velocity deceleration indices from Mondal et al. (2018) to calculate the
angular size at later epochs. The brightness temperatures corresponding to our
flux density measurements are $\sim 10^{5}-10^{8}$ K for an adopted distance
of 2.68${}^{+0.17}_{-0.15}$ Kpc (Gaia DR3; Bailer-Jones et al., 2021). Based
on the estimated spectral indices and high brightness temperatures, we
conclude that the radio emission is predominantly non-thermal at the observed
frequencies.

Table 1: Details of the uGMRT observations of RS Ophiuchi in its 2021
outburst.

Date of observation | Timea | Frequency | Flux density
---|---|---|---
(UT) | (Days) | (GHz) | (mJy)
2021 Sep 02.53 | 25.03 | 1.36 | 85.10 $\pm$ 8.60
2021 Sep 14.73 | 37.23 | 1.37 | 61.00 $\pm$ 6.30
2021 Sep 26.69 | 49.19 | 1.36 | 43.10 $\pm$ 4.70
2021 Sep 31.67 | 54.17 | 1.36 | 45.90 $\pm$ 4.60
2021 Oct 05.34 | 57.84 | 1.36 | 60.40 $\pm$ 6.00
2022 Feb 05.26 | 180.76 | 1.36 | 13.10 $\pm$ 1.31
2022 Apr 25.03 | 259.53 | 1.36 | 09.90 $\pm$ 0.99
2021 Aug 31.67 | 23.17 | 0.69 | 77.50 $\pm$ 8.00
2021 Sep 17.67 | 40.17 | 0.69 | 78.00 $\pm$ 7.90
2021 Sep 21.72 | 44.22 | 0.69 | 63.00 $\pm$ 6.40
2021 Oct 01.58 | 54.08 | 0.69 | 72.10 $\pm$ 7.70
2021 Oct 08.41 | 60.91 | 0.66 | 58.90 $\pm$ 6.10
2021 Nov 15.42 | 98.92 | 0.64 | 43.30 $\pm$ 4.40
2021 Nov 26.25 | 109.75 | 0.69 | 23.90 $\pm$ 2.60
2021 Dec 06.15 | 119.65 | 0.69 | 26.90 $\pm$ 3.00
2021 Dec 18.21 | 131.71 | 0.69 | 25.00 $\pm$ 2.50
2021 Dec 31.26 | 144.76 | 0.69 | 21.30 $\pm$ 2.30
2022 Jan 10.08 | 154.58 | 0.69 | 16.50 $\pm$ 1.80
2022 Feb 08.19 | 183.69 | 0.69 | 16.49 $\pm$ 1.71
2022 Mar 01.04 | 204.54 | 0.69 | 13.36 $\pm$ 1.37
2022 Mar 28.95 | 232.45 | 0.69 | 11.55 $\pm$ 1.17
2022 Apr 24.98 | 259.48 | 0.62 | 10.14 $\pm$ 1.02
2021 Sep 05.51 | 28.01 | 0.40 | 73.10 $\pm$ 11.20
2021 Sep 16.70 | 39.20 | 0.46 | 71.70 $\pm$ 7.80
2021 Sep 21.38 | 43.88 | 0.46 | 85.10 $\pm$ 12.80
2021 Sep 30.62 | 53.12 | 0.45 | 80.00 $\pm$ 8.10
2021 Oct 03.42 | 55.92 | 0.46 | 75.80 $\pm$ 11.50
2021 Nov 15.49 | 98.99 | 0.40 | 48.30 $\pm$ 05.60
2021 Nov 26.33 | 109.83 | 0.44 | 38.43 $\pm$ 3.85
2021 Dec 07.46 | 120.96 | 0.45 | 37.15 $\pm$ 3.74
2021 Dec 18.46 | 131.96 | 0.43 | 29.68 $\pm$ 3.07
2021 Dec 31.34 | 144.84 | 0.44 | 28.06 $\pm$ 2.86
2022 Jan 10.32 | 154.82 | 0.44 | 23.91 $\pm$ 2.41
2022 Feb 08.00 | 183.81 | 0.46 | 18.40 $\pm$ 1.93
2022 Mar 01.12 | 204.62 | 0.44 | 14.46 $\pm$ 1.48
2022 Mar 29.04 | 232.54 | 0.44 | 13.12 $\pm$ 1.34
2022 May 22.85 | 287.35 | 0.44 | 08.93 $\pm$ 1.22
2021 Sep 30.53 | 53.03 | 0.147 | 33.20 $\pm$ 8.20
2021 Oct 03.34 | 55.84 | 0.147 | 31.40 $\pm$ 6.60
2021 Oct 08.38 | 60.88 | 0.146 | 48.90 $\pm$ 10.80
2022 May 22.91 | 287.41 | 0.147 | 26.80 $\pm$ 5.10

a This column indicates the time from the nova outburst i.e., 2021 Aug 08.5 UT
(Munari & Valisa, 2021a).

Table 2: Near simultaneous radio spectral indices of RS Ophiuchi in its 2021
outburst.

Timea |  | Spectral Indices |
---|---|---|---
(Day) | (1.4 and 0.69 GHz) | (0.69 and 0.44 GHz) | (0.44 and 0.15 GHz)
54.13 |  $-$0.67 $\pm$ 0.22 | - | -
182.22 |  $-$0.34 $\pm$ 0.15 | - | -
259.51 |  $-$0.03 $\pm$ 0.18 | - | -
44.05 | - |  $-$0.74 $\pm$ 0.26 | -
53.60 | - |  $-$0.24 $\pm$ 0.32 | -
98.96 | - |  $-$0.23 $\pm$ 0.27 | -
109.79 | - |  $-$1.06 $\pm$ 0.21 | -
120.31 | - |  $-$0.76 $\pm$ 0.27 | -
131.84 | - |  $-$0.36 $\pm$ 0.25 | -
144.80 | - |  $-$0.61 $\pm$ 0.26 | -
154.70 | - |  $-$0.82 $\pm$ 0.24 | -
183.75 | - |  $-$0.27 $\pm$ 0.32 | -
204.58 | - |  $-$0.18 $\pm$ 0.30 | -
232.49 | - |  $-$0.28 $\pm$ 0.28 | -
53.07 | - | - | 0.79 $\pm$ 0.31
55.88 | - | - | 0.77 $\pm$ 0.45
287.38 | - | - |  $-$1.00 $\pm$ 0.26

a This column indicates the time from the nova outburst i.e., 2021 Aug 08.5 UT
(Munari & Valisa, 2021a).

Figure 1: Left panel: The uGMRT light curves of RS Ophiuchi in its 2021
outburst at 1.36, 0.69, 0.44, and 0.15 GHz. Right panel: The near simultaneous
spectral indices between 1.36/0.69 GHz, 0.69/0.44 GHz, and 0.44/0.15 GHz.

##  5 Radio emission model

Low-frequency radio light curves of RS Ophiuchi can be represented by adopting
the synchrotron emission model from Weiler et al. (2002). According to this
model, the outburst ejecta plows into the dense red giant wind driving a
forward shock. Particles are accelerated at the shock via diffusive shock
acceleration and emit synchrotron radiation. Initially, the emission is
suppressed via free-free absorption by a homogeneous and clumpy circumbinary
medium in the line of sight. The spectral and temporal evolution of radio flux
densities can be represented as

$\displaystyle F\rm(\nu,t)=K_{\rm 1}\left(\frac{\nu}{1\,\rm
GHz}\right)^{\alpha}\left(\frac{t}{20\,\rm days}\right)^{\beta}$
$\displaystyle\times\rm exp(-\tau_{\rm homog})\left[\frac{1-exp(-\tau_{\rm
clumps})}{\tau_{\rm clumps}}\right]$ (2)

Here, $K_{\rm 1}$ denotes the flux density normalization constant and $t$
denotes the time since the outburst in days. $\beta$ is the temporal decay
index of flux densities in the optically thin regime. $\tau_{\rm homog}$ and
$\tau_{\rm clumps}$ represent the optical depths due to homogeneous and clumpy
circumbinary medium, respectively.

$\displaystyle\tau_{\rm homog}=K_{\rm 2}\left(\frac{\nu}{1\,\rm
GHz}\right)^{-2.1}\left(\frac{t}{20\,\rm days}\right)^{\delta}$ (3)

$\displaystyle\tau_{\rm clumps}=K_{\rm 3}\left(\frac{\nu}{1\,\rm
GHz}\right)^{-2.1}\left(\frac{t}{20\,\rm days}\right)^{\delta^{{}^{\prime}}}$
(4)

In the above equations, $K_{2}$ and $K_{3}$ are the optical depths due to
homogeneous and clumpy circumbinary medium, respectively at 1 GHz on $t=$ 20
days. $\delta$ and $\delta^{{}^{\prime}}$ denote the temporal indices of
optical depths due to homogeneous and clumpy absorbing medium, respectively.

We model the radio flux density measurements using equations (2, 3, and 4)
keeping $K_{1}$, $K_{2}$, $K_{3}$, $\alpha$, $\beta$, $\delta$, and
$\delta^{{}^{\prime}}$ as free parameters. In addition to the uGMRT flux
measurements, data from de Ruiter et al. (2023) are also used. Since it is
likely that data at $t<$ 5 days are due to a residual synchrotron component
from previous outbursts (de Ruiter et al., 2023), only the data from $t>$ 5
days are included for modeling. The python emcee package (Foreman-Mackey et
al., 2013) that adopts Markov chain Monte Carlo (MCMC) method is used to
execute the fit. We chose 5000 steps and 32 walkers to explore the parameter
space to get the best-fit parameters (68% confidence interval). The best-fit
model along with the observed flux values are plotted in Fig. 2 and the best-
fit parameters are listed in Table 3. The reduced chi-square value for the
best-fit model is $\chi_{\mu}^{2}=3.7$. We also modeled the data including the
flux measurements at t < 5 days and the fit resulted in a much higher reduced
$\chi_{\mu}^{2}=6.4$.

The best-fit optically thin spectral and temporal indices are
$\alpha=-0.37^{+0.05}_{-0.04}$ and $\beta=-1.04^{+0.03}_{-0.03}$,
respectively, roughly similar to the values seen in radio supernovae (Weiler
et al., 2002). The best-fit optical depth normalization parameters due to
uniform and clumpy medium are $K_{\rm 2}=1.68\times 10^{-3}$ and $K_{\rm
3}=0.68$, respectively. These values suggest that a clumpy medium in the line
of sight is dominating the absorption processes. We repeated the modeling
using only GMRT flux measurements and the best-fit parameters are
$K_{1}=171^{+11}_{-9}$ $K_{2}=0.04^{+0.02}_{-0.02}$,
$K_{3}=1.33^{+0.35}_{-0.25}$, $\alpha=-0.50^{+0.06}_{-0.06}$,
$\beta=-1.16^{+0.04}_{-0.04}$, $\delta=-1.74^{+0.52}_{-0.53}$, and
$\delta^{{}^{\prime}}=-3.48^{+0.43}_{-0.38}$. The values of $K_{1}$, $\alpha$,
$\beta$, and $\delta^{{}^{\prime}}$ are roughly similar if we use either GMRT
or GMRT+LOFAR+MeerKAT data. However, the value of optical depth normalization
parameters ($K_{2}$, $K_{3}$) and $\delta$ are significantly different. This
indicates that including early time data is critically important to pin down
the dominant absorption processes.

While the synchrotron model broadly represents the spectral and temporal
behavior of radio light curves, the plateauing towards the peak is not well
fitted in the model.

We thus modeled the data in a two-component scenario and found that the fit
improved with a reduced chi-square $\chi_{\mu}^{2}=1.5$. The best-fit model
along with the observed flux density measurements are shown in Fig. 3. The
best-fit parameters that represent each of the synchrotron components are;
$K_{1}\sim 60.3$, $K_{2}\sim 1.2\times 10^{-4}$, $K_{3}\sim 2.0\times
10^{-3}$, $\alpha\sim-0.12$, $\beta\sim-0.8$, $\delta\sim-4.8$, and
$\delta^{{}^{\prime}}\sim-4.4$ for the first component, and $K_{1}\sim 18.4$,
$K_{2}\sim 4.6\times 10^{-4}$, $K_{3}\sim 1.8,\alpha\sim-1.9$,
$\beta\sim-2.5$, $\delta\sim-1.1$, and $\delta^{{}^{\prime}}\sim-4.5$ for the
second component. The flattening or marginal flux density enhancement close to
the peak of the light curves is well represented in the two-component model.
The first synchrotron component peaks at $t\sim 12-20$ days, whereas the
contribution from the second component is maximum at $t\sim 24-45$ days. The
flux density enhancement in the light curves above the prediction of a single
(first component) synchrotron model during $t\sim 24-45$ days could be due to
a density enhancement in the CBM, as also seen in the 2006 outburst (Kantharia
et al., 2007). Alternatively, it could also be a consequence of the complex
velocity changes of the ejecta as the shock decelerates, as seen for example
during the first 50 days of the 2006 eruption (Fig.2; Das et al., 2006).
Velocity changes would result in fluctuations in the kinetic temperature of
the shocked gas which would, in turn, affect the non-thermal emission.

An alternative scenario could involve the emergence of additional radio
outflow components to the visibility during $t\sim$ 24$-$45 days. The EVN
observations of RS Oph at $t\sim$ 34.3 days (Munari et al., 2022) showed three
components; one compact central component (CC) and two lobes towards the east
and west of the CC (see Fig 2 of Munari et al., 2022). According to the model
proposed by Munari et al. (2022), the eastern arc (EA) remains hidden from the
observer at early times due to the high free-free opacity caused by the
density enhancement in the orbital plane (DEOP). At times $t>t_{\tau}$ (where
$t_{\tau}$ is the time required for EA to expand from the central binary such
that the optical depth due to DEOP is $\sim$ 1), the radio emission from EA
becomes detectable. As more and more of EA component becomes visible to the
observer, the detected flux density will increase. In the 2006 outburst, it
was believed that EA became visible during $t\sim$ 21–51 days (O’Brien et al.,
2006, 2008; Rupen et al., 2008; Sokoloski et al., 2008). In the 2021 outburst,
the EA is in the emerging phase (outer arc is visible) at t $\sim$ 34.3 days
(Munari et al., 2022). Considering the similarities between the 2006 and 2021
outbursts in terms of shell velocity and radio morphology, it is very likely
that the enhanced flux densities during $t\sim$ 24$-$45 days are due to the
eastern lobe becoming visible to the observer.

Kantharia et al. (2007) presented low-frequency radio light curves of the 2006
outburst of RS Ophiuchi. We use these flux density measurements and model them
(in a single component) following the same method as discussed here. The best-
fit values are given in Table 3. We use the results from the single-component
modeling and attempt a detailed comparison of the best-fit parameters with
those of the 2006 outburst in §6.

Table 3: Best fit parameters of the synchrotron emission model.

2021 outburst | 2006 outburst
---|---
$K_{1}=141.25^{+6.97}_{-7.17}$ | $K_{1}=79.77^{+14.27}_{-18.55}$
$K_{2}=1.68^{+0.78}_{-0.55}\times$ 10-3 | $K_{2}=0.07^{+0.07}_{-0.06}$
$K_{3}=0.68^{+0.10}_{-0.10}$ | $K_{3}=0.53^{+0.37}_{-0.29}$
$\alpha=-0.37^{+0.05}_{-0.04}$ | $\alpha=-0.73^{+0.14}_{-0.13}$
$\beta=-1.04^{+0.03}_{-0.03}$ | $\beta=-1.19^{+0.14}_{-0.09}$
$\delta=-5.74^{+0.28}_{-0.30}$ | $\delta=-2.32^{+0.66}_{-0.61}$
$\delta^{{}^{\prime}}=-3.05^{+0.10}_{-0.12}$ | $\delta^{{}^{\prime}}=-3.51^{+0.58}_{-0.65}$
$\chi^{2}=3.70$ | $\chi^{2}=2.72$

$K_{1}$, $K_{2}$, $\alpha$, $\beta$, $\delta$, and $\delta^{{}^{\prime}}$ are
best-fit parameters of radio emission model defined by equations 2, 3 and 4
(see §5). The modeling (of 2021 outburst) was done including flux density
measurements from our uGMRT observations and the ones reported by de Ruiter et
al. (2023) at $t>$ 5 days.

Figure 2: Low-frequency radio light curves of RS Ophiuchi in its 2021
outburst. The data points include flux density measurements from our uGMRT
observations and the measurements reported by de Ruiter et al. (2023) at $t>$
5 days. The purple star symbols denote the uGMRT flux density measurements and
the black circles denote the flux measurements reported by de Ruiter et al.
(2023). The orange solid curves represent the best-fit single component
synchrotron emission model (Weiler et al., 2002) as described by equations 2,
3, and 4 (see §5). The bottom panels show the residuals defined as
$F_{o}-F_{\rm mod}$, where $F_{\rm o}$ is the observed flux density and
$F_{\rm mod}$ is the modeled flux value. Figure 3: Low-frequency radio light
curves of RS Ophiuchi in its 2021 outburst. The data points include flux
density measurements from our uGMRT observations and the measurements reported
by de Ruiter et al. (2023) at $t>$ 5 days. The purple star symbols denote the
uGMRT flux density measurements and the black circles denote the flux
measurements reported by de Ruiter et al. (2023). The grey and blue dashed
curves represent the first and second components of the synchrotron emission
model (Weiler et al., 2002), respectively. The orange solid curves represent
the total fit which is the sum of two components. The second component of the
0.05 GHz light curve is out of the plot range (the modeled flux values are too
low; $F\ll$ 1 mJy).

###  5.1 Mass-loss rate

Mass-loss rate ($\dot{M}$) of the red giant companion can be estimated from
the modeled radio light curves (Weiler et al., 1986, 2002).

$\displaystyle\frac{\dot{M}}{v_{w1}}=3.0\times 10^{-6}\left<\tau_{\rm
eff}^{0.5}\right>m^{-1.5}\left(\frac{v_{\rm i}}{10^{4}\rm
km\,s^{-1}}\right)^{1.5}$ $\displaystyle\times\left(\frac{t_{\rm i}}{45\rm
days}\right)^{1.5}\left(\frac{t}{t_{\rm
i}}\right)^{1.5m}\left(\frac{T}{10^{4}\rm K}\right)^{0.68}$ (5)

Here, $v_{w1}$ is the stellar wind velocity of the red giant companion in
units of 10 km s-1. $v_{\rm i}$ is the ejecta velocity at $t_{\rm i}$ days
post-outburst. $m$ and $T$ represent the shock deceleration parameter
($m=\delta/3$) and wind electron temperature, respectively. $\tau_{\rm eff}$
is the effective optical depth due to a uniform and clumpy absorbing medium
(Weiler et al., 2002).

$\displaystyle\left<\tau_{\rm eff}^{0.5}\right>=0.67\left[(\tau_{\rm
homog}+\tau_{\rm clumps})^{1.5}-\tau_{\rm homog}^{1.5}\right]\tau_{\rm
clumps}^{-1}$ (6)

We estimate the mass-loss rate to be $\dot{M}=$ 7.5${}^{+1.5}_{-1.6}$ $\times$
10-8 $M_{\odot}$ yr-1 for a wind velocity of 20 km s-1 and $T=$ 20,000 K
(Weiler et al., 2002). We use $v_{\rm i}\sim$ 5660 km s-1 at $t_{\rm i}=$ 34
days (see §4). The derived mass-loss rate is in agreement with the typical
mass-loss rates (10-9 $-$ 10-6 $M_{\odot}$ yr-1) of red giant branch (RGB)
stars (van Loon et al., 2006; Origlia et al., 2007; McDonald & van Loon, 2007)
and with that presented in (de Ruiter et al., 2023). The above-mentioned
formulation of mass-loss rate is originally developed for supernovae (Weiler
et al., 2002) assuming a spherical geometry. In the case of RS Oph, the ejecta
is bipolar and the estimated $\dot{M}$ value could be a rough approximation.

##  6 Comparison between the radio evolution of 1985, 2006, and 2021 outbursts
of RS Ophiuchi

Radio emission was detected and monitored from the 1985 (Hjellming et al.,
1986) and 2006 (Kantharia et al., 2007; Eyres et al., 2009) outbursts of RS
Ophiuchi. In this section, we compare the properties of radio emission from
previous outbursts with that of 2021 to understand the evolution of the system
in terms of shock energetics and the circumbinary medium.

Hjellming et al. (1986) presented radio observations spanning $t\sim$ 29 $-$
370 days post the 1985 outburst at frequencies 1.49, 4.85, 14.94, and 22.46
GHz with the Very Large Array (VLA). Radio emission from the 2006 outburst of
RS Ophiuchi was detected at frequencies 0.24, 0.33, 0.61, 1.49, 4.89, 6.0,
14.96, and 22.48 GHz (Narumi et al., 2006; Eyres et al., 2009; Kantharia et
al., 2007). In both outbursts, two spectral components were identified, one
with a negative spectral index between 1.4/5 GHz and below, and the other with
a positive spectral index above 5 GHz. The GMRT observations of the 2006
outburst were carried out by Kantharia et al. (2007) during $t\sim$ 13 $-$ 110
days post-discovery at frequencies 1.39, 1.28, 1.06, 0.61, 0.33, 0.24, and
0.15 GHz. Radio emission was detected at all frequencies except at $\nu=$ 0.15
GHz at $t\sim$ 110 days with a 3$\sigma$ flux density limit of 27 mJy.

We compile all published low-frequency ($\nu\leq$ 1.4 GHz) radio observations
of 1985 and 2006 outbursts to investigate the evolution of shock-driven
synchrotron emission over successive outbursts. Fig. 4 shows the radio light
curves at $\nu\sim$ 0.15 $-$ 1.4 GHz of RS Ophiuchi in its 2006 and 2021
outbursts and Fig. 5 shows the 1.4 GHz light curves of 1985, 2006, and 2021
outbursts.

When we compare the low-frequency radio light curves of the 2006 and 2021
outbursts, two characteristics are particularly striking.

  1. 1.

The radio flux densities during the 2021 outburst are systematically higher
than those of the 2006 outburst.

  2. 2.

Radio emission at 0.15 GHz is detected from $t\sim$ 24 to 287 days at 9 epochs
during the 2021 outburst. Radio observations at 0.15 GHz were carried out only
at a single epoch ($t\sim$ 110 days) in the 2006 outburst which resulted in
non-detection with a 3-sigma flux density limit of $F<$ 27 mJy.

###  6.1 Brighter synchrotron emission during the 2021 outburst

The observed peak flux densities of the bands 5, 4, and 3 light curves in the
2021 outburst are $\sim$ 85, 78, and 85 mJy, respectively. The corresponding
values in the 2006 outburst are $\sim$ 57, 49, and 57 mJy, respectively
(Kantharia et al., 2007), about 1.5–1.6 times lower. If we compare the flux
densities at higher frequencies between 2021 (Sokolovsky et al., 2021) and
2006 outbursts (Eyres et al., 2009), this trend seems to be broadband.

The key dependencies on the synchrotron emissivity [$J(\nu)$] are (equation
8.89; Longair, 2011)

$J(\nu)\propto N_{\rm 0}\,B^{(p+1)/2}\,\nu^{-(p-1)/2}$ (7)

where the power-law distribution of electron energies is represented as
$N(E)dE=$ $N_{\rm 0}E^{-p}dE$. Here, $N(E)dE$ is the number density of
electrons of energies between $E$ and $E+dE$. Thus a higher synchrotron radio
flux density could be due to more synchrotron emitting electrons in the
emitting region and/or due to a higher magnetic field.

Assuming the amplified post-shock magnetic field is in equipartition with the
relativistic electrons in the radio-emitting region (Chevalier & Fransson,
2006; Chugai, 2010), the magnetic field can be written as (Chomiuk et al.,
2012)

$B=\sqrt{8\pi\epsilon_{\rm B}\mu m_{\rm H}n_{\rm CBM}v^{2}}$ (8)

Here, $n_{\rm CBM}$ is the particle number density in the CBM and $v$ is the
shock velocity. Munari et al. (2022) reported the size and expansion velocity
of the radio-emitting region at $t\sim$ 34 days in the 2021 outburst and
mentioned that the size is remarkably similar in the case of 2006 outburst at
similar epochs. The hard X-ray emission and its temporal evolution are also
reported to be broadly similar between the two outbursts (Page et al., 2022).
Both these observations are indicative of similar shock velocities during the
2006 and 2021 outbursts. Thus the higher synchrotron brightness in the 2021
outburst could be due to a higher density of the CBM. The synchrotron radio
flux density can be represented in terms of electron number density in the
synchrotron plasma as $F\propto$ $n_{\rm e}^{(p+5)/4}$ where $p$ is related to
the optically thin spectral index as; $\alpha=(p-1)/2$. Thus a 1.5$-$1.6 times
increase in synchrotron flux densities during the 2021 outburst compared to
that of 2006 could indicate a $\sim$ 30% increase in $n_{\rm e}$ considering
the shock velocity and magnetic field are the same during both outbursts.

We detect 0.15 GHz radio emission at $t\sim 52,55,60,287$ days post the 2021
outburst of RS Ophiuchi. de Ruiter et al. (2023) reported radio emission at
0.15 GHz at $t\sim$ 24, 34, 49, 65, and 223 days post outburst. Observations
at this frequency at a single epoch ($t\sim$ 110 days) during the 2006
outburst resulted in non-detection with a 3$\sigma$ flux density limit of 27
mJy. However, during the 2021 outburst, the best fit model predicts a flux
density of $\sim$ 44 mJy at 110 days, which is significantly higher than the
flux density limit during the 2006 outburst. This could either be indicative
of an early turn-on of the 0.15 GHz radio emission in the 2021 outburst or be
due to the intrisically brighter synchrotron emission in 2021 compared to that
of 2006.

Figure 4: Radio light curves of RS Ophiuchi in its 2006 and 2021 outburst. The
uGMRT radio flux density measurements of the 2021 outburst are marked as red
star symbols. The flux densities of the 2021 outburst reported in de Ruiter et
al. (2023) are marked as orange star symbols. The blue and black filled
circles and downward triangles represent the flux densities and limits,
respectively during the 2006 outburst (Kantharia et al., 2007). Figure 5:
Radio light curves of RS Ophiuchi in its 1985, 2006, and 2021 outburst. The
uGMRT radio flux density measurements of the 2021 outburst are marked as red
star symbols. The flux densities of the 2021 outburst reported in de Ruiter et
al. (2023) are marked as orange star symbols. The blue, cyan, and magenta
filled circles represent the flux densities during the 2006 outburst
(Kantharia et al., 2007). Grey diamonds denote the flux density measurements
from the 1985 outburst (Hjellming et al., 1986).

We now look at the published 1.4 GHz radio observations during the 1985
outburst along with those of 2006 and 2021 in Fig. 5. The radio emission is
brightest in the 2021 outburst whereas the 2006 and 1985 outbursts show
roughly similar flux values. The time of peak of the 1.4 GHz light curve
($t_{\rm p}$) is not well constrained in the 1985 outburst due to sparse
sampling of the light curve at $t<20$ days (see Fig 1 of Hjellming et al.,
1986). We assume the peak of the 1.4 GHz light curve in the 1985 outburst to
be $t_{\rm p}\lesssim$ 20 days, similar to that of the 2006 ($t_{\rm p}\sim$
15 days) and 2021 ($t_{\rm p}\sim$ 15 days) outbursts. The 0.325 GHz emission
was not detected from the 1985 outburst at $t\sim$ 48 days down to a flux
density limit of 5 mJy whereas 0.325 GHz emission was detected at $t\sim$ 38
days in the 2006 outburst with a flux density of 43.7 mJy. The non-detection
of 0.325 GHz emission are indications of higher foreground absorption in 1985
than in 2006 (Kantharia et al., 2007).

###  6.2 Evolution of optical depths due to uniform and clumpy absorbing
medium

We investigate the evolution of optical depth due to uniform and clumpy
components of the absorbing medium during the 2006 and 2021 outbursts from the
best-fit modeled parameters. Table 4 shows the values of $\tau_{\rm homo}$,
$\tau_{\rm clumps}$, and $\tau_{\rm eff}$ in the 2006 and 2021 outbursts at
the time of the peak of the 1.36, 0.69, 0.44, and 0.15 GHz light curves in
2021 outburst. The optical depth due to the inhomogeneous component of the
absorbing medium is systematically higher in 2021 than in 2006 suggesting that
the clumpiness in the circumbinary medium has increased from 2006 to 2021.

Table 4: Optical depth due to uniform and clumpy circumburst medium during the
2021 and 2006 outbursts.

Timea |  | 2021 |  |  | 2006 |
---|---|---|---|---|---|---
(Day) | $\tau_{\rm homo}$ | $\tau_{\rm clumps}$ | $\tau_{\rm eff}$ | $\tau_{\rm homo}$ | $\tau_{\rm clumps}$ | $\tau_{\rm eff}$
| ($\times$10-3) |  |  |  |  |
15.0 | 1.58 $\pm$ 0.46 | 0.38 $\pm$ 0.06 | 0.17 $\pm$ 0.03 | 0.04 $\pm$ 0.03 | 0.30 $\pm$ 0.21 | 0.18 $\pm$ 0.10
24.2 | 0.60 $\pm$ 0.19 | 0.44 $\pm$ 0.07 | 0.20 $\pm$ 0.03 | 0.06 $\pm$ 0.05 | 0.29 $\pm$ 0.21 | 0.20 $\pm$ 0.11
33.1 | 0.59 $\pm$ 0.20 | 0.69 $\pm$ 0.11 | 0.31 $\pm$ 0.05 | 0.11 $\pm$ 0.09 | 0.41 $\pm$ 0.33 | 0.30 $\pm$ 0.18
70.4 | 0.04 $\pm$ 0.02 | 0.52 $\pm$ 0.12 | 0.23 $\pm$ 0.05 | 0.15 $\pm$ 0.18 | 0.21 $\pm$ 0.24 | 0.25 $\pm$ 0.22

a This column indicates the time from the nova outburst i.e., 2021 Aug 08.5 UT
(Munari & Valisa, 2021a).

$\tau_{\rm homo}$ and $\tau_{\rm clumps}$ denote the optical depths due to
uniform and clumpy mediums, respectively. $\tau_{\rm eff}$ is the effective
optical depth (Weiler et al., 2002).

##  7 A comprehensive picture

Shock-driven synchrotron radio emission from multiple outbursts of RS Ophiuchi
shows differences in radio brightness. It is important to tie down these
differences with observational characteristics at other frequencies (optical,
X-ray, $\gamma$-rays) and theoretical models to build a comprehensive picture.

The optical light curves and evolution between 2006 and 2021 outbursts are
almost identical, indicating that the ejecta mass, the mass of the white
dwarf, and ejecta velocities are broadly the same (Page et al., 2022). The
hard X-ray emission ($>$ 1 keV) during the 2021 and 2006 outbursts are similar
in terms of brightness and temperature evolution (Orio et al., 2022; Page et
al., 2022) suggesting the gross properties of the shock interaction to be
similar. The results from VLBI observations also point to the same
conclusions. The shock radius at $t\sim$ 34 days in the 2021 outburst from
VLBI measurements (Munari et al., 2022) is remarkably the same as the shock
radius at similar epochs in the 2006 eruption (O’Brien et al., 2006, 2008;
Sokoloski et al., 2008; Rupen et al., 2008). Contrary to a very similar
evolution in optical, hard X-ray, and VLBI observations, the soft ($0.3-1$
keV) X-ray emission during the SSS phase was found to be very different. The
SSS emission during the 2021 outburst is $4-5$ times fainter than that of the
2006 outburst (Page et al., 2022).

A higher accretion rate can lead to a rapid increase in temperature of the WD
surface for a given accreted mass (Prialnik et al., 1982) and result in a
weaker nova and SSS emission. However, there are no observational signatures
of a difference in accretion rate during the quiescent intervals prior to 2006
and 2021 outbursts (Orio, 1993; Page et al., 2022). The quiescent X-ray
luminosities (0.3 $-$ 10 keV) are $\sim$ 6 $\times$ 1031 erg s-1, orders of
magnitude less compared to the expected value of $\sim$ 1036 erg s-1 (Osborne
et al., 2011; Anupama & Mikołajewska, 1999). A possible explanation for such
suppression in quiescent luminosity is the complex (partially covering)
absorbing medium in the system. Ness et al. (2023) compared the X-ray grating
spectra of 2006 and 2021 outbursts during the SSS phase and concluded that
higher absorption in the line of sight is the main reason for the observed
lower SSS emission in 2021.

The line of sight absorbing material is created due to aspherical mass-loss in
the binary system and could have a dependence on the orbital phase (Orlando et
al., 2009; Drake et al., 2009; Walder et al., 2008). The 1985, 2006, and 2021
outbursts occurred at phases 0.32, 0.26, and 0.72, respectively. Simulation
results from Booth et al. (2016) showed that a spiral accretion during the
quiescent phase builds up material such that there could be more material in
the line of sight during the 2021 outburst. The simulation shows a low density
at an orbital phase $\sim$ 0.26, which is during the 2006 outburst. The
authors also suggest that the simulated material is clumpy.

We re-visit the differences in the evolution of low-frequency radio light
curves in 2006 and 2021 outbursts in light of the results from Ness et al.
(2023). The radio emission is the brightest in the 2021 outburst with a peak
flux density of $\sim$ 87 mJy at 1.4 GHz on $t\sim$ 15 days post outburst (de
Ruiter et al., 2023). The 1.4 GHz light curve is sparsely sampled in the 2006
outburst and the peak flux and time from the best-fit model are $F_{\rm
p}\sim$ 61 mJy and $t\sim$ 15 days, respectively (Kantharia et al., 2007).
Assuming similar ejecta and shock properties, these peak flux densities
indicate that there are more particles in the synchrotron emitting plasma in
2021 than that of 2006 outburst.

The integrated X-ray count from SSS emission in the 2021 outburst during
$t\sim$ 30 to 86 days is 4$-$5 times lower than that of 2006 possibly due to
higher absorption (Ness et al., 2023). These two observational characteristics
point out that the medium in the vicinity of the binary system during the 2021
outburst is relatively denser compared to that of 2006 outburst.

##  8 Summary

We present low-frequency radio observations of the 2021 outburst of RS
Ophiuchi during 23 to 287 days post outburst at frequencies 0.15 $-$ 1.4 GHz.
Our main findings are the following:

  1. 1.

Radio emission at the observed frequencies are characterized by high
brightness temperatures (105 $-$ 108 K) with an average optically thin
spectral index of $\alpha\sim$ $-$ 0.4, consistent with non-thermal emission
process.

  2. 2.

The spectral and temporal evolution of the radio flux measurements are best
represented by synchrotron emission that arises due to the interaction of
external shock from the nova outburst with the dense wind of the red giant
companion. The emission is initially suppressed due to free-free absorption by
the uniform and clumpy foreground circumbinary medium.

  3. 3.

The mass-loss rate of the red giant has been estimated to be $\dot{M}=$ 7.5
${}^{+1.5}_{-1.6}$ $\times$ 10-8 $M_{\odot}$ yr-1 for a wind velocity of 20 km
s-1, consistent with the mass-loss rates of RGB stars (McDonald & van Loon,
2007).

  4. 4.

Comparison of light curves from the 2021 and 2006 outbursts in the 0.15$-$1.4
GHz radio bands shows that the peak flux densities in the 2021 outburst are
1.5 $-$ 1.6 times higher than that of the 2006 outburst values. Considering
similar shock and ejecta properties of these outbursts as supported by
multiwavelength observations (Page et al., 2022; Orio et al., 2022; Munari et
al., 2022), we interpret the brighter radio emission in 2021 to be due to a
higher particle density in the synchrotron emitting plasma.

Our study suggests that rapid multi-frequency radio follow-up observations of
recurrent novae are crucial to tie down various absorption processes and infer
the properties of shock and CBM.

## Acknowledgements

We thank the anonymous referee for the comments that helped to improve our
manuscript. We thank the staff of the GMRT that made these observations
possible. The GMRT is run by the National Centre for Radio Astrophysics of the
Tata Institute of Fundamental Research. Nayana A.J. would like to acknowledge
DST-INSPIRE Faculty Fellowship (IFA20-PH-259) for supporting this research.
KPS thanks the Indian National Science Academy for support under the INSA
Senior Scientist Programme.

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Acciari et al. (2022) Acciari V. A., et al., 2022, Nature Astronomy, 6, 689
  * Anupama & Mikołajewska (1999) Anupama G. C., Mikołajewska J., 1999, A&A, 344, 177
  * Bailer-Jones et al. (2021) Bailer-Jones C. A. L., Rybizki J., Fouesneau M., Demleitner M., Andrae R., 2021, AJ, 161, 147
  * Bode et al. (2007) Bode M. F., Harman D. J., O’Brien T. J., Bond H. E., Starrfield S., Darnley M. J., Evans A., Eyres S. P. S., 2007, ApJ, 665, L63
  * Booth et al. (2016) Booth R. A., Mohamed S., Podsiadlowski P., 2016, MNRAS, 457, 822
  * Cheung et al. (2022) Cheung C. C., et al., 2022, ApJ, 935, 44
  * Chevalier & Fransson (2006) Chevalier R. A., Fransson C., 2006, ApJ, 651, 381
  * Chomiuk et al. (2012) Chomiuk L., et al., 2012, ApJ, 761, 173
  * Chomiuk et al. (2021) Chomiuk L., et al., 2021, ApJS, 257, 49
  * Chugai (2010) Chugai N. N., 2010, Astronomy Letters, 36, 540
  * Das et al. (2006) Das R., Banerjee D. P. K., Ashok N. M., 2006, ApJ, 653, L141
  * De Sarkar et al. (2023) De Sarkar A., Nayana A. J., Roy N., Razzaque S., Anupama G. C., 2023, ApJ, 951, 62
  * Diesing et al. (2023) Diesing R., Metzger B. D., Aydi E., Chomiuk L., Vurm I., Gupta S., Caprioli D., 2023, ApJ, 947, 70
  * Drake et al. (2009) Drake J. J., et al., 2009, ApJ, 691, 418
  * Eyres et al. (2009) Eyres S. P. S., et al., 2009, MNRAS, 395, 1533
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306
  * Greisen (2003) Greisen E. W., 2003, in Heck A., ed., Astrophysics and Space Science Library Vol. 285, Information Handling in Astronomy - Historical Vistas. p. 109, doi:10.1007/0-306-48080-8_7
  * H. E. S. S. Collaboration et al. (2022) H. E. S. S. Collaboration et al., 2022, Science, 376, 77
  * Hjellming et al. (1986) Hjellming R. M., van Gorkom J. H., Taylor A. R., Sequist E. R., Padin S., Davis R. J., Bode M. F., 1986, ApJ, 305, L71
  * Kantharia et al. (2007) Kantharia N. G., Anupama G. C., Prabhu T. P., Ramya S., Bode M. F., Eyres S. P. S., O’Brien T. J., 2007, ApJ, 667, L171
  * Longair (2011) Longair M. S., 2011, High Energy Astrophysics
  * McDonald & van Loon (2007) McDonald I., van Loon J. T., 2007, A&A, 476, 1261
  * Mikolajewska et al. (2021) Mikolajewska J., Aydi E., Buckley D., Galan C., Orio M., 2021, The Astronomer’s Telegram, 14852, 1
  * Mondal et al. (2018) Mondal A., Anupama G. C., Kamath U. S., Das R., Selvakumar G., Mondal S., 2018, MNRAS, 474, 4211
  * Munari & Valisa (2021a) Munari U., Valisa P., 2021a, arXiv e-prints, p. arXiv:2109.01101
  * Munari & Valisa (2021b) Munari U., Valisa P., 2021b, The Astronomer’s Telegram, 14840, 1
  * Munari & Valisa (2022) Munari U., Valisa P., 2022, arXiv e-prints, p. arXiv:2203.01378
  * Munari et al. (2022) Munari U., Giroletti M., Marcote B., O’Brien T. J., Veres P., Yang J., Williams D. R. A., Woudt P., 2022, A&A, 666, L6
  * Narumi et al. (2006) Narumi H., Hirosawa K., Kanai K., Renz W., Pereira A., Nakano S., Nakamura Y., Pojmanski G., 2006, IAU Circ., 8671, 2
  * Nayana et al. (2017) Nayana A. J., et al., 2017, MNRAS, 467, 155
  * Nayana et al. (2021) Nayana A. J., Anupama G. C., Banerjee D., Roy N., Singh K. P., Sonith L. S., 2021, The Astronomer’s Telegram, 14899, 1
  * Ness et al. (2023) Ness J. U., et al., 2023, A&A, 670, A131
  * O’Brien et al. (2006) O’Brien T. J., et al., 2006, Nature, 442, 279
  * O’Brien et al. (2008) O’Brien T. J., et al., 2008, in Evans A., Bode M. F., O’Brien T. J., Darnley M. J., eds, Astronomical Society of the Pacific Conference Series Vol. 401, RS Ophiuchi (2006) and the Recurrent Nova Phenomenon. p. 239
  * Origlia et al. (2007) Origlia L., Rood R. T., Fabbri S., Ferraro F. R., Fusi Pecci F., Rich R. M., 2007, ApJ, 667, L85
  * Orio (1993) Orio M., 1993, A&A, 274, L41
  * Orio et al. (2022) Orio M., et al., 2022, ApJ, 938, 34
  * Orlando et al. (2009) Orlando S., Drake J. J., Laming J. M., 2009, A&A, 493, 1049
  * Osborne et al. (2011) Osborne J. P., et al., 2011, ApJ, 727, 124
  * Page et al. (2022) Page K. L., et al., 2022, MNRAS, 514, 1557
  * Pandey et al. (2022) Pandey R., Habtie G. R., Bandyopadhyay R., Das R., Teyssier F., Guarro Fló J., 2022, MNRAS, 515, 4655
  * Prialnik et al. (1982) Prialnik D., Livio M., Shaviv G., Kovetz A., 1982, ApJ, 257, 312
  * Rupen et al. (2008) Rupen M. P., Mioduszewski A. J., Sokoloski J. L., 2008, ApJ, 688, 559
  * Schaefer (2010) Schaefer B. E., 2010, ApJS, 187, 275
  * Sokoloski et al. (2008) Sokoloski J. L., Rupen M. P., Mioduszewski A. J., 2008, ApJ, 685, L137
  * Sokolovsky et al. (2021) Sokolovsky K., et al., 2021, The Astronomer’s Telegram, 14886, 1
  * Walder et al. (2008) Walder R., Folini D., Shore S. N., 2008, A&A, 484, L9
  * Weiler et al. (1986) Weiler K. W., Sramek R. A., Panagia N., van der Hulst J. M., Salvati M., 1986, ApJ, 301, 790
  * Weiler et al. (2002) Weiler K. W., Panagia N., Montes M. J., Sramek R. A., 2002, ARA&A, 40, 387
  * Williams et al. (2021) Williams D., O’Brien T., Woudt P., Nyamai M., Green D., Titterington D., Fender R., Sivakoff G., 2021, The Astronomer’s Telegram, 14849, 1
  * Zheng et al. (2022) Zheng J.-H., Huang Y.-Y., Zhang Z.-L., Zhang H.-M., Liu R.-Y., Wang X.-Y., 2022, Phys. Rev. D, 106, 103011
  * de Ruiter et al. (2023) de Ruiter I., Nyamai M. M., Rowlinson A., Wijers R. A. M. J., O’Brien T. J., Williams D. R. A., Woudt P., 2023, MNRAS,
  * van Loon et al. (2006) van Loon J. T., McDonald I., Oliveira J. M., Evans A., Boyer M. L., Gehrz R. D., Polomski E., Woodward C. E., 2006, A&A, 450, 339


# Optimization for Secure and Humane Border Operations

Fatemeh Farajzadeh  Data Science Program, Worcester Polytechnic Institute,
Worcester, MA 01609, USA  Andrew C. Trapp  Data Science Program, Worcester
Polytechnic Institute, Worcester, MA 01609, USA  WPI Business School,
Worcester Polytechnic Institute, Worcester, MA 01609, USA

In 2020 tens of millions of migrants worldwide sought safety due to various
push-pull factors such as economic, geographic, demographic, and societal
challenges [1, 2]. Large movements of human flows toward national borders
constitute a complex humanitarian crisis that requires appropriate
preparations of border resources. Motivated by prevailing national agenda
issues for improved border support, we identify common operational challenges
along the border and discuss how they can be addressed using advanced
analytics. Multiple stakeholders, corresponding goals and complex decisions
call for advanced optimization techniques to holistically treat such
challenges. This chapter highlights how optimization can be employed to
improve both the security and humanitarian aspects of border operations in the
context of migration flows.

## Introduction

Shifting and evolving migration patterns call for border plans that feature
proactive mechanisms to improve operations in more effectively using limited
resources. The European Border and Coast Guard Agency, also known as Frontex,
reported the need to deploy additional resources and increase operations in
the Mediterranean for three primary efforts [3]. The US Customs and Border
Protection (CBP) agency has reported that movements of large migrant groups
have resulted in excessive use of border capacities and sacrificing the border
security mission [4]. While border operations must necessarily consider
security and orderliness, a large majority of countries worldwide (over 150)
are signatories on the 1951 Convention Relating to the Status of Refugees [5],
as well as its 1967 Protocol [6]. Thus, the proper management of such
operational resources involves dual aspects of security and humanity, and
offers opportunities for advanced analytics to inform critical border
decisions.

Researchers in the early 2000’s primarily emphasized issues pertaining to
homeland security operations [7]. Emergency preparedness and response, border
and transportation security, threat analysis, and assessment were collectively
positioned as the opportunities with the greatest potential for improvement
via optimization techniques. The use of discrete optimization to improve
homeland security and disaster management has also been discussed [8]. This
study classified emergency and security operations into four pre-defined
components of a disaster life cycle, including mitigation, preparedness,
response, and recovery. Moreover, basic optimization models were presented
according to the disaster life cycle sequence: checkpoint screening of
arrivals and network fortification, prepositioning of supplies and evacuation
planning, rescue and aid distribution, and infrastructure restoration. These
efforts support at-risk migrants through securing EU borders by the deployment
of high-tech surveillance systems, disrupting human smuggling and trafficking
networks by effective patrolling strategies, and humanitarian relief and
processing [9, 10]. With the fluctuating need to accommodate waves of
displaced peoples, additional facilities are needed along borders where the
timely supply and delivery of essential products and services are key
decisions [11].

The management of human flows at national borders constitutes a complex
emergency that requires a time- and cost-efficient, fair and secure planning
approach. The planning characteristics of complex emergencies differ from
natural disaster response operations that typically are sudden-onset relief
efforts and involve heightened vulnerabilities that, with time, tend to
subside. Rather, border operations in relation to migration flow management is
an ongoing process, and it takes place in an ad hoc manner across multiple
agencies. Real-time planning decisions are affected by the evolution of
stochastic flows over time, including high levels of uncertainty on human
needs, variability of supply cycle duration, and availability of accurate
information regarding the engagement of abusive groups such as smugglers and
human traffickers. To date, relatively few studies have leveraged analytical-
based methods to improve border operations in response to shifting human
flows.

The aforementioned considerations motivate new structural components and
modeling choices for the design of appropriate mathematical models. Fairness
considerations such as deprivation costs must be balanced alongside
traditional objectives like cost, time, and distance [12]. Other important
modeling components include the presence of multiple decision-makers that
influence any one party’s control level in making unilateral decisions, as
well as a high level of uncertainty and unreliability associated with planning
parameters. Mathematical optimization appears to be the most widely used
methodology in managing relief resources while addressing the aforementioned
challenges under a dynamic planning environment [13]. And, it can be adapted
to address the problems arising from global immigration crises with the goal
of properly allocating security and support resources. This chapter proceeds
to discuss border operations in two areas, the staging of security and
shelter-based resources through location analysis; and search and rescue
challenges in borderlands both for saving migrants and protecting against bad
actors.

## Location and Allocation Challenges in Border Operations

Decisions on locating facilities have perhaps the greatest impact on the
performance of relief operations in humanitarian crises, due to the direct
relationship with response times and distribution costs. Key operational
issues have been studied by researchers and linked to the broader operations
research literature by referring to well-known location decision problems in
collaboration with the Turkish Red Crescent [14]. A recent study thoroughly
summarized various types of facility location challenges according to aspects
such as the location space (continuous, network-based, discrete), the type of
facilities (suppliers, distribution centers, shelters), the critical decision-
makings criteria (equity, cost, reliability), and source of uncertainty
(supply, demand, network) [15]. Although facility location problems have
theoretically been addressed in isolation, in practice, location decisions and
positioning challenges are coupled with allocation, inventory, distribution,
and routing decisions in disaster relief operation settings [15]. As it
relates to accommodating and processing displaced populations at international
borders, scholars must address the same positioning challenges while
accounting for the specific aspects of the associated intricacies and
corresponding planning complexities.

The major services of border control agencies, whenever human flows arrive,
are administrative and essential aid arrangements that require three primary
types of decisions. These decisions include the location of service bases in
harmony with human flow directions and volumes, capacity and inventory levels
according to functions of the respective facilities, and timely service and
aid distribution. The aim of establishing borderland service bases is to focus
on essential operations: safety and security, health and sanitation, temporary
settlements, and processing activities for identification and registration.
The aforementioned operations will provide decision-making authorities with
the basic steps of strategic and tactical actions when responding to the
humanitarian crises at the borders.

A limited number of studies have taken a similar optimization perspective in
refugee camp management and administration, seeking a cost-effective and
humane planning approach. To optimally allocate aid to refugee camps under
demand and replenishment cycle uncertainty, an integer optimization model with
the objective of minimizing the expected total cost of holding excess
inventory, referring external refugees elsewhere, and depriving internal
refugees of critical aid is proposed [16]. Another study investigated the
distribution of cash and e-vouchers to refugees as well as the routing problem
of trucks providing child-friendly spaces to vulnerable refugee children. The
authors presented a bi-objective model based on covering vehicle routing
problems with integrated tours where the aim is to minimize the unsatisfied
demand and traveled distance [17]. A multi-objective mixed-integer program was
developed to optimally choose the locations and capacity level of short-term
refugee camps in Syria to minimize maintenance and startup costs and ensure
camp safety [18]. A combination of machine learning and integer optimization
techniques are used to assign refugees to initial locations with the aim of
improving integration for both refugees and destinations [19].

While border operations involve issues related to flow management, they also
involve security operations to protect the border from bad actors. Not
surprisingly, these two missions are intertwined. At the same time, in the era
of emerging technologies such as surveillance towers, underground sensing
technologies, drones, and other AI-based solutions, border control agencies
are increasingly transforming their operations through investing in the
deployment of smart tools. From the security perspective, the use of
technology enhances situational awareness, agent effectiveness, and safety by
supporting border patrol agents to cover more territory for the purposes of
interdiction. From the humanitarian perspective, these technologies help
agents efficiently discriminate items of interest, protect vulnerable people
from criminal networks, and proactively respond to detected events. Location
optimization models can improve the process of designing an efficient
surveillance system to serve both purposes. This is a relatively new area of
research that has recently attracted the attention of operations researchers.
A maximum coverage optimization model was developed to determine hubs for
unmanned aerial vehicles via two subproblems, namely, the selection of hub
locations and demand point allocation [20]. Similarly, a multi-objective bi-
level optimization model was developed to address the problem of sensor
relocation for the purpose of maximizing the minimal expected exposure of
intruders when traversing a defended border region [21]. Allocation
optimization problems over a rolling horizon are also of great importance in
this highly variable planning setting. These sets of optimization problems
fall under dynamic resource allocation problems that are often modeled as
dynamic assignment and stochastic knapsack problems [22].

## Search and Rescue Challenges in Border Operations

Official reports have documented the risks of illegal border crossings in
remote areas arising from a lack of processing at lawful ports of entry [23,
24]. Migrants on the move face substandard living conditions and many types of
risks, including getting lost, dangerous terrain, natural hazards, and man-
made dangers such as trafficking, fraud, and smuggling. Border agents are
responsible for patrolling the assigned borderlands searching for any
potential hazards, notably, people in need of humanitarian support. Search and
rescue practices must be agile and able to track autonomous behaviors and
movements in a situation the clock is ticking. These practices need skilled
individuals with different capabilities to operate as a multi-talented team.
In addition, a disaster-impacted zone with a time-dependent demand rate urges
rescuers to be flexible in dealing with changing circumstances. There exist
specialized units that are trained primarily to assist injured or stranded
migrants, such as Border Patrol Search, Trauma, and Rescue Unit (BORSTAR) in
the United States and and Search and Rescue (SAR) in Europe [25, 26]. Thus,
the use of optimization for team planning and scheduling of various rescue
missions is of particular importance. Although rescue operations have already
risen to prominence in disaster emergency management, a significant gap exists
to address such problems using optimization techniques.

Search and rescue strategies must adapt to newly revealed information of team
efforts, changing circumstances, and the progress of ongoing rescue operation.
Beyond sharing prompt information of ongoing operations updated by central
planners, real-life constraints must be reflected in rescue teams tours
planning, including hierarchical skill levels, collaborative rescue service,
time-dependent service duration, and idling rest periods [27]. In addition to
operational constraints considered in rescue operations of natural disasters,
border operations management is characterized by dynamics associated with
various uncertain service requests or unforeseen events that typically took
place in close relation with border security operations. Such problem dynamics
may be tackled through myopic approaches, by look-ahead procedures that rely
only on available information at the time of decision-making or perhaps take
into account probabilistic information concerning future events.

Search and rescue tasks can be classified by the type of terrain, such as air-
sea rescue, urban rescue, hard-to-access rescue such as mountainous terrain,
and general ground rescue. Borderland search and rescue may feature additional
unique characteristics, such as incident type, level of complexity,
geographical factors, and immediate survival resources, all in a traumatic
situation. Respecting geographical factors, the decision of rescue team
location in sea rescue operations differs based upon the given required
operation capacity and boat workload [28]. This study developed an incident-
based allocation optimization model to assign responsibility zones of the
rescue team considering historical demand locations to account for demand
uncertainty. A relocation possibilities of resources, using a scenario-based
multi-objective programming model to generate a balanced workload distribution
in the scheme of location-allocation analysis also has been addressed [29].
Another study used a two-stage stochastic programming framework to fairly
cover the demand of affected zones for the purpose of maximizing the utility
of deployed rescue teams while also considering the transfer of resources to
other affected areas or discharging decisions so as minimize the completion
times of search and rescue teams [30]. The challenges of determining an
optimal strategy for dispatching rescue teams to post-disaster sites was also
another subject of interest. Relatively, a two-stage stochastic programming
formulated to first seek a set of tours while maximizing the total expected
rescued people and subsequently, second-stage variables update the allocation
decisions over time by the accessibility of newly released resources and
revealed planning information [31]. To enhance the performance of rescue
operations, multimodal interaction is presented as an approach to incorporate
a mixed-use. A goal programming approach is proposed to dynamically optimize
the location and allocation of rescuers, boats, and helicopters considering
the seasonality of demand [32]. The author attempted to achieve a balanced
workload for all boats and helicopters such that the sum of shortage and
excess operation hours are minimized via deviations from the critical
incidents response time and budget limitation. An integrated platform of
ground and aerial vehicles is also introduced to support search and rescue
activities [33].
The widespread implementation of cutting-edge technologies for surveillance
and monitoring makes it possible to quickly assess key planning parameters,
follow changes in these conditions as they evolve, and identify at-risk
persons that may need rescuing. A variety of search techniques exist to
effectively cover the target area. These algorithms include Bayesian
approaches for tracking probable target locations, Markov decision processes
for unobserved underlying states, and greedy heuristics where minimizing the
solving time is the optimization criterion [34]. A study on creating a good
drone presented various successes and practical factors adapted to more
effectively detect migrants, which can help researchers familiarize themselves
with the operational challenges of drones [35]. A study explored simulation-
optimization techniques to address the potential use of drones in searching
for and locating victims of search and rescue providers in a mountain
environment [36].
We believe this chapter outlines key borderland operations challenges in
location, allocation, and search and rescue. Guided by dual humanitarian and
security purposes, in conjunction with increasing technological capabilities,
we sincerely hope operations research and analytics researchers consider the
excellent opportunity to further explore this complex planning environment and
underlying optimization challenges.

## Acknowledgements

We are grateful for the assistance of many in making this chapter successful.
In particular, the authors are grateful for the support of the National
Science Foundation (Operations Engineering) grant CMMI-1825348.

## References

  * UNHCR [2021] UNHCR.  Global Trends: Forced Displacement in 2020.  Technical report, United Nations High Commissioner for Refugees, 2021\.  URL https://www.unhcr.org/en-us/60b638e37.pdf.
  * IOM, UN Migration [2021] IOM, UN Migration.  World migration report.  Technical report, 2021.  Accessed: 2021-12-01.
  * Council of the European Union [2020] Council of the European Union.  Saving lives at sea and targeting criminal networks.  https://www.consilium.europa.eu/en/policies/eu-migration-policy/saving-lives-at-sea/#, 2020\.
  * Customs and Border Protection [2020] Customs and Border Protection.  Border crisis: CBP’s response.  www.cbp.gov/frontline/border-crisis-cbp-s-response, 2020.
  * Assembly [1951] UN General Assembly.  Convention relating to the status of refugees.  Technical report, 1951.  Accessed: 2021-09-22.
  * UNHCR [2011] UNHCR.  The 1951 convention relating to the status of refugees and its 1967 protocol.  https://www.unhcr.org/en-us/about-us/background/4ec262df9/1951-convention-relating-status-refugees-its-1967-protocol.html, 2011\.  Accessed: 2021-09-22.
  * Wright et al. [2006] P. Daniel Wright, Matthew J. Liberatore, and Robert L. Nydick.  A survey of operations research models and applications in homeland security.  _INFORMS Journal on Applied Analytics_ , 36(6):514–529, 2006.
  * Albert McLay [2015] Laura Albert McLay.  Discrete optimization models for homeland security and disaster management.  _INFORMS TutORials in Operations Research_ , pages 111–132, 2015\.
  * Walters [2010] William Walters.  Foucault and Frontiers: Notes on the birth of the humanitarian border.  In _Governmentality_ , page 27. Taylor & Francis, 2010.
  * Csernatoni [2018] Raluca Csernatoni.  Constructing the EU’s high-tech borders: FRONTEX and dual-use drones for border management.  _European Security_ , 27(2):175–200, 2018.
  * Customs and Border Protection [2019] Customs and Border Protection.  New temporary facilities in Texas and Arizona expand CBP holding capacity.  www.cbp.gov/newsroom/spotlights/new-temporary-facilities-expand-cbp-holding-capacity, 2019.
  * Holguín-Veras et al. [2013] José Holguín-Veras, Noel Pérez, Miguel Jaller, Luk N. Van Wassenhove, and Felipe Aros-Vera.  On the appropriate objective function for post-disaster humanitarian logistics models.  _Journal of Operations Management_ , 31(5):262–280, 2013.
  * Galindo and Batta [2013] Gina Galindo and Rajan Batta.  Review of recent developments in OR/MS research in disaster operations management.  _European Journal of Operational Research_ , 230(2):201–211, 2013.
  * Karsu et al. [2019] Ozlem Karsu, Bahar Y. Kara, and Bayram Selvi.  The refugee camp management: a general framework and a unifying decision-making model.  _Journal of Humanitarian Logistics and Supply Chain Management_ , 9:131–150, 2019.
  * Dönmez et al. [2021] Zehranaz Dönmez, Bahar Y. Kara, Özlem Karsu, and Francisco Saldanha-da Gama.  Humanitarian facility location under uncertainty: Critical review and future prospects.  _Omega_ , 102:102–393, 2021.
  * Azizi et al. [2021] Shima Azizi, Cem Deniz Caglar Bozkir, Andrew C. Trapp, O. Erhun Kundakcioğlu, and Ali Kaan Kurbanzade.  Aid allocation for camp-based and urban refugees with uncertain demand and replenishments.  _Production and Operations Management_ , 0:1–20, 2021.
  * Buluç, Elfe Naz [2018] Buluç, Elfe Naz.  Covering vehicle routing problem: applications for refugee related services, 2018.  Master thesis, Bilkent University.
  * Vatasoiu et al. [2015] Sorin Vatasoiu, Martin Carlsen, and Julia Romanski.  Optimizing the size and location of short-term refugee camps in Syria.  Technical report, 2015.
  * Ahani et al. [2021] Narges Ahani, Tommy Andersson, Alessandro Martinello, Alexander Teytelboym, and Andrew C. Trapp.  Placement optimization in refugee resettlement.  _Operations Research_ , 69:131–150, 2021.
  * Yildirim [2011] Mehmet Fatih Yildirim.  Optimization of hub locations of UAVs for border monitoring as coverage problem.  2011\.  Master thesis, State University of New York at Binghamton.
  * Lessin et al. [2019] Aaron M. Lessin, Brian J. Lunday, and Raymond R. Hill.  A multi-objective, bilevel sensor relocation problem for border security.  _IISE Transactions_ , 51(10):1091–1109, 2019\.
  * Kleywegt and Papastavrou [1998] Anton J. Kleywegt and Jason D. Papastavrou.  The dynamic and stochastic knapsack problem.  _Operations Research_ , 46:17–35, 1998.
  * Custom and Border Protection [2019] Custom and Border Protection.  Border patrol warns immigrants of life-threatening risks of crossing Rio Grande river.  https://www.cbp.gov/newsroom/local-media-release/border-patrol-warns-immigrants-life-threatening-risks-crossing-rio, 2019\.
  * UNHCR [2017] UNHCR.  Refugees and migrants face heightened risks while trying to reach Europe.  https://www.unhcr.org/en-us/news/press/2017/2/58b458654/refugees-migrants-face-heightened-risks-trying-reach-europe-unhcr-report.html, 2017\.
  * Customs and Border Protection [2014] Customs and Border Protection.  Border Patrol Search, Trauma, and Rescue (BORSTAR).  https://www.cbp.gov/document/fact-sheets/border-patrol-search-trauma-and-rescue-borstar, 2014\.
  * FRONTEX [2016] FRONTEX.  Search and Rescue (SAR).  https://frontex.europa.eu/we-support/types-of-operations/search-rescue/, 2016\.
  * [27] Felix Wex, Guido Schryen, Stefan Feuerriegel, and Dirk Neumann.  Emergency response in natural disaster management: Allocation and scheduling of rescue units.  _European Journal of Operational Research_ , 235(3):697–708.
  * Razi and Karatas [2016] Nasuh Razi and Mumtaz Karatas.  A multi-objective model for locating search and rescue boats.  _European Journal of Operational Research_ , 254:279–293, 2016.
  * Karatas [2021a] Mumtaz Karatas.  A dynamic multi-objective location-allocation model for search and rescue assets.  _European Journal of Operational Research_ , 288:620–633, 2021a.
  * Ahmadi et al. [2020] Ghazaleh Ahmadi, Reza Tavakkoli‑Moghaddam, Armand Baboli, and Mehdi Najafi.  A decision support model for robust allocation and routing of search and rescue resources after earthquake: a case study.  _Operational Research_ , 2020.  doi: https://doi.org/10.1007/s12351-020-00591-5.
  * Chen and Miller-Hooks [2012] Lichun Chen and Elise Miller-Hooks.  Optimal team deployment in urban search and rescue.  _Transportation Research Part B: Methodological_ , 46:984–999, 2012.
  * Karatas [2021b] Mumtaz Karatas.  A dynamic multi-objective location-allocation model for search and rescue assets.  _European Journal of Operational Research_ , 288(2):620–633, 2021b.
  * Cacace et al. [2016] Jonathan Cacace, Alberto Finzi, and Vincenzo Lippiello.  Multimodal interaction with multiple co-located drones in search and rescue missions.  _Computing Research Repository_ , abs/1605.07316, 2016.  URL :http://arxiv.org/abs/1605.07316.
  * Waharte and Trigoni [2010] Sonia Waharte and Niki Trigoni.  Supporting search and rescue operations with UAVs.  In _2010 International Conference on Emerging Security Technologies_ , pages 142–147, 2010.
  * Maria Gabrielsen Jumbert [2016] Maria Gabrielsen Jumbert.  Creating the EU drone: control, sorting, and search and rescue at sea.  In _The Good Drone_ , page 212. Routledge, 2016.
  * Karaca et al. [2018] Yunus Karaca, Mustafa Cicek, Ozgur Tatli, Aynur Sahin, Sinan Pasli, Muhammed Fatih Beser, and Suleyman Turedi.  The potential use of unmanned aircraft systems (drones) in mountain search and rescue operations.  _The American Journal of Emergency Medicine_ , 36(4):583–588, 2018.