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

We report the Fermi-LAT gamma-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 sub-intervals, are consistent with $\pi^{\rm 0}$-decay emission from shocks in the ejecta as proposed by Tatischeff&Hernanz (2007) for its previous 2006 outburst. The LAT light-curve displayed a fast rise to its peak>0.1 GeV flux of $\simeq$6x10^-6 ph cm^-2 s^-1 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 AAVSO, we measured a constant ratio of ~2.8x10^-3 between the gamma-ray and optical luminosities except for a ~5x smaller ratio within the first day of the eruption likely indicating attenuation of gamma 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.


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 years 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 (2009Schaefer ( , 2018, Montez et al. (2022), andMAGIC Collaboration et al. (2022).
Its 2006 outburst was well studied (Evans et al. 2008) with observations of hard X-ray emission from 2-25 keV (Sokoloski et al. 2006) and 14-50 keV ) that indicated shocked emission in the nova ejecta, and highresolution radio imaging resolved the shocked regions 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 1 , 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 RG wind model and helped to solidify the predictions of GeV emission in RS Oph specifically (Hernanz & Tatischeff 2012).
Thus the next outburst of RS Oph was highly anticipated by multi-wavelength 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 years, close to its average recurrence interval of 14.7 years (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) 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 π 0 -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 0 = 2021 August 8.5 (Munari & Valisa 2021, and § B.1) are used to describe the RS Oph outburst.

Fermi -LAT OBSERVATIONS AND DATA ANALYSIS
For the LAT analysis, we used Pass 8 (P8R3) SOURCE class data 3 (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 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 was 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. 2020Abdollahi et al. , 2022 based on twelve 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 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. component (iso P8R3 SOURCE V3 v1.txt) 4 . 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 ) −Γ , and fitted the normalization (N 0 ) and photon index (Γ), 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 ≥ 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 one-year dataset 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 × 10 −8 ph cm −2 s −1 (<1.5 × 10 −8 ph cm −2 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 one-year analysis. For point sources with free parameters found to have TS ≥ 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) where a is the fitted exponential factor 5 . 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) ≥ 1 × 10 −6 ph cm −2 s −1 based on preliminary results (Cheung et al. 2021b) 6 . We 4 Both files are available for download, https://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html 5 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 Ecut = a −1 . 6 Following the initial LAT detection of RS Oph on August 8.75-9.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 dataset used to fit the background, and a single PL spectral model was adopted for RS Oph throughout.
performed a binned maximum likelihood analysis on this 10-day period over the 0.05-300 GeV energy range, similar to what was done for the 1-year 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 TS curve = −2(L ECPL − L PL ) as a measure of significance of curvature in the spectrum, where 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), Γ = 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-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 § 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).

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 § 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 ( § 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 dataset described in § 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 ( § 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 hrs (see § A.2 for details).
Composite results from the different time-binned results are shown Figure 2 and Figure 3. Throughout, we consider significant detections at TS ≥ 12 (≥3σ significance for two degrees of freedom, d.o.f.), 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 4 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 § 3 and § 4, respectively.

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 ∼ √ TS for 1 d.o.f.). 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.
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 γ (>0.1 GeV) = (5.8 ± 1.2) × 10 −6 ph cm −2 s −1 , while the adjacent orbit-bin fluxes are (3.4 ± 1.0) × 10 −6 ph cm −2 s −1 , 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 π(N p ) −γ with γ = 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σ.

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 from the optical position of RS Oph. The resultant list of photons is given in the Supplement ( § 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 dataset described in § 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 § 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 ( § 2.1), a single 12.8 GeV photon with smaller gtsrcprob weight value = 0.85 was found at day 55.

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 γ = (2.1 ± 0.7) × 10 −6 ph cm −2 s −1 ; (b) the source rises to a peak F γ = (5.7 ± 1.2) × 10 −6 ph cm −2 s −1 at day 0.745; (c) the peak is flat, with an approximately constant (slope of −0.18 ± 0.17), and average F γ 5 × 10 −6 ph cm −2 s −1 through to day 2.5; (d) the flux declined by a factor of two 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 ± 0.11 fit from day 2.5 to 45, and the average peak flux from day 0.745-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 highestenergy (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-1.940, and the subsequent flux decline. The fitted LAT spectra have a wide range of spectral slopes, Γ = 1.5 to 2.1, with typical errors of 0.2. In the rising portion of the light-curve, the γ-ray spectral slopes were Γ ∼ 2.0-2.3, while the spectra hardened (Γ ∼ 1.7) during the peak. The hardest spectrum (Γ = 1.4 ± 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 ∼3month 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 ( § 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 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 ± 0.69) is consistent with the optical V -band one (1.28 ± 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-1.3 (Visual = 4.5-4.6 mag), within the wider timespan of the observed γ-ray peak fluxes from day 0.75-1.67 (see Table 3). Fitting a broken power-law to the Visual measurements from days 0.37-3.0, the best-fit peak was at day 1.09 (±0.12), consistent with the day 1.08 (±0.05) estimated by Munari & Valisa (2021). Although the γ-ray peak fluxes appear constant over a ∼1-day interval, assuming a broken 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 7 as a function of the observed optical luminosity estimated in 4-day bins using the AAVSO data (see § B.2). Both the 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 MeV)/L opt = (2.81 ± 0.15) × 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 § 4 for a discussion), while the other ratio values are compatible with the constant value.
7 The γ-ray luminosities were calculated by fitting the fluxes and the proton spectrum slopes (without energy cutoff) in the π 0 model, to the Fermi-LAT data (see § 4).
Optical luminosity (erg s 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 § 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 § 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 ( § 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.

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 (cf., 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 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 distributions (SEDs; Figure 5). The corresponding best-fit parameters obtained for each Energy (MeV) observation period are presented in Table 1 along with the difference between the TS values resulting from the best-fit hadronic and ECPL models (∆TS). According to the ∆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 (Γ = 1.66 ± 0.05, E cut = 6.0 ± 1.2 GeV; § 2), suggests a modest spectral change in the fit for the subset of data for the peak period (E cut = 3.16 +1.10 −0.07 GeV). The slopes of the proton spectra are compatible with a constant value of ∼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 π 0 -decay produced by high-energy protons interacting with the material of the ejecta as proposed by Tatischeff & Hernanz (2007) to 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 of 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 −1 for t < 6 days and proportional to t −0.43 at t > 6 day (which is close to the t −0.5 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 4 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 to 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 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 π 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 −6 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 inj ∼ 4 × 10 −6 , lower than the value of 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 8 (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 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.

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 γ (>0.1 GeV) 6 × 10 −6 ph cm −2 s −1 , 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 γ (>0.1 GeV) 4 × 10 −6 ph cm −2 s −1 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 ( § 2. 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 two fluctuations in the flux on timescales of 30 to 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 of 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 γ (>0.1 GeV) 6 × 10 −6 ph cm −2 s −1 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 year −1 starting in 2012 (Ackermann et al. 2014;Cheung et al. 2016); see Chomiuk 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). 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 hadronic processes from the external shock discussed in 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 2010Luna et al. 2020), and if the distance-scaled fluxes hold, should be remarkably bright with fluxes, F γ (>0.1 GeV) > ∼ 1 × 10 −5 ph cm −2 s −1 thus 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 nine years (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 γ-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 a 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.  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. The full results of the LAT >0.1 GeV light-curve spectral analysis in 6-hr, 1-day, and 4-day bins described in § 2.1 are presented in Figure 7.
The list of LAT photons with E > 5 GeV found in the analysis described in § 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 (θ) between the center of our ROI and the LAT z-axis and exclude from our final analysis those orbits where θ 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 fermitool gtltcube and then performed an unbinned likelihood analysis, with only the isotropic diffuse emission component and the normalization and Γ 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 five minute bins, made using the fermitools gtbin and gtexposure. All detections were TS > 25 (>4σ for 2 d.o.f.) 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 to 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 which yields equal exposure. The results are presented in Table 4. The most-significant variability was a factor of two drop at day 1.41 as described in Section 2.2. For a single orbit (days 0.745), the photon index Γ appears to increase by 1.0 ± 0.5 (∼2σ).
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; § 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 (±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 9 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 zeropoint on these unfiltered images (labelled CV in the AAVSO) was performed using the VaST code (Sokolovsky & Lebedev 2018) assuming V = 3.34 for ν 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 temperaturedependent 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 of 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.
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 10 , 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 multi-wavelength 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 ± 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 14000 K after the outburst date with an average of T c ≈ (11000 ± 1700) K. Before the outburst, the color temperature was T c = (5900 ± 400) K. Figure 9 shows the thus-obtained variation of the observed optical luminosity of nova RS Oph 2021 before and after its outburst obtained from daily-averaged V -and B-band measurements 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.

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, 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 Windowed Timing (WT) mode, since the XRT count rate always exceeded 1 count s −1 . Despite the high count rate, however, there were a number of separate observations between September 24 and October 8 performed in 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 11 (Evans et al. 2007(Evans et al. , 2009), using only grade 0 (single pixel) events to help minimize optical loading 12 . Relevant to the comparison to the γ-ray emission, we used the higher-energy rates at 2-10 keV energies (Figures 2 & 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 dataset, including the detailed discussion of the supersoft emission around day 27 (using the t 0 defined here; Page 2021a,b), is presented elsewhere (Page et al. 2022).
In our modeling of the 2021 outburst presented in § 4, we used temperatures measured from spectral analysis of Swift-XRT data to estimate the time-variation of 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 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 hot values to derive velocity data (only values of kT 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,  . 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 t0 by adding 1.1 day).
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).