# Fermi LAT Gamma-ray Detection of the Recurrent Nova RS Ophiuchi during its 2021 Outburst

C. C. Cheung, T. J. Johnson, P. Jean, M. Kerr, K. L. Page, J. P. Osborne, A. P. Beardmore, K. V. Sokolovsky, F. Teyssier, S. Ciprini, et al.

► To cite this version:

C. C. Cheung, T. J. Johnson, P. Jean, M. Kerr, K. L. Page, et al., Fermi LAT Gamma-ray Detection of the Recurrent Nova RS Ophiuchi during its 2021 Outburst. The Astrophysical Journal, 2022, 935, 10.3847/1538-4357/ac7eb7. insu-03863537

HAL Id: insu-03863537

https://insu.hal.science/insu-03863537

Submitted on 21 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Fermi LAT Gamma-ray Detection of the Recurrent Nova RS Ophiuchi during its 2021 Outburst

C. C. Cheung¹, P. Jean³,⁴, M. Kerr¹, K. L. Page⁵, J. P. Osborne⁵, A. P. Beardmore⁵, K. V. Sokolovsky⁶,⁷, F. Teyssier⁸, S. Ciprinı⁹,¹⁰, G. Martí-Devesa¹¹, I. Mereu¹²,¹³, S. Razzaque¹⁴, K. S. Wood¹⁵, S. N. Shore¹⁶, S. Korotkiy¹⁷, A. Levina¹⁸, and A. Blumenzweig¹⁸

¹Space Science Division, Naval Research Laboratory, Washington, DC 20375, USA; Teddy.Cheung@nrl.navy.mil

²College of Science, George Mason University, Fairfax, VA 22030, resident at Naval Research Laboratory, Washington, DC 20375, USA
Tyrel.Johnson.ctr@nrl.navy.mil

³CNRS, IRAP, F-31028 Toulouse cedex 4, France; Pierre.Jean@irap.omp.eu

⁴GAHEC, Université de Toulouse, UPS-OMP, IRAP, Toulouse, France

⁵School of Physics and Astronomy, University of Leicester, Leicester LE1 7RH, UK

⁶Center for Data Intensive and Time Domain Astronomy, Department of Physics and Astronomy, Michigan State University, 567 Wilson Rd, East Lansing, MI 48824, USA

⁷Sternberg Astronomical Institute, Moscow State University, Universitetskii pr. 13, 119992 Moscow, Russia

⁸Astronomical Ring for Amateur Spectroscopy Group, F-76100 Rouen, France

⁹Istituto Nazionale di Fisica Nucleare, Sezione di Roma “Tor Vergata”, I-00133 Roma, Italy

¹⁰Space Science Data Center—Agenzia Spaziale Italiana, Via del Politecnico, snc, I-00133, Roma, Italy

¹¹Institut für Astro- und Teilchenphysik, Leopold-Franzens-Universität Innsbruck, A-6020 Innsbruck, Austria

¹²Dipartimento di Fisica, Università degli Studi di Perugia, I-06123 Perugia, Italy

¹³Istituto Nazionale di Fisica Nucleare, Sezione di Perugia, I-06123 Perugia, Italy

¹⁴Centre for Astro-Particle Physics (CAPP) and Department of Physics, University of Johannesburg, PO Box 524, Auckland Park 2006, South Africa

¹⁵Praxis Inc., Alexandria, VA 22303, resident at Naval Research Laboratory, Washington, DC 20375, USA

¹⁶Dipartimento di Fisica “Enrico Fermi”, Università di Pisa, and INFN—Sezione di Pisa, largo B. Pontecorvo 3, I-56127 Pisa, Italy

¹⁷Astrology, Ka-Dar astronomy foundation, Kuzminki, P.O. Box 82, 142717 Moscow region, Russia

¹⁸Israeli Astronomical Association, Tel Aviv, Israel

Received 2022 April 28; revised 2022 June 22; accepted 2022 July 3; published 2022 August 12

## Abstract

We report the Fermi LAT γ-ray detection of the 2021 outburst of the symbiotic recurrent nova RS Ophiuchi. In this system, unlike classical novae from cataclysmic binaries, the ejecta from the white dwarf form shocks when interacting with the dense circumstellar wind environment of the red giant companion. We find the LAT spectra from 50 MeV to ~20–23 GeV, the highest-energy photons detected in some subintervals, are consistent with π⁰-decay emission from shocks in the ejecta as proposed by Tatischeff & Hernanz for its previous 2006 outburst. The LAT light curve displayed a fast rise to its peak >0.1 GeV flux of ~6 × 10⁻⁶ ph cm⁻² s⁻¹ beginning on day 0.745 after its optically constrained eruption epoch of 2021 August 8.50. The peak lasted for ~1 day and exhibited a power-law decline up to the final LAT detection on day 45. We analyze the data on shorter timescales at early times and found evidence of an approximate doubling of emission over ~200 minutes at day 2.2, possibly indicating a localized shock-acceleration event. Comparing the data collected by the American Association of Variable Star Observers, we measured a constant ratio of ~2.8 × 10⁻³ between the γ-ray and optical luminosities except for a ~5×smaller ratio within the first day of the eruption likely indicating attenuation of γ rays by ejecta material and lower high-energy proton fluxes at the earliest stages of the shock development. The hard X-ray emission due to bremsstrahlung from shock-heated gas traced by the Swift-XRT 2–10 keV light curve peaked at day ~6, later than at GeV and optical energies. Using X-ray derived temperatures to constrain the velocity profile, we find the hadronic model reproduces the observed >0.1 GeV light curve.

Unified Astronomy Thesaurus concepts: Gamma-ray transient sources (1853); Recurrent novae (1366)

Supporting material: data behind figure

## 1. Introduction

RS Ophiuchi (RS Oph) is one of the best-studied recurrent novae because of its numerous outbursts since the first detection in 1898 (Pickering 1905). It has recurred at irregular intervals (with inferred outbursts missed due to solar conjunction; Oppenheimer & Mattei 1993; Schaefer 2004) of 9 to 27 yr up to its previous outburst on 2006 February 12 (Narumi et al. 2006; Evans et al. 2008; Schaefer 2010). RS Oph is a symbiotic binary system with

a 453.6 ± 0.4 day orbital period consisting of a massive white dwarf (1.2–1.4 $M_☉$) and a red giant (RG) commonly identified as type M0 III (Dobrzycka & Kenyon 1994; Barry et al. 2008; Brandi et al. 2009). Its widely adopted distance of 1.6 ± 0.3 kpc (Hjellming et al. 1986; Bode 1987) is assumed here to facilitate direct comparison to previous work; this distance is consistent with the value of 1.4$^{+0.6}_{-0.2}$ kpc obtained considering various methods (Barry et al. 2008)—but see arguments for greater distances presented by Schaefer (2009, 2018), Montez et al. (2022), and MAGIC Collaboration et al. (2022).

Its 2006 outburst was well studied (Evans et al. 2008) with observations of hard X-ray emission from 2 to 25 keV (Sokoloski et al. 2006) and 14–50 keV (Bode et al. 2006), which indicated

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
shocked emission in the nova ejecta, and high-resolution radio imaging resolved the shocked regions (O'Brien et al. 2006; Rupen et al. 2008; Sokoloski et al. 2008). In γ rays, Tatischeff & Hernanz (2007) predicted high-energy particle acceleration in the nova ejecta from interactions with the dense RG wind that could have been observed in 2006 at GeV energies,<fnref n="19" /> but that explosion preceded the launch of Fermi in 2008. Instead, the first >0.1 GeV detection of a nova by the Fermi Large Area Telescope (LAT; Atwood et al. 2009) in 2010 was of a less-known symbiotic binary V407 Cyg (Abdo et al. 2010; Cheung et al. 2010). This discovery demonstrated the viability of the nova explosion in the RG wind model and helped to solidify the predictions of GeV emission in RS Oph in particular (Hernanz & Tatischeff 2012).

Thus the next outburst of RS Oph was highly anticipated by multiwavelength observers, particularly its >0.1 GeV observation by the Fermi LAT. Indeed, a new optical outburst from RS Oph was discovered in 2021, independently by A. Amorim (August 8.913), E. Muyllaert (August 8.920), and K. Geary (August 8.931); see Beck (2021). The time since its 2006 detection is 15.5 yr, close to its average recurrence interval of 14.7 yr (Schaefer 2010). Following the Geary (2021) announcement, we reported the independent Fermi LAT detection of the nova while performing its normal all-sky monitoring during the last 6 hr interval of 2021 August 8, overlapping with the optical discovery epoch (Cheung et al. 2021a).<fnref n="20" />

Here, we present the details of the Fermi LAT observations of RS Oph 2021. In the following, Section 2 describes the LAT observations and analysis. Section 3 describes the LAT results, and comparison to optical and X-ray data. Section 4 presents the results of the LAT spectral variability analysis for four defined emission phases and examines the originally proposed π⁰-decay model (Tatischeff & Hernanz 2007; Hernanz & Tatischeff 2012) to reproduce the γ-ray spectra and light curves. The results are summarized in Section 5. All times are UTC, while days relative to the optical eruption epoch, t₀ = 2021 August 8.5 (Munari & Valisa 2021, and Section B.1) are used to describe the RS Oph outburst.

## 2. Fermi LAT Observations and Data Analysis

For the LAT analysis, we used Pass 8 (P8R3 SOURCE class data<fnref n="21" /> (Atwood et al. 2013; Bruel et al. 2018), as defined under the P8R3_SOURCE_V3 instrument response functions. Photons with 0.05–300 GeV energies, within 25° of R.A., decl. (J2000) = 267.7, −6.7, and with a maximum zenith angle of 90° were selected. The center of the region of interest (ROI) was chosen to put RS Oph near the center but not on a pixel edge in the binned analysis described below. We filtered the events to include only good time intervals when the LAT data were flagged as good, and the instrument was in nominal science observations mode. All data processing and analysis were done using version 2.0.8 of fermitools (Fermi Science Support Development Team 2019).

We constructed a spatial and spectral model of our ROI starting from the 4FGL catalog incremental data release 3

(DR3; Abdollahi et al. 2020, 2022) based on 12 years of LAT data. We included all DR3 sources within 35° of the ROI center as well as a model for the Galactic diffuse emission (gll_iem_v07.fits) and a diffuse isotropic emission component (iso_P8R3_SOURCE_V3_v1.txt).<fnref n="22" /> Note the source of emission from RS Oph as seen in the counts map (Figure 1) is relatively isolated from other point sources (the closest 4FGL-DR3 source is 4FGL J1752.4–0758, offset by 1.38) and the Galactic diffuse emission, but is close enough to the latter to potentially cause low-level contamination at the nova position, especially at the lowest energies.

We added a point source to our model at the optical position of the nova (R.A., decl. = 267.5550, −6.7079), initially assuming a single power-law (PL) spectral shape, $(dN/dE) = N_0(E/E_0)^{-\Gamma}$, and fitted the normalization ($N_0$) and photon index ($\Gamma$), while keeping the scale energy fixed to $E_0$=1 GeV. We performed a binned maximum likelihood analysis on a 35.3 × 35.3 square region, binning data into 0.1 × 0.1 pixels. For our starting model, the spectral parameters of sources within 15° of the ROI center were allowed to vary if they were found in the 4FGL-DR3 analysis with test statistic (Mattox et al. 1996), $TS \ge 100$, while the parameters of all other sources were held fixed. We left the normalization and index of the Galactic diffuse spectrum free in the fit, and the isotropic component also had a free normalization. To refine our starting model for the ROI before analyzing the outburst, we fit a 1 yr data set spanning 2020 July 19 to 2021 July 19, which ends 20 days before $t_0$. The fit was done over the energy range 0.05–300 GeV, including the effect of energy dispersion in our analysis (except for the isotropic component), and resulted in no significant detection at the position of the nova with a 95% confidence flux upper limit, $< 5.3 \times 10^{-8} \text{ ph cm}^{-2} \text{ s}^{-1}$ ($< 1.5 \times 10^{-8} \text{ ph cm}^{-2} \text{ s}^{-1}$ from 0.1–300 GeV). Using the results of this fit, we created a new model of the region for the following 10 day spectral analyses of the nova outburst. In this updated model, we fixed the spectral parameters of all the point sources in the ROI that had been free but had $TS < 100$ in the 1 yr analysis. For point sources with free parameters found to have $TS \ge 100$, we left the normalization parameters free but fixed all other spectral parameters. We also fixed the spectral parameters of all extended sources and the Galactic diffuse emission while leaving the normalization of the isotropic component free to vary.

Analyses of previous γ-ray novae found significant curvature in the spectrum (e.g., Ackermann et al. 2014) with the best-fit spectra using an exponentially cutoff power-law (ECPL) shape, $(dN/dE) = N_0(E/E_0)^{-\Gamma} \exp(-aE)$, where $a$ is the fitted exponential factor.<fnref n="23" /> To test for curvature in the γ-ray spectrum of RS Oph, we analyzed the time period from 2021 August 8.5 to 18.5, corresponding to the main activity when there were consecutive power-law fluxes, $F_γ$ ($>0.1$ GeV) $\ge 1 \times 10^{-6} \text{ ph cm}^{-2} \text{ s}^{-1}$ based on preliminary results (Cheung et al. 2021b).<fnref n="24" /> We performed a binned maximum likelihood analysis on this 10 day period over the 0.05–300 GeV energy

<footnote n="22">²² Both files are available for download at https://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html.</footnote>

<footnote n="23">²³ This functional form of an ECPL is more stable than other available models in https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/source_models.html, and a cutoff energy can be derived as $E_{cut} = a^{-1}$.</footnote>

<footnote n="24">²⁴ Following the initial LAT detection of RS Oph on August 8.759.00 (§ 1), we continued to monitor it through an automated daily analysis pipeline started by the LAT team in March 2021 aimed at detecting the anticipated nova outburst independent of other instruments. The main differences between this preliminary analysis and the final results presented are the background source model (4FGL), the time range of the data set used to fit the background, and a single PL spectral model was adopted for RS Oph throughout.</footnote>

<footnote n="19">¹⁹ See the talk given by Tatischeff & Hernanz (2008) at https://www.astro.keele.ac.uk/rsoph/pdfs/tatischeff.pdf.</footnote>

<footnote n="20">²⁰ The LAT automated science processing (ASP; Atwood et al. 2009, see Section 2.6.3 therein) pipeline for these data completed on 2021 August 9, 02:02 UT and was publicly reported on August 9, 05:05 UT.</footnote>

<footnote n="21">²¹ See https://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Pass8_usage.html.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Adaptively smoothed LAT count maps of 0.2–5 GeV energy photons centered on the optical position of RS Oph (marked with a 1° radius circle). For visualization purposes, the time intervals in the panels were selected to have similar exposure before (a: July 29.0 to August 6.5) and after (b: August 9.0 to 11.8) its 2021 outburst. Note the structured Galactic diffuse emission in the bottom portion of the images.</figcaption></figure>

range, similar to what was done for the 1 yr period prior to the outburst and starting from the updated model resulting from that fit. The nova was modeled using both an ECPL and a single PL, deriving a TS curve = -2($\mathcal{L}_{ECPL} - \mathcal{L}_{PL}$) as a measure of significance of curvature in the spectrum, where $\mathcal{L}$ are the respective best-fit–log likelihood values of the two fits. The results show ~8σ (TS$_{curve}$ = 66) evidence for curvature in the spectrum with best-fit ECPL parameters (TS = 2856), $\Gamma$ = 1.66 ± 0.05, and $E_{cut}$ = 6.0 ± 1.2 GeV, with 0.05–300 GeV and 0.1–300 GeV fluxes = (3.60 ± 0.20) × 10$^{-6}$ ph cm$^{-2}$ s$^{-1}$ and (2.16 ± 0.09) × 10$^{-6}$ ph cm$^{-2}$ s$^{-1}$, respectively. A significant detection in the lowest-energy bin from 50 to 100 MeV (TS=22, or 4.7σ) helped to constrain the curvature. The spectral curvature is likely due to the hadronic origin of the emission (see Section 4).

Though the correlated variability provides a firm identification of the LAT source with the nova, we used the *fermitools gtfindsrc* to localize the γ-ray emission during the first 10 day main activity interval. We selected 1–300 GeV photons because these events have the best per-photon resolution (Abdollahi et al. 2020) while also excluding the energy range where the Galactic diffuse emission is dominant. The resulting R.A., decl. = (267°.558, -6°.736) is offset by 0°.028 from the optical position, and within the 95% confidence-level containment radius of $r_{95}$ = 0°.029 (statistical only).

## 2.1. LAT Light Curves

We generated γ-ray flux light curves starting 20 days before the optical eruption epoch, $t_0$ = 2021 August 8.5 (Munari & Valisa 2021, see Section B.1) and ending at $t_0$ + 92 days in 6 hr, 1 day, and 4 day bins. The end date of day 92 (2021 November 8.5) probes ~ 2 × further than the last significant (≥ 3σ) LAT detection at day ~45 (see below) and coincides with the end of the Swift observing season (Section C). We produced light curves at >0.1 GeV energies to facilitate comparison to previous LAT-detected novae (e.g., Ackermann et al. 2014; Cheung et al. 2016); thus we restricted the events to 0.1 to 300 GeV photons within 15° of the ROI center for these analyses. In each time bin, we performed a binned maximum likelihood analysis on a 21°.2 × 21°.2 region, binned into pixels

0°.1 on a side. We adopted the best-fit background model from the analysis of the 10 day data set described in Section 2, while freezing the spectral normalizations of the field point sources with TS < 80 in that analysis, and removed sources more than 25° from the ROI center. For the nova, we used an ECPL model with a fixed $a$ parameter (corresponding to the best-fit $E_{cut}$ = 6 GeV from the 10 day interval) while fitting the normalization and the photon index. The resultant LAT 6 hr, 1 day, and 4 day bin >0.1 GeV flux light curves and fitted photon indices are presented in the Appendix (Section A.1).

For the first three days of the outburst, we performed further analyses on individual Fermi spacecraft orbits (~1.6 hrs). The orbit-binned analysis more accurately reflects the details of the LAT exposure profile at the RS Oph position at early times when the source was brightest (see Figure 2, top panel). Specifically, each of the analyzed orbital intervals consist of ~1 hr (~0.04 days) of exposure every ~3.2 hr (see Section A.2 for details).

Composite results from the different time-binned results are shown Figures 2 and 3. Throughout, we consider significant detections at TS ≥ 12 (≥ 3σ significance for two degrees of freedom (dof)), as indicated with black points in each LAT data panel. To help visualize fluxes at lower significance, particularly in the 4 day analysis, we use gray points with error bars to indicate intervals with TS = 6–12 (2–3σ). We report 95% confidence-level flux upper limits for time bins in which the nova was found with TS < 6 (<2σ), less than four predicted counts ($N_{pred}$), or had determined uncertainties greater than the fitted fluxes. The results derived from the light curves and spectral analysis of four defined γ-ray emission phases are given in Sections 3 and 4, respectively.

## 2.2. Shortest-timescale LAT Analysis

We searched for variability on the shortest timescales in the first 10 days of LAT data by applying the methods of Kerr (2019) to estimate the likelihood as a function of only the nova flux in 10 minute time bins (Figure 3, top panel). In this analysis, relative fluxes, $F_{γ,rel}$ (relative to the mean flux over the first 10 days set at 1.0) with error bars are shown for points with TS ≥ 4 corresponding to ≥ 2σ significances (significance
<figure><img src="image_5.png" /><figcaption>**Figure 2.** LAT γ-ray (top), AAVSO optical (middle), and Swift-XRT (bottom) light curves of RS Oph 2021 relative to $t_0$. The LAT light curve shows various binning as follows: orbit-timescale with typical exposure of 60 minutes up to day 3.0, then 6 hr up to day 6.0, 1 day up to day 12.0, and 4 day thereafter; data with error bars (black) are >3σ detections, while 2–3σ data and upper limits are shown in gray. The optical data are V-band (green) and visual observations (black) from the AAVSO; to reduce scatter, the visual points from day 1.3 onward are 5 point median values. Additional V-band measurements from days 0.233–0.381 and 1.227–1.375 are from observations described in Section B.1. The Swift-XRT data are in the 2–10 keV band.</figcaption></figure>

~√TS for 1 dof). We subsequently grouped these intervals into longer partitions with the Bayesian blocks (BB) algorithms (Scargle et al. 2013). Additionally, we estimated the likelihood for each spacecraft orbit bin (not shown) and confirmed the results of the *gtlike* analysis reported above.

The 10 minute analysis indicates a potential additional ~20 minute duration feature at day 2.21 with $F_{\gamma,rel} = 3.45 \pm 0.54$. This is >2×brighter than the measurements on the preceding day 1.94–2.20 and subsequent day 2.22–3.80 bins with respective $F_{\gamma,rel} = 1.35 \pm 0.21$ and $1.50 \pm 0.10$. This feature also appears in the orbit-binned analysis at day 2.207 (60 minute bin), albeit at a slightly less pronounced level, with a flux, $F_{\gamma}(>0.1 \text{ GeV}) = (5.8 \pm 1.2) \times 10^{-6} \text{ ph cm}^{-2} \text{ s}^{-1}$, while the adjacent orbit bin fluxes are ≈(3.4 ± 1.0) × 10⁻⁶ ph cm⁻² s⁻¹, amounting to ~1.5σ differences. Assessing the significance of features obtained with the BB algorithm is challenging. There are 431 individual 10 minute intervals in the analyzed range, and we adopted an exponential prior on the number of change points $\pi(N_p)^{-\gamma}$ with $\gamma = 7$, yielding a “false-positive rate” for a spurious change point of 0.4. However, if we restrict attention to the first three days when variations can be more easily detected when the source was brightest, the false-positive rate drops to ~0.1. Thus we estimate the total significance of this ~20 minute duration feature at day 2.21 (involving two additional change points) to be about 2σ.

We additionally split select orbits into two bins and performed the *gtlike* analysis (see Section A.2). In only one orbit bin was there observed flux variability between the

two 30 minute intervals, indicating an observed peak flux, $F_{\gamma}(>0.1 \text{ GeV}) = (8.1 \pm 1.8) \times 10^{-6} \text{ ph cm}^{-2} \text{ s}^{-1}$ at day 1.40, halving to (4.0 ± 1.2) × 10⁻⁶ ph cm⁻² s⁻¹ (~2σ difference). This feature was not observed in the BB analysis above.

## 2.3. Highest-energy LAT Photons

We searched for the highest-energy LAT events from $t_0$–20 to +92 days by selecting >5 GeV photons within 0.5°5 from the optical position of RS Oph. The resultant list of photons is given in the Supplement (Section A.1). The times and energies of the photons detected during the first 10 day main activity interval are shown in Figure 3.

There are eight >10 GeV LAT photons detected with *gtsrcprob* weight value (probability of association calculated using the best-fit ECPL model for the 10 day data set described in Section 2) of > 0.90 observed within 10 days of $t_0$. Curiously, the first of these events (10.5 GeV) at day 0.222 is in the orbit bin prior to the first orbit detection on day 0.352 (± 0.023). The highest energy photon is 23 GeV on day 2.861, just after the flux peak.

At later times, after the main 10 day activity interval (see the Figure in Section A.1), the highest-energy photons detected with *gtsrcprob* weight values > 0.90 were on days 13.8 and 19.2 (E = 11 and 19 GeV, respectively). Outside the 45 day LAT emission duration (Section 2.1), a single 12.8 GeV photon with smaller *gtsrcprob* weight value = 0.85 was found at day 55.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Selected data for the first 10 days of activity for RS Oph 2021. From top: LAT >0.1 GeV light curve in 10 minute bins and units of relative flux (to the mean flux over the first 10 days, set to 1.0) with upper limits indicated with gray arrows, and Bayesian Block partitions indicated in red, LAT >0.1 GeV composite light curve (see Figure 2, top), LAT photon indices, Γ, energies of the detected *E* > 5 GeV photons, optical light curve (see Figure 2, middle; except all individual visual observations are shown here), and X-ray (2–10 keV) light curve (see Figure 2, bottom).</figcaption></figure>

**3. Fermi LAT and Multiwavelength Results**

The LAT light curve for RS Oph 2021 shown in Figure 2 (top panel) is a composite of the orbit-binned, 6 hr, 1 day, and 4 day analysis relative to $t_0$, choosing increasing size time bins at later times. The LAT 6 hr and 1 day analyses show significant consecutive detections up to day 10. The longer, 4 day integration increased the sensitivity to lower fluxes and helped to define the light curve at later times, when the shorter-timescale analysis resulted in more intermittent and lower-significance detections.

The main LAT-based results determined from the LAT >0.1 GeV light curve are as follows: (a) the observed γ-ray onset is constrained to the orbit centered on day 0.35, with a flux, $F_{\gamma} = (2.1 \pm 0.7) \times 10^{-6}$ ph cm$^{-2}$ s$^{-1}$; (b) the source rises to a peak $F_{\gamma} = (5.7 \pm 1.2) \times 10^{-6}$ ph cm$^{-2}$ s$^{-1}$ at day 0.745; (c) the peak is flat, with an approximately constant (slope of $-0.18 \pm 0.17$), and average $F_{\gamma} \simeq 5 \times 10^{-6}$ ph cm$^{-2}$ s$^{-1}$ through to day 2.5; (d) the flux declined by a factor of 2 from the day 0.745 flux value at approximately day 3; (e) considering the last significant 6 hr detection on day 45.125 and the first orbital-bin detection at day 0.352, the total γ-ray duration was approximately 45 days; and (f) taking the power-law slope of $-1.53 \pm 0.11$ fit from day 2.5 to 45, and the average peak flux from day 0.745 to 2.5, the source declined by a factor of 10 on day 10.

Results from additional LAT analysis performed—the BB analysis, the spectral photon indices, and the highest-energy (*E* > 5 GeV) photons—are also presented in Figure 3 for the main 10 day activity phase. The above results (a) and (b) are consistent with the BB analysis with the first BB detection at day 0.350, the peak from 0.746 to 1.940, and the subsequent flux decline. The fitted LAT spectra have a wide range of spectral slopes,

$\Gamma = 1.5$ to 2.1, with typical errors of 0.2. In the rising portion of the light curve, the γ-ray spectral slopes were $\Gamma \sim 2.0-2.3$, while the spectra hardened ($\Gamma \sim 1.7$) during the peak. The hardest spectrum ($\Gamma = 1.4 \pm 0.2$) was observed on day 2.87 when the highest-energy LAT photon with $E = 23.3$ GeV was detected (along with an 8.6 GeV photon only 42.4 s later).

To compare to the LAT data, optical and X-ray (2–10 keV) light curves are shown in Figure 2 for the entire ~3 month interval studied while Figure 3 details the main activity interval during the first 10 days. Overall, the γ-ray and optical light curves are similar (see also Figure 3), both peaking early (around day 1; see below), while the X-rays peak later (day ~6).

The optical data (Section B.1) were mainly taken from the American Association of Variable Star Observers (AAVSO) database, with additional serendipitous measurements on the first two nights (August 8, 9) obtained by a Global Meteor Network (Vida et al. 2021) camera IL0003 (by A.L. and A.B.) newly presented here. During the rise in γ-ray flux, the optical increased by four magnitudes in the V band from observations spanning about days 0.07 to 0.80. The power-law slope fitted to the earliest γ-ray data up to day 0.80 ($1.58 \pm 0.69$) is consistent with the optical V-band one ($1.28 \pm 0.01$). The γ-ray onset observed on 2021 August 8.852 was delayed by ~0.35 day after the optical eruption ($t_0$), but earlier than the visual discovery epochs by about 0.08 day. The optically brightest emission was observed from day 0.9 to 1.3 (visual = 4.5–4.6 mag), within the wider time span of the observed γ-ray peak fluxes from day 0.75 to 1.67 (see Table 3). Fitting a broken power law to the visual measurements from days 0.37 to 3.0, the best-fit peak was at day 1.09 ($\pm 0.12$), consistent with the day 1.08 ($\pm 0.05$) estimated by Munari & Valisa (2021). Although the γ-ray peak fluxes appear constant over a ~1 day interval, assuming a broken
<figure><img src="image_7.png" /><figcaption>**Figure 4.** (Left): Gamma-ray luminosity as a function of the observed optical luminosity (see Section B.2) for RS Oph 2021 extracted in 4 day time bins. The black solid line shows the best-fit ratio with its uncertainty (dashed lines) obtained by excluding the largest luminosity values (see text). (Right): Ratio of the γ-ray luminosity to the observed optical luminosity extracted in 4 day time bins. The black solid line shows the best-fit ratio, $(2.81 \pm 0.15) \times 10^{-3}$, with its uncertainty (dashed lines), obtained without the ratio in the bin at day 0–4. The insets in each panel show the results from reanalysis of the highest-luminosity, earliest 4 day bin in smaller time intervals (blue data points; see text).</figcaption></figure>

power-law parameterization of the composite LAT light curve (Figure 2) gives a best-fit peak at day 1.64 (± 0.11) that is formally delayed with respect to the optical peak. The slopes of the γ-ray and optical declines, -1.35 ± 0.07 and -1.395 ± 0.006, respectively, are consistent given the uncertainties.

We examined the best-fit γ-ray luminosity in 4 day bins obtained from the LAT data analysis<fnref n="25" /> as a function of the observed optical luminosity estimated in 4 day bins using the AAVSO data (see Section B.2). Both luminosities are proportional to each other, except for the largest luminosity value (Figure 4, left). The ratios of the γ-ray to optical luminosities in RS Oph (Figure 4, right) are similar to those derived for classical novae (Metzger et al. 2015; Li et al. 2017; Aydi et al. 2020). Excluding the first 4 day bin, we calculated a luminosity ratio $L_γ(E > 50 \text{ MeV})/L_{\text{opt}} = (2.81 \pm 0.15) \times 10^{-3}$, while the ratio in the first 4 day bin is lower by a factor of ~1.5. Reanalysis of the first 4 day bin in smaller time intervals (day 0–1, 1–2, and 2–4; blue data points in panel insets in Figure 4) found a ratio lower by a factor of ~5 in the first 1 day bin (see Section 4 for a discussion), while the other ratio values are compatible with the constant value.

The X-ray observations obtained with the X-ray Telescope (XRT; Burrows et al. 2005) on the Neil Gehrels Swift Observatory (Gehrels et al. 2004) are summarized in Section C, and their detailed analysis is presented elsewhere (Page et al. 2022). The X-ray light curve for the 2021 outburst is overall similar to that seen following the 2006 eruption (see Bode et al. 2006; Osborne et al. 2011), except that the 2006 data count rate maximum was about double that measured in 2021 (Page et al. 2022). Here, the XRT 2–10 keV light curves are used to compare to the LAT γ-ray and optical light curves. The early X-ray emission at >2 keV energies is dominated by

Bremsstrahlung from shocked gas in the nova ejecta (e.g., Sokoloski et al. 2006), and the results of the X-ray temperature spectral fits of the 2021 XRT data presented by Page et al. (2022) are used to constrain the temporal evolution of the ejecta velocity in our modeling (see Section 4 and Appendix C). The main feature of the 2021 XRT 2–10 keV light curve is that the peak at day 6.4 ± 0.1 is consistent with the approximate time of the break in the X-ray derived velocity profile (Section 4). Thereafter, the 2–10 keV light curve decline can be parameterized as a broken power law with a slope = -1.24 ± 0.03 to day 19.1 (± 0.2) and a steeper decline (slope = -3.11 ± 0.03) to the last data at day 87.6.

## 4. Pion-decay Gamma-ray Emission

In the context of the hadronic model, we define four emission phases for spectral study with the >50 MeV LAT data (Figure 5)—the rise from day 0 to 1.0, a peak from day 1.0 to 2.75, decline-a from day 3.0 to 9.0 (2–10×smaller than the peak), and decline-b from day 9.0 to 46.0. The aim is to compare the evolution of the flux and spectral properties of the γ-ray emission in defined observation periods that are long enough to have sufficient statistics to determine accurate values of the spectral parameters. The Rise phase in particular was extended to include the initial portion of the peak to increase statistics in that bin (see, Figure 6, right).

The γ-ray spectrum of the hadronic model is calculated following the method described in Kamae et al. (2006), assuming that the energy distribution of the high-energy protons is a power law in proton momentum multiplied by an exponential cutoff (see the supplementary material of Ackermann et al. 2014). The parameters fit to the LAT data are the normalization, the slope ($s_p$) and the cutoff energy ($E_{\text{cp}}$) of the high-energy proton spectrum. The fits for the four intervals were performed by maximizing the likelihood and the best-fit normalization is used to calculate the γ-ray spectral energy

<footnote n="25">²⁵ The γ-ray luminosities were calculated by fitting the fluxes and the proton spectrum slopes (without energy cutoff) in the $\pi^0$ model to the Fermi LAT data (see Section 4).</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 5.** LAT γ-ray spectra and best-fit hadronic models for the rise (top left), peak (top right), decline-a (bottom left), and decline-b (bottom right) phases for RS Oph 2021—see Table 1. Vertical lines indicate $1\sigma$ uncertainties when $TS \ge 4.0$ and arrows indicate $2\sigma$ upper limits when $TS < 4.0$.</figcaption></figure>

(The data used to create this figure are available.)

distributions (SEDs; Figure 5). The corresponding best-fit parameters obtained for each observation period are presented in Table 1 along with the difference between the TS values resulting from the best-fit hadronic and ECPL models ($\Delta$TS). According to the $\Delta$TS values, the hadronic model provides a fit as good as the one obtained with the ECPL model for most of the considered observation periods except for the decline-a period where the former is preferred over the latter. Between these four observation periods, the best-fit ECPL spectral parameters do not change significantly within the uncertainties which suggests that there is no significant variation in the spectral shape. Comparing the individual values to the averaged best-fit ones from the first 10 days after the outburst

($\Gamma = 1.66 \pm 0.05$, $E_{\text{cut}} = 6.0 \pm 1.2$ GeV; Section 2), suggests a modest spectral change in the fit for the subset of data for the peak period ($E_{\text{cut}} = 3.16_{-0.07}^{+1.10}$ GeV). The slopes of the proton spectra are compatible with a constant value of $\sim 2.4$ and there is no significant energy cutoff (i.e., lower limits are derived).

The hadronic model can explain the measured LAT spectra of RS Oph 2021 (see Figure 5) and has also been proposed to explain its very-high-energy (VHE; $>0.1$ TeV) γ-ray emission (H.E.S.S. Collaboration et al. 2022; MAGIC Collaboration et al. 2022). Therefore, we constructed a simplified light-curve model assuming that the γ-ray emission results from $\pi^0$-decay produced by high-energy protons interacting with the material of the ejecta as proposed by Tatischeff & Hernanz (2007) to
<figure><img src="image_9.png" /><figcaption>**Figure 6.** (Left): Maximum proton energy with time for the hadronic model for RS Oph 2021. (Right): Comparison of the hadronic model with the observed γ-ray flux light curve presented in Figure 2. The dashed line shows the date at which the ejecta reach the radius of 1.5×the binary separation (see text).</figcaption></figure>

model the RS Oph 2006 outburst; see also Hernanz & Tatischeff (2012). The protons are accelerated via the Fermi process in the shock between the nova ejecta as it propagates through the RG wind. We estimated the speed of the ejecta with time using an analytical model fitted to the velocities derived from the X-ray temperatures (see Bode et al. 2006) measured with the Swift-XRT data by Page et al. (2022)—see details in Appendix C. This velocity model is representative of the average variation in the shock velocity with time and is used to compute the time evolution of the nova shell radius. The model can be roughly described by a constant value of 2470 km s⁻¹ for *t* < 6 days and proportional to *t*⁻⁰.⁴³ at *t* > 6 day (which is close to the *t*⁻⁰.⁵ variation used by Tatischeff & Hernanz 2007). The density profile of the RG wind is taken from Tatischeff & Hernanz (2007). For simplicity, the RG wind density is assumed to be uniform within a radius of <1.5×the binary separation of 1.5. au (Fekel et al. 2000) because of the complexity of the matter distribution inside and around the binary system.

We calculated the maximum energy of accelerated protons as a function of time (Figure 6, left) by integrating the sum of energy loss and gain rates. The energy gain rate is computed as in Tatischeff & Hernanz (2007), with a compression ratio of the shock of 4 (strong adiabatic shock) and a magnetic field estimated assuming equipartition in the compressed gas (from the RG wind) with a wind temperature of 10⁴ K. The energy loss rate takes into account Coulomb collisions and inelastic p–p collisions in the compressed gas. They were estimated from the methods described in Mannheim & Schlickeiser (1994) and Dermer & Powale (2013), respectively. In our model, the maximum proton energy as a function of time is similar to that in (Tatischeff & Hernanz 2007; Figure 2, therein) except the maximum energy is ~7.4 TeV instead of ~3 TeV because in our case the acceleration started before day 1, while it started at day 1 in Tatischeff & Hernanz (2007). The maximum proton energies estimated are >1 TeV starting at 0.3 days after the outburst, in agreement with the lower limits on the cutoff energy of the proton spectrum derived from the LAT data (see Table 1). In our model, the maximum proton energy changes from ~2 to ~5 TeV between days 1.3 and 5.4, the

times spanned by the early H.E.S.S. VHE observation periods of RS Oph 2021. This could explain the hardening observed in the VHE spectra and the increasing high-energy photon detections from ~0.3 to 1 TeV between these two epochs (H. E.S.S. Collaboration et al. 2022).

To calculate γ-ray light curves for RS Oph 2021 in our model, we assumed that an injection fraction $f_{\text{inj}}$ of the RG wind protons crossed by the shock is accelerated toward the expanding ejecta with a Fermi-type spectrum with a slope of –2.1, up to the maximum energy. This spectrum is used as a source function of a diffusion equation that estimates the time evolution of the proton spectrum, taking into account losses via inelastic collisions in the expanding ejecta. We neglect escape losses, which would have the effect of softening the proton spectrum as the highest-energy protons would be the first to escape the shock and would also require a larger injection fraction to fit the measured γ-ray light curve. The resulting spectra are then used to compute the rate of $\pi^0$ production (with the energy-dependent cross section of Dermer 1986) produced in collisions with the hadrons in the expanding ejecta, which is assumed to have a uniform density. The γ-ray fluxes were computed for a distance of 1.6 kpc and an ejecta mass of ~1.1 × 10⁻⁶ M☉ (see Bode et al. 2006). The resulting γ-ray model fluxes compare well with the observed 2021 LAT light-curve data (Figure 6, right), and is similar to that obtained by Hernanz & Tatischeff (2012) for the 2006 outburst. Our model was obtained for an injection fraction $f_{\text{inj}} \sim 4 \times 10^{-6}$, lower than the value in Tatischeff & Hernanz (2007). Considering the ejecta mass and velocity model adopted in our calculations, we estimated the attenuation of γ rays due to Compton scattering and pair creation from interactions with nuclei in the expanding ejecta during the earliest stages when the gas densities are highest<fnref n="26" /> (Martin et al. 2018, Section 3.2 therein). According to our model, the attenuation amounts to a factor of ~2 smaller observed flux at ~400 MeV (the peak of the SED) at day 0.2 and is negligible after day 0.6, so does not fully explain the factor of ~5 smaller γ-ray to optical luminosity ratio observed

<footnote n="26">²⁶ The Compton scattering and pair production processes dominate for ~100 MeV and ~100 MeV photons, respectively.</footnote>
<table><caption>**Table 1**
Fermi LAT γ-Ray Spectral and Hadronic Model Fit Results for RS Oph 2021</caption><thead><tr><th>Phase<br>Days</th><th>Rise<br>0.0–1.0</th><th>Peak<br>1.0–2.75</th><th>Decline-a<br>3.0–9.0</th><th>Decline-b<br>9.0–46.0</th></tr></thead><tbody><tr><td colspan="5" style="text-align: center; font-weight: bold;">Exponentially Cutoff Power Law</td></tr><tr><td>Photon Index, Γ</td><td>1.89 ± 0.15</td><td>1.50<sup>+0.09</sup><sub>-0.14</sub></td><td>1.89 ± 0.05</td><td>1.48 ± 0.34</td></tr><tr><td>E<sub>cut</sub> (GeV)</td><td>11.4 ± 12.1</td><td>3.16<sup>+1.10</sup><sub>-0.07</sub></td><td>78.6 ± 62.8</td><td>3.39 ± 2.69</td></tr><tr><td>Flux(0.05–300 GeV)</td><td>3.66 ± 0.64</td><td>7.36<sup>+0.37</sup><sub>-0.29</sub></td><td>2.45 ± 0.20</td><td>0.17 ± 0.06</td></tr><tr><td>TS</td><td>212.2</td><td>1585.1</td><td>1012.3</td><td>135.2</td></tr><tr><td colspan="5" style="text-align: center; font-weight: bold;">Hadronic Model</td></tr><tr><td>Slope, s<sub>p</sub></td><td>2.70<sup>+0.16</sup><sub>-0.27</sub></td><td>2.30<sup>+0.23</sup><sub>-0.37</sub></td><td>2.35<sup>+0.19</sup><sub>-0.23</sub></td><td>2.15<sup>+0.32</sup><sub>-1.14</sub></td></tr><tr><td>E<sub>cp</sub> (GeV)</td><td>>46</td><td>>28</td><td>>63</td><td>>8</td></tr><tr><td>Flux(0.05–300 GeV)</td><td>2.56<sup>+0.18</sup><sub>-0.23</sub></td><td>6.02<sup>+0.26</sup><sub>-0.33</sub></td><td>1.60<sup>+0.09</sup><sub>-0.08</sub></td><td>0.15<sup>+0.02</sup><sub>-0.01</sub></td></tr><tr><td>ΔTS</td><td>-2.3</td><td>8.9</td><td>38.0</td><td>-3.5</td></tr></tbody></table>

Note. The observation phases are defined in Section 4. For the ECPL model fit results, Γ and E<sub>cut</sub> are the best-fit photon index and the cutoff energy of the photon spectrum, respectively. For the hadronic model, the best-fit slope, s<sub>p</sub> and cutoff energy, E<sub>cp</sub> (2σ lower limits) of the high-energy proton spectrum are given. The fluxes at >50 MeV energies are in units of × 10<sup>-6</sup> ph cm<sup>-2</sup> s<sup>-1</sup>. Values are given with their 1σ statistical uncertainties. ΔTS is the difference between the TS values of the hadronic model and the ECPL model.

in the first day (see Figure 4). The smaller ratio could be more fully accounted for by the lower high-energy proton fluxes at the earliest acceleration phase in the shock development (see Figure 6, left) resulting in a lower production of γ rays.

## 5. Discussion and Summary

The anticipated γ-ray detection of the outburst of the recurrent nova RS Ophiuchi with Fermi LAT was realized in 2021 (Hernanz & Tatischeff 2012). With a peak flux, F<sub>γ</sub> (>0.1 GeV)≈6 × 10<sup>-6</sup> ph cm<sup>-2</sup> s<sup>-1</sup>, the RS Oph 2021 outburst is the brightest nova detected thus far in γ rays by the LAT. The previous brightest LAT-detected nova was the classical nova V906 Car 2018 with a peak F<sub>γ</sub> (>0.1 GeV)≈4 × 10<sup>-6</sup> ph cm<sup>-2</sup> s<sup>-1</sup> (Jean et al. 2018; Aydi et al. 2020), although its early γ rays were missed due to an anomaly with a Fermi solar panel assembly motor (Abdollahi et al. 2022). The detection of >5 GeV and up to ~20–23 GeV energy photons (Section 2.3) help constrain the high-energy portion of the SEDs in our modeling of select intervals. RS Oph’s brightness in the LAT band extended to higher energies, resulting in the first VHE detection of a nova (H.E.S.S. Collaboration et al. 2022; MAGIC Collaboration et al. 2022). The VHE observations constrain the maximum energy of the accelerated protons in this nova to be ~10 TeV (H.E.S.S. Collaboration et al. 2022), in agreement with the value derived in our semi-analytical model (see Section 4).

The LAT-observed >0.1 GeV γ rays peaked at ~1 day after the optical outburst and could reflect the time needed for the shock to accelerate enough protons to significant energy while it propagates through the inner part of the binary system. However, we cannot exclude that this timescale is also due to a change in the γ-ray attenuation produced by the material in which the shock propagates during the first day after optical eruption. The LAT light curve showed evidence, albeit at low significances, of factors of 2 fluctuations in the flux on timescales of 30–200 minutes during the first few days. Such variability may be expected because of density variations of the swept-up material by the shock in the inner part of the binary system and/or changes in the outflow (see O’Brien et al. 2006; Rupen et al. 2008; Sokoloski et al. 2008; Walder et al. 2008;

Orlando et al. 2009). The decline in the γ-ray emission starting at days 2–3 could correspond to the timescale at which the shock enters a region where the density of the RG wind decreases with the shock radius.

The γ-ray detection of RS Oph can be compared to previous novae in symbiotic binaries detected by the LAT, particularly the first LAT-detected nova V407 Cyg 2010 (Abdo et al. 2010). The RS Oph LAT detection was expected based on external shocks from the nova ejecta evolving in the dense RG wind (Tatischeff & Hernanz 2007; Hernanz & Tatischeff 2012), while the V407 Cyg detection was unanticipated because it was a less-studied symbiotic binary. The 6 hr LAT peak flux, F<sub>γ</sub> (>0.1 GeV)≈6 × 10<sup>-6</sup> ph cm<sup>-2</sup> s<sup>-1</sup> in RS Oph was ~4× brighter than in V407 Cyg, consistent with the factor of ~2 distance scaling if a 2.7 kpc distance is adopted for the latter. V407 Cyg was not detected in the VHE band despite past observations performed by Cherenkov telescopes (see Aliu et al. 2012; Ahnen et al. 2015; López-Coto et al. 2015). Other symbiotic recurrent novae with optical outbursts observed with the Fermi LAT were V745 Sco 2014 (consisting of two low-significance 1 day detections coincident with its optical peak; Cheung et al. 2014; Franckowiak et al. 2018) and the >5σ LAT detection of V3890 Sgr 2019 (Buson et al. 2019). Interestingly, the nova V1535 Sco 2015 was a low-significance (2–3σ) LAT detection as well (Franckowiak et al. 2018) and is proposed to be a symbiotic system (Srivastava et al. 2015; Linford et al. 2017). The latter three systems have distances at ~6 kpc and greater, so are broadly consistent with distance-scaled fluxes.

The majority of LAT detections are of classical novae, involving less-dense circumbinary environments from their main-sequence companions, at a rate of ~1 yr<sup>-1</sup> starting in 2012 (Ackermann et al. 2014; Cheung et al. 2016); see Chomik et al. (2021) for a review. These observations demonstrated the importance of internal shocks inside the nova ejecta in producing γ rays (e.g., Martin et al. 2018; Vurm & Metzger 2018). The fact that symbiotic recurrent novae ejecta must interact with the RG wind makes the case that both internal and external shocks are important in novae. In RS Oph 2021, it appears that the hadronic processes from the external shock discussed in Section 4 are the dominant
component (see H.E.S.S. Collaboration et al. 2022; MAGIC Collaboration et al. 2022).

With all the known symbiotic recurrent novae detected in outburst during the LAT era, the only known system remaining is T CrB, at a distance of only 0.8 kpc. Its next outburst could be as soon as in 2023–2026 (Schaefer 2010, 2019; Luna et al. 2020), and if the distance-scaled fluxes hold, it should be remarkably bright with fluxes, $F_γ$ (>0.1 GeV) $\gtrsim 1 \times 10^{-5}$ ph cm$^{-2}$ s$^{-1}$; thus it will be studied in remarkable detail with the LAT. It is also conceivable that the LAT could detect the next RS Oph outburst because its recurrence interval has been as short as 9 yr (thus 2030 or later). In these cases (see Hernanz & José 2008, for a discussion of RS Oph), the radioactive decay emission at MeV energies could also be observed with the Compton Spectrometer and Imager (Tomsick et al. 2022), which will operate as an all-sky survey similar to Fermi, and is expected to launch in 2026. This could build the most-complete MeV–GeV picture of the different $\gamma$-ray components of a nova evolution.

The Fermi LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat à l'Energie Atomique and the Centre National de la Recherche Scientifique / Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council and the Swedish National Space Board in Sweden.

Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d'Études Spatiales in France. This work performed in part under DOE Contract DE-AC02-76SF00515.

Research at the Naval Research Laboratory is supported by NASA DPR S-15633-Y. K.L.P. and A.P.B. acknowledge funding from the UK Space Agency.

This work made use of data supplied by the UK Swift Science Data Centre at the University of Leicester. We acknowledge with thanks the variable star observations from the AAVSO International Database contributed by observers worldwide and used in this research.

We thank the anonymous referee and A. Azzollini for their comments on the manuscript.

Facilities: Fermi, Swift, AAVSO.

# Appendix A
## Additional Details of the LAT Analysis

### A.1. Light-curve Data and Highest-energy LAT Photon List

The full results of the LAT >0.1 GeV light-curve spectral analysis in 6 hr, 1 day, and 4 day bins described in Section 2.1 are presented in Figure 7.

The list of LAT photons with $E > 5$ GeV found in the analysis described in Section 2.3 is given in Table 2 and plotted in Figure 7 (bottom panel).

### A.2. LAT Orbit-binned and Split-orbit Spectral Analysis

Fermi achieves its all-sky exposure profile by alternating between north and south scans of the sky during each spacecraft orbit (~1.6 hr). The bulk of the exposure at the RS Oph position comes from the alternating southern sky scans. For the LAT orbit-binned analysis, we define individual orbits by calculating the zenith angle of the center of our ROI, in 30 s steps, and looking for points where the angle reaches a maximum and then begins decreasing. We refer to these maxima as “orbit midnights” and use them to define orbit bin endpoints. We also calculate the angle ($\theta$) between the center of our ROI and the LAT z-axis and exclude from our final analysis those orbits where $\theta$ is always greater than 60° and those orbits with fewer than 10 events. In each orbit, we used the same radius ROI and energy selections described previously in Section 2.1. When calculating the exposure, we considered the instrument azimuth angle, using phibins=5 in the fermi-tool gtltcube and then performed an unbinned likelihood analysis, with only the isotropic diffuse emission component and the normalization and $\Gamma$ value of the nova free to vary. This approach resulted in 22 orbital bins (Table 3), over the time period from $t_0$ to $t_0 + 3$ days, with average exposure lengths of 57.5 minutes (0.04 days) and average offsets of the bin centers of 190.4 minutes (0.13 days). The exposure lengths were estimated from an exposure light curve for a 5° radius selection around our ROI center, in 5 minute bins, made using the fermitools gtbin and gtexposure. All detections were $TS > 25$ (>4$\sigma$ for 2 dof) with 95% confidence upper limits derived for the first two bins (days 0.088 and 0.221).

We additionally conducted a “split orbit” analysis with gtlike, selecting nine orbit bins of particular interest. The orbits include seven high-flux bins spanning the observed peak (days 0.745 and days 1.012–1.668), and two later high-flux bins (days 2.207 and 2.594). The exposure profile over an orbit is not generally uniform, so we divide it by examining the exposure accumulation and determining the time that yields equal exposure. The results are presented in Table 4. The most significant variability was a factor of 2 drop at day 1.41 as described in Section 2.2. For a single orbit (days 0.745), the photon index $\Gamma$ appears to increase by $1.0 \pm 0.5$ (~2$\sigma$).
<figure><img src="image_10.png" /><figcaption>**Figure 7.** From top to bottom: LAT >0.1 GeV 6 hr light curve and photon index, 1 day light curve and photon index, 4 day light curve and photon index, and energies of detected >5 GeV photons. Fluxes with error bars (black) are shown for TS ⩾ 12 (3σ) points, while those with TS from 6.0 up to 12.0 (2–3σ) and upper limits are shown in gray. Photon indexes are shown only for TS ⩾ 12 points.</figcaption></figure>

<table><caption>**Table 2**
Highest-energy LAT Photons Selected at $E > 5$ GeV</caption><thead><tr><th>Day</th><th>MET (sec)</th><th><i>E</i> (GeV)</th><th>gtsrcprob</th></tr></thead><tbody><tr><td>0.222</td><td>650136003.9</td><td>10.54</td><td>0.957</td></tr><tr><td>0.887</td><td>650193475.7</td><td>8.90</td><td>0.998</td></tr><tr><td>1.146</td><td>650215820.9</td><td>5.07</td><td>0.999</td></tr><tr><td>1.151</td><td>650216236.4</td><td>5.19</td><td>0.999</td></tr><tr><td>1.155</td><td>650216600.8</td><td>6.01</td><td>0.999</td></tr><tr><td>1.411</td><td>650238679.9</td><td>12.86</td><td>0.998</td></tr><tr><td>1.546</td><td>650250356.7</td><td>12.12</td><td>0.999</td></tr><tr><td>1.681</td><td>650262005.6</td><td>8.61</td><td>0.999</td></tr><tr><td>2.597</td><td>650341199.0</td><td>6.55</td><td>0.968</td></tr><tr><td>2.740</td><td>650353560.9</td><td>5.65</td><td>1.000</td></tr><tr><td>2.861</td><td>650363956.4</td><td>23.32</td><td>0.999</td></tr><tr><td>2.861</td><td>650363998.8</td><td>8.57</td><td>0.997</td></tr><tr><td>2.863</td><td>650364171.6</td><td>5.18</td><td>0.969</td></tr><tr><td>3.260</td><td>650398494.0</td><td>9.53</td><td>0.999</td></tr><tr><td>3.391</td><td>650409765.4</td><td>14.12</td><td>0.999</td></tr><tr><td>3.518</td><td>650420734.2</td><td>6.84</td><td>0.999</td></tr><tr><td>3.645</td><td>650431699.5</td><td>7.48</td><td>0.996</td></tr><tr><td>3.859</td><td>650450260.5</td><td>18.80</td><td>0.971</td></tr><tr><td>4.049</td><td>650466613.6</td><td>6.69</td><td>0.992</td></tr></tbody></table>

<table><caption>**Table 2**
(Continued)</caption><thead><tr><th>Day</th><th>MET (sec)</th><th><i>E</i> (GeV)</th><th>gtsrcprob</th></tr></thead><tbody><tr><td>4.721</td><td>650524733.6</td><td>9.83</td><td>1.000</td></tr><tr><td>5.774</td><td>650615677.8</td><td>9.17</td><td>0.996</td></tr><tr><td>6.295</td><td>650660671.4</td><td>5.48</td><td>0.999</td></tr><tr><td>6.692</td><td>650694990.2</td><td>14.03</td><td>1.000</td></tr><tr><td>6.703</td><td>650695986.1</td><td>5.78</td><td>0.998</td></tr><tr><td>7.034</td><td>650724541.4</td><td>16.54</td><td>0.881</td></tr><tr><td>7.815</td><td>650791983.0</td><td>18.02</td><td>0.916</td></tr><tr><td>9.090</td><td>650902208.1</td><td>5.95</td><td>0.999</td></tr><tr><td>9.471</td><td>650935065.6</td><td>5.70</td><td>0.999</td></tr><tr><td>13.776</td><td>651307020.6</td><td>11.35</td><td>0.978</td></tr><tr><td>16.147</td><td>651511940.3</td><td>5.96</td><td>1.000</td></tr><tr><td>19.205</td><td>651776094.4</td><td>19.20</td><td>0.996</td></tr><tr><td>54.945</td><td>654864060.0</td><td>12.83</td><td>0.847</td></tr></tbody></table>

**Note.** Days are relative to $t_0$. The mission elapsed time (MET) for Fermi is defined as seconds since 2001.0 UTC.
<table><caption>**Table 3**
Orbit-binned LAT Analysis in the First Three Days</caption><thead><tr><th>Day</th><th>ΔT (day)</th><th>TS</th><th>N<sub>pred</sub></th><th>Flux (0.1–300 GeV) (× 10<sup>-6</sup> ph cm<sup>-2</sup> s<sup>-1</sup>)</th><th>Γ</th></tr></thead><tbody><tr><td>0.088</td><td>0.023</td><td>0</td><td>0.0</td><td>&lt;0.5</td><td>—</td></tr><tr><td>0.221</td><td>0.023</td><td>5.3</td><td>1.4</td><td>&lt;0.6</td><td>—</td></tr><tr><td>0.352</td><td>0.023</td><td>33</td><td>18.0</td><td>2.1 ± 0.8</td><td>2.1 ± 0.3</td></tr><tr><td>0.484</td><td>0.023</td><td>28</td><td>17.5</td><td>2.1 ± 0.8</td><td>2.1 ± 0.3</td></tr><tr><td>0.611</td><td>0.017</td><td>36</td><td>20.2</td><td>3.2 ± 1.1</td><td>2.0 ± 0.3</td></tr><tr><td>0.745</td><td>0.019</td><td>95</td><td>39.6</td><td>5.7 ± 1.2</td><td>2.3 ± 0.2</td></tr><tr><td>0.880</td><td>0.023</td><td>92</td><td>38.0</td><td>4.0 ± 1.0</td><td>1.7 ± 0.2</td></tr><tr><td>1.012</td><td>0.023</td><td>140</td><td>48.4</td><td>5.5 ± 1.1</td><td>1.9 ± 0.2</td></tr><tr><td>1.146</td><td>0.021</td><td>230</td><td>49.1</td><td>4.7 ± 0.9</td><td>1.4 ± 0.1</td></tr><tr><td>1.278</td><td>0.021</td><td>120</td><td>43.7</td><td>4.7 ± 1.0</td><td>1.7 ± 0.2</td></tr><tr><td>1.410</td><td>0.021</td><td>180</td><td>54.9</td><td>5.9 ± 1.0</td><td>1.6 ± 0.1</td></tr><tr><td>1.535</td><td>0.014</td><td>120</td><td>29.0</td><td>5.1 ± 1.3</td><td>1.7 ± 0.2</td></tr><tr><td>1.668</td><td>0.016</td><td>150</td><td>42.5</td><td>6.1 ± 1.2</td><td>1.7 ± 0.2</td></tr><tr><td>1.804</td><td>0.019</td><td>110</td><td>32.2</td><td>3.6 ± 0.9</td><td>1.5 ± 0.2</td></tr><tr><td>1.939</td><td>0.023</td><td>130</td><td>43.1</td><td>4.9 ± 1.0</td><td>1.8 ± 0.2</td></tr><tr><td>2.071</td><td>0.023</td><td>49</td><td>27.3</td><td>3.4 ± 1.0</td><td>2.1 ± 0.3</td></tr><tr><td>2.207</td><td>0.019</td><td>150</td><td>44.9</td><td>5.8 ± 1.2</td><td>1.9 ± 0.2</td></tr><tr><td>2.338</td><td>0.019</td><td>58</td><td>26.5</td><td>3.4 ± 0.9</td><td>1.8 ± 0.2</td></tr><tr><td>2.467</td><td>0.023</td><td>82</td><td>32.2</td><td>3.9 ± 0.9</td><td>1.7 ± 0.2</td></tr><tr><td>2.594</td><td>0.014</td><td>99</td><td>34.6</td><td>6.7 ± 1.5</td><td>2.0 ± 0.2</td></tr><tr><td>2.727</td><td>0.016</td><td>62</td><td>23.7</td><td>3.5 ± 1.0</td><td>1.9 ± 0.2</td></tr><tr><td>2.865</td><td>0.021</td><td>93</td><td>24.0</td><td>2.5 ± 0.8</td><td>1.4 ± 0.2</td></tr></tbody></table>

Note. Days are relative to $t_0$, and $\Delta T$ represents the half bin width.

<table><caption>**Table 4**
Split Orbit-binned LAT Analysis Results For Select Orbits</caption><thead><tr><th>Day</th><th>ΔT (day)</th><th>TS</th><th>Flux (0.1–300 GeV) (× 10<sup>-6</sup> ph cm<sup>-2</sup> s<sup>-1</sup>)</th><th>Γ</th></tr></thead><tbody><tr><td>0.7378</td><td>0.0121</td><td>57</td><td>4.5 ± 1.5</td><td>1.9 ± 0.3</td></tr><tr><td>0.7569</td><td>0.0070</td><td>42</td><td>6.4 ± 1.8</td><td>2.9 ± 0.5</td></tr><tr><td>1.0026</td><td>0.0130</td><td>93</td><td>5.8 ± 1.4</td><td>1.8 ± 0.2</td></tr><tr><td>1.0251</td><td>0.0096</td><td>45</td><td>5.6 ± 1.6</td><td>2.3 ± 0.3</td></tr><tr><td>1.1362</td><td>0.0113</td><td>114</td><td>4.4 ± 1.2</td><td>1.3 ± 0.2</td></tr><tr><td>1.1571</td><td>0.0096</td><td>119</td><td>5.0 ± 1.3</td><td>1.5 ± 0.2</td></tr><tr><td>1.2683</td><td>0.0113</td><td>43</td><td>4.2 ± 1.4</td><td>1.8 ± 0.3</td></tr><tr><td>1.2891</td><td>0.0095</td><td>80</td><td>5.6 ± 1.4</td><td>1.8 ± 0.2</td></tr><tr><td>1.4003</td><td>0.0114</td><td>125</td><td>8.1 ± 1.8</td><td>1.6 ± 0.2</td></tr><tr><td>1.4211</td><td>0.0095</td><td>60</td><td>4.0 ± 1.2</td><td>1.5 ± 0.2</td></tr><tr><td>1.5304</td><td>0.0095</td><td>61</td><td>6.7 ± 2.2</td><td>1.8 ± 0.3</td></tr><tr><td>1.5443</td><td>0.0044</td><td>62</td><td>4.2 ± 1.6</td><td>1.5 ± 0.3</td></tr><tr><td>1.6632</td><td>0.0104</td><td>48</td><td>5.6 ± 1.8</td><td>1.8 ± 0.3</td></tr><tr><td>1.6788</td><td>0.0052</td><td>103</td><td>6.5 ± 1.6</td><td>1.6 ± 0.2</td></tr><tr><td>2.1970</td><td>0.0095</td><td>73</td><td>6.9 ± 1.7</td><td>2.0 ± 0.2</td></tr><tr><td>2.2161</td><td>0.0096</td><td>75</td><td>5.8 ± 1.8</td><td>1.8 ± 0.3</td></tr><tr><td>2.5894</td><td>0.0096</td><td>50</td><td>6.7 ± 2.1</td><td>1.8 ± 0.3</td></tr><tr><td>2.6033</td><td>0.0043</td><td>51</td><td>7.1 ± 2.1</td><td>2.3 ± 0.4</td></tr></tbody></table>

Note. Days are relative to $t_0$, and $\Delta T$ represents the half bin width.

## Appendix B
Optical Data

### B.1. Optical Monitoring Data

We utilized visual-based measurements and optical photometry collected by the AAVSO for comparison to the LAT γ-ray light curves (Figure 2). The optical nova discoveries from days 0.41–0.43 (2021 August 8.91–8.93; Section 1) were at visual magnitudes of ~4.9–5.0. Two earlier visual observations

at ~5.1–5.2 mag from day 0.37 were later reported by J. Manzorro and A. Kosa-Kiss, and subsequent monitoring helped to detail the overall rising trend in the optical flux up to the peak on day ~1.1.

Pre-discovery photometry data were obtained serendipitously with commercial cameras and reported to the AAVSO. The earliest observations showing evidence of brightening were from four DSLR images obtained by Wang (2021; AAVSO Observer Code WBIA, including observers L.P. Lou and D.Y. Chen) indicating the initial rapid rise in brightness from $V=9.1$ to $7.0$ mag over a ~100 minute span from day 0.075 to 0.144; a limit of $V > 9.1$ was obtained from an observation on day 0.033. From a linear extrapolation of these observations back to the quiescent brightness of $V=11.1$ mag, Munari & Valisa (2021) derived a time of eruption of 2021 August 8.50 ($\pm 0.01$), which we adopt here as the reference epoch ($t_0$) in the LAT analysis and discussion. For the three weeks before the eruption, the quiescent brightnesses in the AAVSO database typically ranged from 11.1–11.3 mag, so the variations have only a small effect on the eruption epoch constraint.

The nova was also captured by the Global Meteor Network (Vida et al. 2021) camera IL0003 located in Israel (by A.L. and A.B.) with a 4 mm f/0.95 lens and a Sony IMX291 sensor. The 58 individual measurements from August 8, 17:35:46 to 21:09:05 (day 0.233 to 0.381) detailed further brightening<fnref n="27" /> following the first four observations from Wang (2021) up to the discovery observations. The individual measurements for August 8 as well as further observations on August 9 (37 measurements from 17:26:19 to 21:00:29; day 1.227 to 1.375) are newly presented in this paper (see Figure 8 in particular). Aperture photometry with $V$ zero-point on these unfiltered images (labeled CV in the AAVSO) was performed using the VaST code (Sokolovsky & Lebedev 2018) assuming $V=3.34$ for $\nu$ Oph (Ducati 2002), and these are simply referred to as $V$-band measurements throughout this paper. The second night measurements overlapping in time with the visual observations from the AAVSO are in good agreement considering the large scatter (~0.5 mag) in the meteor-camera observations due to the noisy and variable background structure visible across the CMOS images.

### B.2. Optical Luminosities From Observed Magnitudes

We also used available photometry collected by the AAVSO for RS Oph to calculate observed optical luminosities by assuming a Planck function spectral form with gray opacity to fit the photometric data at different epochs. Magnitudes were corrected for extinction due to the interstellar medium (ISM), and we applied approximate blackbody temperature-dependent bolometric corrections (Weidemann & Bues 1967). The temperatures were estimated with two methods, similar to those used by Li et al. (2017) and Aydi et al. (2020) to analyze the time variation in the ratio of γ-ray to optical luminosities of classical novae. Our aim was to use the fitted temperature as a parameter to estimate total optical energy fluxes, rather than deriving temperatures of the nova pseudo-photosphere.

<footnote n="27">²⁷ These are shown in place of the four averaged measurements reported by one of us (K.V.S.) to the AAVSO.</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Selected panels from Figure 3— LAT >0.1 GeV light curve in 10 minute bins in units of relative flux and Bayesian Block partitions indicated in red, LAT >0.1 GeV orbit bin light curve, and optical light curve, but for the first three days of activity. This version of the figure helps show more detail in the optical observations obtained from the Global Meteor Network camera IL0003 on the nights of August 8 and 9.</figcaption></figure>

The first method involves fitting a functional blackbody temperature at several dates to the SED composed of observed B-, V-, R-, and I-band magnitudes available in the AAVSO,<fnref n="28" /> after correction of extinction. The extinction correction for each band was calculated with a relative ISM extinction model from optical to mid-infrared bands obtained by Wang & Chen (2019) with photometric data collected by several instruments. This method requires quasi-simultaneous measurements in each of the spectral bands. Such measurements were performed at only six dates from day 4 to day 34. For all the observation periods, except the first one, the best-fit temperatures are ≈8000 K, with uncertainties of ≈1000 K. The quality of the SED fit in the first period (day 4) was poor; thus the best-fit temperature was not used. Instead, we adopted a temperature of 8500 K, which is intermediate between the temperatures of ~7000–10000 K obtained by Skopal (2015) from modeling the multiwavelength SEDs measured during the first days of the 2006 outburst of RS Oph (see Figure 9).

The second method is to calculate color temperatures, $T_c$, using quasi-simultaneous V- and B-band magnitude measurements and two-color indices corrected for extinction, as done by Li et al. (2017). Extinction was derived from the reddening $E(B-V) = 0.73 \pm 0.06$ (Snijders 1987) and the extinction ratio, $R_V = 3.1$. The V- and B-band measurements were obtained at 35 epochs from days ~-12 to ~55 (from AAVSO observers WGR, BDG, and FJQ). This method provides a range of $T_c = 8000$ to 14,000 K after the outburst date with an average of $T_c \approx (11,000 \pm 1700)$ K. Before the outburst, the color temperature was $T_c = (5900 \pm 400)$ K.

Figure 9 shows the obtained variation in the observed optical luminosity of nova RS Oph 2021 before and after its outburst obtained from daily averaged V- and B-band measurements

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Optical luminosities for the 2021 pre- and post-outburst emission from RS Oph derived from AAVSO data (red circles; see text). The best-fit power-law model for the decline (black line) has a slope of -1.28 ± 0.04. The optical luminosities estimated by Skopal (2015) for times near the peak of the previous 2006 outburst are shown for comparison (blue triangles; converted to our definition of $t_0$ by adding 1.1 day).</figcaption></figure>

from the AAVSO database. The values close to the luminosity peak are similar to the ones obtained by Skopal (2015) for the first days after the 2006 outburst. The luminosity appears to vary as a power law in time from day 2 onward with a best-fit slope of -1.28 ± 0.04.

## Appendix C: Swift X-ray Monitoring Data

Swift began observing RS Oph on 2021 August 9 (Page et al. 2021), continuing until the season ended on November 4,

<footnote n="28">²⁸ The approximate wavelengths for the filters are listed at https://www.aavso.org/filters.</footnote>
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Shock velocities for the nova ejecta from the 2021 outburst of RS Oph derived from Swift-XRT temperature measurements (red circles; see text). The dashed line corresponds to the profile of the 2021 outburst that we adopted in our modeling, consisting of a constant value of 2470 km s⁻¹ for *t* < 6 days and being proportional to *t*⁻⁰.⁴³ for *t* > 6 day (see Section 4). The velocities derived for the previous 2006 outburst are shown for comparison (blue triangles; converted to our definition of *t*₀ by adding 1.1 day).</figcaption></figure>

after which RS Oph became Sun constrained (observable again on 2022 February 1). Early observations were obtained on a daily basis; the cadence was increased to twice a day from the start of September, and then to every 8 hr between September 15 and October 1. Daily cadence then resumed until the observations ended. In addition, a very-high-cadence observing campaign (every Swift 1.5 hr orbit) was performed on September 12 to investigate potential supersoft source variability. All of these observations were taken using the windowed timing (WT) mode, as the XRT count rate always exceeded 1 count s⁻¹. Despite the high count rate, however, there were a number of separate observations between September 24 and October 8 performed in the photon counting (PC) mode specifically to investigate the point-spread function of the source emission; only the WT data are discussed in this paper.

The data were analyzed using HEASoft 6.29 and the latest calibration available in 2021 November. The light curve was produced using the online XRT product generator<sup>29</sup> (Evans et al. 2007, 2009), using only grade 0 (single pixel) events to help minimize optical loading.<sup>30</sup> Relevant to the comparison to the γ-ray emission, we used the higher-energy rates at 2–10 keV energies (Figures 2 and 3). In total, we found 121 significant detections in this band spanning days 1.4 to 87.6 from typical exposure times of 0.5–1 ks. The full XRT calibration and data set, including the detailed discussion of the supersoft emission around day 27 (using the *t*₀ defined here; Page 2021a, 2021b), are presented elsewhere (Page et al. 2022).

In our modeling of the 2021 outburst presented in Section 4, we used temperatures measured from spectral analysis of

Swift-XRT data to estimate the time variation in the shock velocity of the nova ejecta (see Equation 3 of Bode et al. (2006)). The results of the spectral fits of the 0.3–10 keV XRT data were taken from Page et al. (2022), who adopted a model using two optically thin plasma components to represent the shock emission, with an additional blackbody component at the lowest energies starting at day ~20. To estimate the velocity, we used the temperature of the hotter of the two shock-related components ($kT_{\text{hot}}$ in Page et al. 2022), which dominates the emission at >2 keV energies (see Figures 9 and 11, therein). The resultant velocities shown in Figure 10 only take the statistical uncertainties of the fitted temperatures into consideration; we do not account for systematic uncertainties related to the fits. Page et al. (2022) also reanalyzed the Swift-XRT data for the 2006 outburst in the same way, and we used those fitted $kT_{\text{hot}}$ values to derive the velocity data (only values of $kT_{\text{hot}}$ constrained to a factor of 2 were used). The obtained velocity profiles for the two explosions based on the Swift-XRT data are similar (Figure 10). Moreover, the XRT-derived velocity profile for the 2006 data is similar to the one used by Tatischeff & Hernanz (2007; Figure 1 therein) to model that explosion based on the early Swift-XRT results of Bode et al. (2006) and the Rossi X-Ray Timing Explorer (RXTE) results from Sokoloski et al. (2006).

## ORCID iDs

C. C. Cheung https://orcid.org/0000-0002-4377-0174

T. J. Johnson https://orcid.org/0000-0002-2771-472X

P. Jean https://orcid.org/0000-0002-1757-9560

M. Kerr https://orcid.org/0000-0002-0893-4073

K. L. Page https://orcid.org/0000-0001-5624-2613

J. P. Osborne https://orcid.org/0000-0002-1041-7542

A. P. Beardmore https://orcid.org/0000-0001-5186-5950

K. V. Sokolovsky https://orcid.org/0000-0001-5991-6863

F. Teyssier https://orcid.org/0000-0002-9107-9791

S. Ciprini https://orcid.org/0000-0002-0712-2479

G. Martí-Devesa https://orcid.org/0000-0003-0766-6473

I. Mereu https://orcid.org/0000-0003-0219-4534

S. Razzaque https://orcid.org/0000-0002-0130-2460

K. S. Wood https://orcid.org/0000-0002-7376-3151

S. N. Shore https://orcid.org/0000-0003-1677-8004

## References

Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, *Sci*, 329, 817

Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, *ApJS*, 247, 33

Abdollahi, S., Acero, F., Baldini, L., et al. 2022, *ApJS*, 260, 53

Ackermann, M., Ajello, M., Albert, A., et al. 2014, *Sci*, 345, 554

Ahnen, M. L., Ansoldi, S., Antonelli, L. A., et al. 2015, *A&A*, 582, A67

Aliu, E., Archambault, S., Arlen, T., et al. 2012, *ApJ*, 754, 77

Atwood, W., Albert, A., Baldini, L., et al. 2013, *arXiv:1303.3514*

Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, *ApJ*, 697, 1071

Aydi, E., Sokolovsky, K. V., Chomiuk, L., et al. 2020, *NatAs*, 4, 776

Barry, R. K., Mukai, K., Sokoloski, J. L., et al. 2008, in *ASP Conf. Ser.* 401, RS Ophiuchi (2006) and the Recurrent Nova Phenomenon, ed. A. Evans et al. (San Francisco, CA: ASP), 52

Beck, S. 2021, *AAVSO Alert Notice* 752: Rare Outburst of Recurrent Nova RS Ophiuchi (Cambridge: AAVSO), https://www.aavso.org/aaenso-alert-notice-752

Bode, M. F. 1987, in Proc. of the Manchester Conf. 241, RS Ophiuchi (1985) and the Recurrent Nova Phenomenon, ed. M. F. Bode (Utrecht: VNU Science Press), 241

Bode, M. F., O’Brien, T. J., Osborne, J. P., et al. 2006, *ApJ*, 652, 629

Brandi, E., Quiroga, C., Mikołajewska, J., et al. 2009, *A&A*, 497, 815

Bruel, P., Burnett, T. H., Digel, S. W., et al. 2018, *arXiv:1810.11394*

Burrows, D. N., Hill, J. E., Nousek, J. A., et al. 2005, *SSRv*, 120, 165

<footnote n="29">29 https://www.swift.ac.uk/user_objects/</footnote>

<footnote n="30">30 See https://www.swift.ac.uk/analysis/xrt/optical_loading.php.</footnote>
Buson, S., Jean, P., & Cheung, C. C. 2019, ATel, 13114, 1

Cheung, C. C., Ciprinini, S., & Johnson, T. J. 2021a, ATel, 14834, 1

Cheung, C. C., Donato, D., Wallace, E., et al. 2010, ATel, 2487, 1

Cheung, C. C., Jean, P., & Shore, S. N. 2014, ATel, 5879, 1

Cheung, C. C., Jean, P., Shore, S. N., et al. 2016, ApJ, 826, 142

Cheung, C. C., Johnson, T. J., Mereu, I., et al. 2021b, ATel, 14845, 1

Chomiuk, L., Metzger, B. D., & Shen, K. J. 2021, ARA&A, 59, 391

Dermer, C. D. 1986, ApJ, 307, 47

Dermer, C. D., & Powale, G. 2013, A&A, 553, A34

Dobrzycka, D., & Kenyon, S. J. 1994, AJ, 108, 2259

Ducati, J. R. 2002, yCat, 2227, 0

Evans, A., Bode, M. F., O'Brien, T. J., Darnley, M. J., et al. 2008, in ASP Conf. Ser. 401, RS Ophiuchi (2006) and the Recurrent Nova Phenomenon, ed. A. Evans et al. (San Francisco, CA: ASP), 313

Evans, P. A., Beardmore, A. P., Page, K. L., et al. 2007, A&A, 469, 379

Evans, P. A., Beardmore, A. P., Page, K. L., et al. 2009, MNRAS, 397, 1177

Fekel, F. C., Joyce, R. R., Hinkle, K. H., et al. 2000, AJ, 119, 1375

Fermi Science Support Development Team 2019, FermiTools: Fermi Science Tools, Astrophysics Source Code Library, ascl:1905.011

Franckowiak, A., Jean, P., Wood, M., et al. 2018, A&A, 609, A120

Geary, K. 2021, VSNET-alert notice, 26131, 1, http://ooruri.kusastro.kyoto-u.ac.jp/mailarchive/vsnet-alert/26131

Gehrels, N., Chincarini, G., Giommi, P., et al. 2004, ApJ, 611, 1005

Hernanz, M., & José, J. 2008, in ASP Conf. Ser. 401, RS Ophiuchi (2006) and the Recurrent Nova Phenomenon, ed. A. Evans (San Francisco, CA: ASP), 313

Hernanz, M., & Tatischeff, V. 2012, BaltA, 21, 62

H.E.S.S. Collaboration, Aharonian, F., Ait Benkhali, F., et al. 2022, Sci, 376, 77

Hjellming, R. M., van Gorkom, J. H., Taylor, A. R., et al. 1986, ApJL, 305, L71

Jean, P., Cheung, C. C., Ojha, R., et al. 2018, ATel, 11546, 1

Kamae, T., Karlsson, N., Mizuno, T., Abe, T., & Koi, T. 2006, ApJ, 647, 692

Kerr, M. 2019, ApJ, 885, 92

Li, K.-L., Metzger, B. D., Chomiuk, L., et al. 2017, NatAs, 1, 697

Linford, J. D., Chomiuk, L., Nelson, T., et al. 2017, ApJ, 842, 73

López-Coto, R., Sitarek, J., Bednarek, W., et al. 2015, ICRC (The Hague), 34, 731

Luna, G. J. M., Sokoloski, J. L., Mukai, K., et al. 2020, ApJL, 902, L14

MAGIC Collaboration, Acciari, V. A., Ansoldi, S., Antonelli, L. A., et al. 2022, NatAs, 6, 689

Mannheim, K., & Schlickeiser, R. 1994, A&A, 286, 983

Martin, P., Dubus, G., Jean, P., et al. 2018, A&A, 612, A38

Mattox, J. R., Bertsch, D. L., Chiang, J., et al. 1996, ApJ, 461, 396

Metzger, B. D., Finzell, T., Vurm, I., et al. 2015, MNRAS, 450, 2739

Montez, R., Luna, G. J. M., Mukai, K., et al. 2022, ApJ, 926, 100

Munari, U., & Valisa, P. 2021, arXiv:2109.01101

Narumi, H., Hirosawa, K., Kanai, K., et al. 2006, IAU Circ., 8671, 2

O'Brien, T. J., Bode, M. F., Porcas, R. W., et al. 2006, Natur, 442, 279

Oppenheimer, B. D., & Mattei, J. A. 1993, JAAVSO, 22, 105

Orlando, S., Drake, J. J., & Laming, J. M. 2009, A&A, 493, 1049

Osborne, J. P., Page, K. L., Beardmore, A. P., et al. 2011, ApJ, 727, 124

Page, K. L. 2021a, ATel, 14885, 1

Page, K. L. 2021b, ATel, 14894, 1

Page, K. L., Beardmore, A. P., Osborne, J. P., et al. 2022, MNRAS, 514, 1557

Page, K. L., Osborne, J. P., & Aydi, E. 2021, ATel, 14848, 1

Pickering, E. C. 1905, ApJ, 22, 90

Rupen, M. P., Mioduszewski, A. J., & Sokoloski, J. L. 2008, ApJ, 688, 559

Scargle, J. D., Norris, J. P., Jackson, B., et al. 2013, ApJ, 764, 167

Schaefer, B. E. 2004, IAU Circ., 8396, 2

Schaefer, B. E. 2009, ApJ, 697, 721

Schaefer, B. E. 2010, ApJS, 187, 275

Schaefer, B. E. 2018, MNRAS, 481, 3033

Schaefer, B. E. 2019, AAS Meeting, 234, 122.07

Skopal, A. 2015, NewA, 36, 128

Snijders, M. A. J. 1987, Ap&SS, 130, 243

Sokoloski, J. L., Luna, G. J. M., Mukai, K., et al. 2006, Natur, 442, 276

Sokoloski, J. L., Rupen, M. P., & Mioduszewski, A. J. 2008, ApJL, 685, L137

Sokolovsky, K. V., & Lebedev, A. A. 2018, A&C, 22, 28

Srivastava, M. K., Ashok, N. M., Banerjee, D. P. K., et al. 2015, MNRAS, 454, 1297

Tatischeff, V., & Hernanz, M. 2007, ApJL, 663, L101

Tatischeff, V., & Hernanz, M. 2008, in ASP Conf. Ser. 401, RS Ophiuchi (2006) and the Recurrent Nova Phenomenon, ed. A. Evans et al. (San Francisco, CA: ASP), 328

Tomsick, J. A., Boggs, S. E., Zoglauer, A., et al. 2022, ICRC (Berlin), 395, 652

Vida, D., Šegon, D., Gural, P. S., et al. 2021, MNRAS, 506, 5046

Vurm, I., & Metzger, B. D. 2018, ApJ, 852, 62

Walder, R., Folini, D., & Shore, S. N. 2008, A&A, 484, L9

Wang, B. 2021, Sky and Telescope, https://skyandtelescope.org/astronomy-news/recurrent-nova-rs-ophiuchi-just-blew-its-top/

Wang, S., & Chen, X. 2019, ApJ, 877, 116

Weidemann, V., & Bues, I. 1967, ZAp, 67, 415