Evidence for multiple shocks from the $\gamma$-ray emission of RS Ophiuchi

In August of 2021, Fermi-LAT, H.E.S.S., and MAGIC detected GeV and TeV $\gamma$-ray emission from an outburst of recurrent nova RS Ophiuchi. This detection represents the first very high energy $\gamma$-rays observed from a nova, and opens a new window to study particle acceleration. Both H.E.S.S. and MAGIC described the observed $\gamma$-rays as arising from a single, external shock. In this paper, we perform detailed, multi-zone modeling of RS Ophiuchi's 2021 outburst including a self-consistent prescription for particle acceleration and magnetic field amplification. We demonstrate that, contrary to previous work, a single shock cannot simultaneously explain RS Ophiuchi's GeV and TeV emission, particularly the spectral shape and distinct light curve peaks. Instead, we put forward a model involving multiple shocks that reproduces the observed $\gamma$-ray spectrum and temporal evolution. The simultaneous appearance of multiple distinct velocity components in the nova optical spectrum over the first several days of the outburst supports the presence of distinct shocks, which may arise either from the strong latitudinal dependence of the density of the external circumbinary medium (e.g., in the binary equatorial plane versus the poles) or due to internal collisions within the white dwarf ejecta (as powers the $\gamma$-ray emission in classical novae).


INTRODUCTION
A major discovery by the Fermi Large Area Telescope (LAT) was that novae−multi-wavelength transients produced by non-terminal thermonuclear explosions on the surface of white dwarfs accreting hydrogen-rich material from a donor star−are sources of luminous ∼GeV γ-ray emission (e.g., Abdo et al. 2010;Ackermann et al. 2014;Cheung et al. 2016;Franckowiak et al. 2018). Novae come in two varieties: "classical" novae, when the donor is a main-sequence or moderately-evolved star overflowing its Roche Lobe onto the white dwarf; and "embedded" or "symbiotic" novae, when the donor is instead a giant star with a dense wind.
The first nova detected in γ-rays, V407 Cyg, was of the embedded type, with a Mira-like red giant donor (Abdo et al. 2010). The γ-rays from this event, which started around the time of optical maximum and lasted about two weeks, were interpreted by Abdo et al. (2010) as non-thermal emission from relativistic particles (ions or electrons) accelerated to high energies at the shock wave generated as the nova ejecta collided with the dense circumbinary environment of the giant companion (e.g., Chomiuk et al. 2012). Martin & Dubus (2013) modeled the diffusive acceleration of particles from the resulting shock wave, finding that inverse Compton emission by relativistic electrons interacting with the red giant optical light are the dominant γ-ray emission source; they also found that a "density enhancement" in the binary plane was needed to match the data in addition to the standard red giant wind (Booth et al. 2016 describe a possible origin for such an equatorial density enhance-ment generated by the white dwarf's accretion of the giant wind).
Given the relative rarity of novae with giant companions, γ-ray detections of novae were predicted to be rare (see the discussion in Abdo et al. 2010). However, this expectation was upended when Fermi LAT began to detect additional Galactic novae, even of the more common classical variety (Ackermann et al. 2014;Cheung et al. 2016;Franckowiak et al. 2018; a total of 15 LAT-detected classical novae to date). Given the lack of a dense wind from a main-sequence star, the shocks in classical novae are likely "internal", i.e., occurring between distinct components of nova ejecta. Indeed, multi-wavelength observations from X-ray (e.g., Mukai & Ishida 2001;Takei et al. 2009;Nelson et al. 2019;Gordon et al. 2021) to radio (e.g., Chomiuk et al. 2014;Weston et al. 2016), show internal shocks to be common if not ubiquitous in classical novae (see Chomiuk et al. 2021 for a recent review). Resolved radio imaging (e.g., Chomiuk et al. 2014) supports a picture in which the shock interaction is caused by a quasi-spherical highvelocity outflow or wind from the white dwarf which collides with a slower outflow concentrated in the plane of the binary and released earlier in the eruption. Additional evidence for this interpretation comes from timeresolved optical spectra, which show evidence for multiple velocity components which exist simultaneously, the faster of which first appears around optical maximum and near the onset of the γ-ray emission (e.g., Aydi et al. 2020a;. The gas surrounding the internal shocks in classical novae is sufficiently dense to act as a "calorimeter" for converting their thermal UV/X-ray emission (Metzger et al. 2014) and non-thermal accelerated particles (Vurm & Metzger 2018;Martin et al. 2018) into reprocessed optical and γ-ray emission, respectively. This should result in a contribution to the nova optical light curve from shocks that tracks the γ-ray light curve (Metzger et al. 2015), consistent with that observed in the few novae for which such a measurement has been possible (Li et al. 2017;Aydi et al. 2020b).
RS Ophiuchi (RS Oph) is a binary consisting of a red giant orbiting a white dwarf with a 454 day period (Brandi et al. 2009), which has undergone a nova eruption every ∼ 10 − 20 years for over the last century (Schaefer 2010). Its 2006 eruption occurred before the launch of Fermi LAT, but signatures of shock interaction were observed in impressive detail at X-ray, infrared, and radio wavelengths (Sokoloski et al. 2006;Bode et al. 2006;O'Brien et al. 2006;Das et al. 2006;Tatischeff & Hernanz 2007;Evans et al. 2007;Rupen et al. 2008). Its most recent eruption, beginning on 2021 Aug 8, qualified RS Ophiuchi as the second embedded nova detected at high-significance by Fermi LAT (Cheung et al. 2022), and also as the first nova detected at ∼ TeV energies, by the Atmospheric Cherenkov Telescopes H.E.S.S. (Aharonian et al. 2022) and MAGIC (Acciari et al. 2022). No classical nova has yet been detected at TeV energies, though relatively few upper limits have been reported in the literature (e.g., by MAGIC and HAWC; Ahnen et al. 2015;Albert et al. 2022) and there are theoretical reasons to believe that the high neutral gas densities surrounding the shock may limit the maximum particle energy (Reville et al. 2007;Metzger et al. 2016).
Interestingly, the TeV light curve peak in RS Ophiuchi is delayed by several days relative to the peak of the GeV light curve, with the Fermi LAT light curve peaking on 2021 Aug 9-10 (Cheung et al. 2022) and the HESS light curve peaking on 2021 Aug 12 (Aharonian et al. 2022; see their Figure 2). The HESS collaboration interprets the entire GeV-TeV emission in terms of hadronic particle acceleration and emission at a single shock, attributing the observed temporal delay between the peaks at different energies to the finite timescale required to accelerate ions to TeV energies. However, in this paper, we show that the γ-rays observed during RS Ophiuchi's 2021 outburst cannot be produced by a single, spherically symmetric shock. Instead, we put forward a scenario involving multiple shocks. These shocks may be generated as the result of distinct velocity components of the nova ejecta interacting with the aspherical external environment. This scenario can reproduce key features of the observed γ-rays without any ad-hoc modifications to the shape or maximum energy of the accelerated particle spectrum.
In Section 2 we describe our model for shock evolution, particle acceleration, and magnetic field amplification. We also introduce constraints on our model from optical spectroscopy. We apply our model to RS Ophiuchi in Section 3 and demonstrate how a single, spherical shock cannot reproduce the γ-ray observations. We show the results of our best-fit multi-shock model in Section 4, and discuss physical scenarios that might correspond to this model in Section 5.
Throughout this work, we assume a distance to RS Ophiuchi of 1.4 kpc (Bode et al. 2008), consistent with that used in Aharonian et al. (2022).

METHOD
To calculate RS Ophiuchi's expected γ-ray emission we use a multi-zone model of particle acceleration and photon production. A detailed description of this model can be found below.

Shock Hydrodynamics
To estimate shock evolution, we use a self-similar formalism similar to that described in Diesing & Caprioli (2018). Namely, we assume that both the material ejected by the nova and the material swept up during expansion are confined to a thin shell behind the shock (see, e.g., Bisnovatyi-Kogan & Silich 1995;Ostriker & McKee 1988;Bandiera & Petruk 2004, for examples of this thin-shell approximation). The evolution of the shock is thus set by the density profile of the ambient medium, chosen to reproduce the γ-ray observations and to be consistent with the presence of a RG wind (see Section 3).
More specifically, we model RS Ophiuchi during two stages of evolution: the ejecta-dominated stage, in which the mass of swept-up material is less than the ejecta mass and the shock expands freely, and the Sedov stage, in which the swept-up mass exceeds the ejecta mass and the nova expands adiabatically. Energy is conserved throughout both stages such that, given a nova with (1) Here, v sh is the forward shock velocity, M ej is the ejecta mass, and M SU is the swept up mass, given by, The above formalism applies to a single shock. However, as we shall discuss in later section, the γ-ray emission and optical spectra of RS Ophiuchi both suggest the presence of multiple shock components. These shocks may arise due to distinct mass ejection events from the white dwarf (as in the case of classical novae; see discussion in Sec. 1), and/or a single ejecta shell interacting with non-spherically symmetric external medium. In the latter case, we will in effect be applying the above formalism separately to distinct angular sectors (e.g., the polar versus equatorial region) over which we assume the density is approximately uniform.
The above model assumes a stationary external medium, which is a good approximation for the slowlyexpanding RG wind or circumbinary disk expected around RS Ophiuchi (e.g., Booth et al. 2016). However, it would less well approximate the dynamics of internal shocks between distinct ejecta components from the white dwarf (e.g., Metzger et al. 2014), as are believed to power the GeV γ-ray emission from classical novae (Sec. 1). Nevertheless, the above framework may still provide a rough description of the internal shock case, provided that the initial shock velocity is interpreted as the relative velocity between the ejecta components.
It is also worth noting that nova outbursts may also be described by a continuous wind as opposed to an instantaneous energy injection (e.g., Kato et al. 2022). In this case, the hydrodynamic model used in our work does not apply and, importantly, a reverse shock may contribute to the observed γ-ray luminosity. However, the fact that the optical emission begins to decay before the GeV peak suggests that such extended energy injection is unlikely (e.g., Metzger et al. 2014;Li et al. 2017). More specifically, after its maximum, the optical luminosity decays as t −1.3 , implying that a substantial portion of the optical energy is released near its peak. By the time of the GeV maximum, we do not expect significant ongoing energy injection from a wind, meaning that our instantaneous ejection model provides a good approximation for the hydrodynamics of the system.
Admittedly, it may still be possible for a reverse shock to contribute to the GeV emission at early times. However, such a contribution does not change the main conclusion of our work, that multiple shocks-be they forward shocks or forward and reverse shocks-are required to reproduce the observed gamma-ray emission.

Particle Acceleration
We model particle acceleration using a semi-analytic model of nonlinear diffusive shock acceleration that selfconsistently accounts for magnetic field amplification and the dynamical back-reaction of accelerated particles on the shock (see Caprioli et al. 2009Caprioli et al. , 2010aCaprioli 2012;Diesing & Caprioli 2019, and references therein, in particular Malkov (1997Malkov et al. (2000); Blasi (2002Blasi ( , 2004; Amato & Blasi (2005, 2006). We assume that protons with momenta above p inj ≡ ξ inj p th are injected into the acceleration process, where p th is the thermal momentum and we choose ξ inj to produce CR pressure fractions ∼ 10% (though we note that the non-thermal acceleration efficiency in classical novae is typically measured to be closer to ∼ 1%; e.g., Metzger et al. 2015;Aydi et al. 2020b). We also calculate the proton maximum energy self-consistently by requiring that the diffusion length (assuming Bohm diffusion) of the most energetic particles accelerated be equal to 5% of the shock radius.
This model produces an instantaneous distribution of protons accelerated at each timestep of nova evolution. Instantaneous electron distributions are calculated from these proton distributions following the analytical approximation in Zirakashvili & Aharonian (2007). We then shift and weight each instantaneous distribution to account for adiabatic and, in the case of electrons, synchrotron losses (see Caprioli et al. 2010b;Morlino & Caprioli 2012;Diesing & Caprioli 2019, for more details). We then sum these weighted distributions to yield the cumulative, multi-zone spectrum of non-thermal particles accelerated by our model novae.
We also account for proton-proton losses by calculating the collision rate for each instantaneous distribution (i.e., each expanding shell of protons) at each timestep, assuming the collisional cross-section parameterized in Kafexhiu et al. (2014) and a target proton density equal to the adiabatically expanded post-shock density of that shell. We further assume that a proton loses half its energy in a single collision (i.e., we assume an inelasticity κ = 0.5, consistent with Martin & Dubus 2013), and modify each instantaneous proton distribution accordingly.
Note that, unless the ambient density is quite large, we do not expect significant proton-proton losses on timescales of a few days 1 since, for an interaction cross section of ∼ 30 mb and κ = 0.5, the proton-proton loss time can be approximated as, where n 0 is the number density of the ambient medium in front of the shock.

Photon Production
To estimate photon spectra from our cumulative proton and electron distributions, we use the radiative processes code naima (Zabalza 2015). Naima computes the emission due to synchrotron, Bremsstrahlung, inverse Compton (IC) and neutral pion decay processes assuming arbitrary proton and electron distributions, as well as our chosen density profile(s). While the IC luminosity depends also on the radiation field chosen, we find that leptonic emission is subdominant regardless of our assumptions (see Section 3).
The main sources of opacity for GeV -TeV photons are pair production-on soft radiation fields in the RG wind and on nuclei in the nova ejecta. The opacity throughout the nova ejecta (Bethe-Heitler process) for hydrogen-rich material is given by (e.g. Zdziarski & Svensson 1989) where x ≡ E γ /m e c 2 , α fs 1/137 and σ T 6.6 × 10 −25 cm −2 is the Thomson cross section. Thus, the Bethe-Heitler optical depth through the RG wind, is negligible. The Bethe-Heitler optical depth is also likely to be irrelevant even if the density is enhanced by a factor ∼ 10 3 on small scales AU. γ-ray photons can also be attenuated due to γ −γ pair production with ambient IR/optical/UV photons, again either from the RG or from the nova outburst itself.
We calculate the optical depth at the location of the nova shock (r sh ) for γ-rays due to absorption on the IR/optical radiation field as, where n opt is the number density of target photons with energy t and σ γγ is the interaction cross section. The target radiation field is normalized to the bolometric optical luminosity, L opt given by Cheung et al. (2022) (their figure 9), The target radiation field is assumed to have a blackbody shape with T BB = 10 4 K. Our lack of knowledge of the actual broadband spectral shape is unlikely to introduce a significant error, since the primary targets for γ-rays of interest are within or close to the optical band (e.g., at 1 eV for 500 GeV γ-rays), hence the γ − γ opacity is less sensitive to the precise shape of the target spectrum outside the observed range than it would be e.g. for γ 1 TeV photons. We show the results of this calculation in Figure 1. r sh is taken to be that of our best-fit model involving a single, external shock, as described in Section 3. While we do expect modest attenuation of TeV γ-rays (by a factor of ∼ 2) at t 1 day, this attenuation is not sufficient to account for the order of magnitude rise in the TeV luminosity observed between days 1 and 4. Furthermore, since the γ − γ opacity is negligible at the radius corresponding to the TeV luminosity peak, we neglect absorption in our emission estimates.
It is also worth noting that eq. (6) employs the angleaveraged γ −γ cross-section, or equivalently, an isotropic target radiation field. If the optical radiation is produced at smaller radii than the VHE γ-rays, its beaming  Table 1). The target radiation is assumed to have blackbody shape of TBB = 10 4 K and a bolometric luminosity given by Cheung et al. (2022). The vertical black and magenta dashed lines indicate the time of Fermi-LAT and H.E.S.S. luminosity peaks, respectively.
around the radial direction tends to suppress the attenuation of radially-directed γ-rays. In the limit of perfectly collimated targets, an escape cone exists for the γ-rays for any target radiation density; the (angle-averaged) attenuation is no longer exponential and a fraction of γ-rays can escape even if the isotropically-averaged τ γγ would indicate essentially complete suppression.

Magnetic Field Amplification
The streaming of energetic particles ahead of the shock is expected to excite various instabilities (Skilling 1975;Bell 2004;Amato & Blasi 2009;Bykov et al. 2013), which amplify magnetic fields and enhance CR diffusion (Caprioli & Spitkovsky 2014a,b,c). This amplification has been inferred observationally in supernova remnants (e.g., Parizot et al. 2006;Bamba et al. 2005;Morlino et al. 2010;Ressler et al. 2014), and is expected to proceed in a similar manner in nova shocks.
We model magnetic field amplification as in Diesing & Caprioli (2021) by assuming saturation of both the resonant streaming instability (e.g., Kulsrud & Pearce 1968;Zweibel 1979;Skilling 1975;Bell 1978;Lagage & Cesarsky 1983a), and the non-resonant hybrid instability (Bell 2004). This prescription reproduces the magnetic fields inferred from X-ray observations of young supernova remnants (Völk et al. 2005;Caprioli et al. 2008). Bell (2004) derives the saturation point of the nonresonant instability to be Here, γ CR = 4/3 is the CR adiabatic index. This saturation has been validated with hybrid simulations in Zacharegkas et al. (2021Zacharegkas et al. ( , 2022. For fast shocks ( 100 km s −1 for typical nova parameters), the non-resonant instability dominates amplification (see Diesing & Caprioli 2021, for a detailed discussion), so we assume that the total magnetic pressure immediately upstream of the shock is then P B,1 P B1,res ; moving further upstream, this pressure is taken to scale with the local P CR (x). Thus, for a fast shock capable of accelerating TeV particles, we expect the magnetic field in front of the shock to scale as v 3/2 sh , assuming the CR pressure scales with the ram pressure, ρ 0 v 2 sh . Note that, given the strong magnetic field amplification in our model, the ambient magnetic field-a relatively unknown quantity-has a negligible impact on our results.
The analysis above assumes fully ionized gas surrounding the shock. However, this may not be justified for sufficiently dense upstream gas, because the recombination timescale may be longer than the ionization timescale by UV/X-ray radiation from the shock (indeed, absorption and reprocessing of the shock thermal emission is key to powering nova optical emission from the shocks; Metzger et al. 2014). Upstream gas with a substantial neutral component can suppress the growth rate of the non-resonant instability (Bell 2004;Reville et al. 2007), thus reducing the maximum energy of the ions accelerated at nova shocks (Metzger et al. 2016).

Constraints from Optical Spectroscopy
We now consider constraints on the properties of the ejecta from RS Ophiuchi based on optical spectroscopy.
We make use of publicly available high-resolution spectroscopy from the Astronomical Ring for Access to Spectroscopy (ARAS; Teyssier 2019) database, covering the first month of the eruption, starting from 0.5 days after t 0 (HJD 2459435.0745). We present the evolution of Hβ and Fe II 5169Å line profiles during the first 30 days in Figures 7, 8, 9, and 10. The line shows multiple absorption components during the first week of the eruption. We identify at least 3 components: (1) an initial component that is present in the first spectrum taken 0.5 days after eruption discovery, and is at a velocity v 1 = −2700 km s −1 ; (2) a faster component characterized with a velocity v 2 = −3700 km s −1 that appears a day later and co-exists with the initial component for a few days. This component is only resolved in the Balmer lines; (3) a third component characterized by a velocity v 3 = −1900 km s −1 and appearing around days 3 -4 in the Balmer lines, but which is prominent earlier in the Fe II lines (possibly as early as days 1 -2). All these components show variations in their velocity/strength as the eruption evolves. Since these components co-exist in the same lines, they possibly originate in distinct regions of the nova ejecta or CSM. Figure 2 shows the evolution of these components relative to the optical and high-energy light curves. One interpretation of these data is that the 2700 km s −1 component originates in a comparatively "slow" initial mass ejection phase concentrated in the binary equatorial plane, while the 3700 km s −1 component−which accelerates to 4000 km s −1 over a few days−originates in a faster, radiation driven, wind that potentially expands more freely in the polar direction, in a scenario similar to that suggested for the early spectral evolution of classical novae (e.g., Aydi et al. 2020a). These components have also been reported in Molaro et al. (2022) and were observed in different species. The slowest component with a velocity of 1900 km s −1 could then arise from swept-up shell formed by the nova ejecta interacting with the slowly expanding (≈ 50 km s −1 ) circumbinary material.
Another possible interpretation of these velocity components is elucidated in . In this case, the components correspond to three different regions surrounding a single shock. However, this interpretation implies a shock velocity of 1900 km s − 1, which is insufficient to accelerate TeV particles.
In what follows, we use the velocities of these ejecta components to motivate those assumed in our shock models. To keep matters simple, we will only consider systems with a maximum of two velocity components. However, we choose shock velocities that are reasonably consistent with those inferred from optical data. Our main objective here is to demonstrate that the optical observations of distinct ejecta components motivate the presence of multiple shocks which, as we shall argue, is needed to interpret the γ-ray data.

THE SINGLE-SHOCK SCENARIO
In this section we apply our model to RS Ophiuchi assuming its emission arises from a single, external shock. Broadly speaking, a viable model must reproduce three key features observed by Fermi, H.E.S.S., and MAGIC (Aharonian et al. 2022;Acciari et al. 2022;Cheung et al. 2022): 1. An initial rise in both the GeV and TeV luminosities. 2. An eventual decay in both the GeV and TeV luminosities that goes as t −α where α 1.3 − 1.4.

3.
A delay in the TeV luminosity peak with respect to the GeV luminosity peak of roughly 2-3 days.
It is worth noting that, while H.E.S.S. observes a clear delay between GeV and TeV peaks, as well as a rise in the TeV emission at early times, MAGIC does not (Acciari et al. 2022). This discrepancy may be due to the fact that MAGIC began observing slightly later; we will therefore take the H.E.S.S. results at face value for the remainder of our analysis. Regardless, as we will show, the shape of the combined GeV-TeV spectrum alone requires the presence of multiple shock components.
The first two items on the above list can be reproduced with hadronic emission arising from a single external shock expanding into a medium of constant density that transitions to a red giant (RG) wind profile with density ρ 0 (r) ∝ r −2 at large radii ( 3 AU). Figure 3 shows an example of shock evolution in such a density profile and the corresponding GeV and TeV light curves (more specifically, the light curves in the 0.06-500 GeV and 250-2500 GeV bands, consistent with the data displayed in Aharonian et al. 2022).  Table 1. Gray lines indicate the approximate time of transition from the ejecta-dominated to the Sedov-Taylor stage.
For the overlaid observational data, t0 corresponds to the peak of the optical light curve. The modeled light curves match the Fermi data quite well (bottom; blue line), but substantially overestimate the H.E.S.S. data (bottom; red line). This overestimation comes from the fact that the combined Fermi and H.E.S.S. data cannot be described as a single power law with an exponential cutoff, implying a more complex picture than a single, external shock.
Shock parameters, listed in Table 1, are chosen to be broadly consistent with observations and to provide a good fit to the Fermi (GeV) data. In particular, we adopt a RG wind velocity and mass loss rate that are comparable to (within a factor of 2 of) those inferred observationally (e.g., Iijima 2009, finds a wind velocity of 33 km s −1 and a mass loss rate of ∼ 10 −6 M yr −1 ). It is worth noting that we adopt a smaller ejecta mass than that inferred by Pandey et al. (2022). Since estimates of the ejecta mass are highly model-dependent, we prioritize a value that yields good agreement with Fermi observations (i.e., yields deceleration around the time of the Fermi luminosity peak). It is also worth noting that our ejecta mass is roughly consistent with those adopted elsewhere in the literature (e.g, Booth et al. 2016).
A light curve peak that occurs ∼ 1 day after the shock is launched corresponds to a constant density region ∼ 3 AU, roughly twice the orbital distance of the RS Ophiuchi system (e.g., Booth et al. 2016). The model GeV and TeV light curves predicted by this model are plotted in the bottom panel of Figure 3. The model produces a good fit to the observed GeV light curve. However, this model substantially overestimates the observed TeV emission, due to the fact that the outburst's observed spectrum is inconsistent with the theoreticallymotivated power-law with an exponential cutoff. A detailed discussion of this inconsistency and its implications can be found later in this section.

Parameter
Quantity (single shock) Mej (ejecta mass) 2 × 10 −7 M v sh,init (initial v sh ) 4500 km s −1 n0,init (initial ambient density) 7 × 10 8 cm −3 rcrit (homogeneous region size) 3 AU M wind (RG mass loss rate) 5 × 10 −7 M yr −1 v wind (RG wind velocity) 30 km s −1 Table 1. Model parameters used to fit the Fermi light curve, assuming a single, external shock. The resulting hydrodynamic evolution and light curves can be found in Figure 3. Note that a smoothing function is applied to the density profile, such that the density-related parameters shown here differ slightly from the ultimate profile that enters into our model. Without including a smoothing function, our model still reproduces key features of the GeV data, with the exception of a sharp jump in the light curve at the transition from the homogeneous region to the wind profile.
The luminosity rise, decay, and corresponding density profile requirements shown in Figure 3 are best understood in terms of simple scaling relations. Namely, the luminosity of accelerated protons scales with the energy flux across the shock multiplied by the area of the shock surface L p ∝ ρ 0 v 3 sh R 2 sh . For hadronic emission (i.e., pion decay), the γ-ray luminosity, L γ , scales with the proton luminosity, the target density, and the shock age: L γ ∝ ρ 2 0 v 2 sh R 3 sh . Assuming a power-law scaling of the ambient density, during the ejecta-dominated stage, when the shock velocity is roughly constant and R sh ∝ t. Thus, a rise in luminosity at early times is only possible if s < 3/2, suggesting that the material in the region closest to the white dwarf does not follow a wind profile (s = 2). Practically speaking, a flat profile (s = 0) at small radii reproduces the observed luminosity rise quite well. The actual CSM profile surrounding the white dwarf, as fed by the wind and Roche-lobe overflow of the giant, is expected to be complex (Martin & Dubus 2013;Booth et al. 2016). During the Sedov stage, R sh ∝ t 2/(5−s) yielding, Thus, if the shock enters the Sedov stage while the ambient density is still relatively constant, L γ is roughly constant. However, if the shock enters the Sedov stage around or after the ambient material has transitioned to a wind, L γ ∝ t −4/3 , consistent with the decays observed by both Fermi and H.E.S.S. Together, these scaling relations suggest that RS Ophiuchi's outburst expanded first into a roughly homogeneous medium on scales comparable to the binary orbital separation, followed by the ρ 0 ∝ r −2 spherical red-giant wind on larger radial scales. The fact that the observed GeV and TeV light curves rapidly transition from rising to power-law decay with α 0.3 − 0.4 suggests that the onset of the Sedov stage roughly coincides with the start of the wind. This coincidence allows us to place constraints on the extent and density of the homogeneous region. Namely, the total mass in this region must be approximately equal to the ejecta mass. This requirement informs the model shown in Figure 3.
To avoid sharp jumps in our modeled light curve, we use a smoothing function to connect the uniform and wind-profile regions of our models. In reality, the CSM distribution in this messy transition, from the "disklike" CSM near the white dwarf to the RG wind, is likely to be complex and not spherically symmetric (Booth et al. 2016); our adopted smoothing thus also crudely accounts for variations in the shock evolution across different solid angle sectors.
While the single-shock model can describe the overall shapes of both the GeV and TeV light curves, it is inconsistent with observations in two key ways. First, and most obviously, it overestimates the very high energy (VHE) γ-ray flux by more than an order of magnitude. This arises from the fact that the combined Fermi and H.E.S.S. data are inconsistent with a power-law γ-ray spectrum with an exponential cutoff. In particular, the H.E.S.S. spectrum five days after the optical peak follows a steep power-law (∝ E −q where q 2) or the very beginning (i.e., the low-energy, slowly falling portion) of an exponential cutoff. However, the overall normalization of this spectrum falls well below the overall normalization of the Fermi spectrum on the same day. As a result, theoretically-motivated models (i.e., power law spectra with exponential cutoffs) that fit the Fermi data will either overestimate the H.E.S.S. observations or, if  . γ-ray spectra from our single-shock model with parameters described in Table 1. The color scale of the lines and data points denotes the shock age and, in the case of the data points, the time after optical peak. A single shock, which produces a power-law particle distribution with an exponential cutoff, cannot simultaneously describe the Fermi and H.E.S.S. data.
a sufficiently small v sh is chosen to reduce the maximum energy, produce VHE spectra that is properly normalized but much steeper than that observed. To illustrate this issue, Figure 4 shows γ-ray spectra from our single shock model, which fit the Fermi data and have maximum energies in the needed range, but overestimate the VHE flux. We note that MAGIC (Acciari et al. 2022), HESS (Aharonian et al. 2022), andZheng et al. (2022) interpret the combined Fermi and VHE spectra as arising from a single shock. To fit these combined spectra, H.E.S.S. evokes a slow exponential cutoff (∝ e −(E/Emax) β ), where β < 1. This modification results in a good fit to the data but is not theoretically motivated. Meanwhile MAGIC does fit their spectra with a theoretically motivated power-law with exponential cutoff. However, they achieve this fit by invoking arbitrary normalizations and maximum energies that do not evolve in a physical manner.
A second key tension regarding the single-shock model is that it cannot reproduce the delay between the GeV and TeV peaks. As our model demonstrates, a single shock yields a luminosity peak that occurs at approximately the same time for all energies. However, as put forward in Aharonian et al. (2022), the VHE peak may be modulated by the maximum energy or, equivalently, by the finite acceleration time for TeV particles.
To illustrate why such a modulation cannot resolve the time delay issue, let us consider some simple scaling relations. Assuming E max is set by requiring that the acceleration time be approximately equal to the diffusion time, and that this time be less than or equal to the age of the system, we find that E max ∝ B 2 v 2 sh t for Bohm diffusion (Caprioli & Spitkovsky 2014c). Here, B 2 is the magnetic field behind the shock, ∝ ρ sh if the non-resonant streaming instability dominates magnetic field amplification. Thus, we find E max ∝ ρ 1/2 0 v 7/2 sh t. This framework for E max is broadly equivalent to that in our model, in which the diffusion length of the highestenergy particles is a fixed fraction of the system size.
Alternatively, since the non-resonant instability is driven by escaping particles, one can assume, as in Aharonian et al. (2022), that the number of escaping particles with energy E max is sufficient to drive the nonresonant instability and thereby inhibit particle propagation. This requirement yields E max ∝ ρ 1/2 0 v 2 sh R sh (see, e.g., Bell et al. 2011, for a detailed derivation). In either case, the maximum energy would be predicted to increase prior to the GeV luminosity peak, when the shock velocity and|given arguments presented earlier in this section|the density, are roughly constant. However, after the GeV luminosity peak, the shape of the light curve demands that the shock enter the Sedov stage. During this stage, both of the assumptions outlined above give a maximum energy that decreases with time, regardless of whether the ambient density is constant or follows a wind profile. Thus, a rise in the maximum particle energy cannot account for the delayed VHE luminosity peak.
It is also worth noting that, to produce good agreement with the GeV data, we require the onset of the Sedov-Taylor phase occur at t 1 day (see the gray lines in Figure 3). This result is inconsistent with the results of Pandey et al. (2022) and Cheung et al. (2022), which find that the free expansion phase lasts until approximately day 4. The fact that the TeV emission also does not peak until this time provides further evidence for multiple shock components.
Finally, we note that IC emission cannot resolve the issues mentioned here or, more to the point, contribute significantly to the γ-ray spectrum. Cosmic ray (CR) electrons suffer strong synchrotron losses in the amplified fields calculated in our model (∼ 1 G), severely reducing their ability to produce substantial VHE emission at any point during the shock's evolution. More importantly, the large radiation fields expected near the forward shock (∼ 1 erg cm −3 , see the supplementary materials of Aharonian et al. 2022) give an IC loss time 15 seconds for TeV electrons, less than their acceleration time.

A MULTI-SHOCK SCENARIO
The previous section demonstrates that a single external shock cannot describe both the GeV and TeV emission from RS Ophiuchi's 2021 outburst, unless spectra  Table 2. and maximum energy scaling are chosen ad-hoc rather than calculated based on the DSA theory. In this section, we instead explore a scenario involving two shocks which initially expand into different, roughly homogeneous media, but eventually probe the same RG wind on large scales. Specifically, we aim to test whether multiple shocks, be they internal interactions between distinct ejecta components or a manifestation of a single ejecta running into an aspherical medium, could feasibly explain the observed γ-ray emission. Conversely, this exercise shows how the observed γ-ray emission places a constraint on the environment within and surrounding the nova system. We shall focus our efforts only the simplest scenario-in terms of number of shock components-that reproduces the main features of the γ-ray observations. There may exist more complex scenarios that also yield a good fit to the GeV and TeV data. Our goal here is simply to provide further evidence for a complex outburst involving multiple shock components.  Right: Modeled spectra after one and five days. As in Figures 3 and 4, all times are given in terms of shock age for the model and relative to t0, the time of optical peak, for the observational data. This two-shock model reproduces both the γ-ray light curves and spectra observed by Fermi and H.E.S.S.
The simplest scenario that fits both the outburst's GeV and TeV light curves consists of two shocks: a slow, highly luminous component (i.e., a component that probes a relatively high density), and a fast, less luminous component (i.e., a component that probes a relatively low density).
The evolution of these components are shown in Figure 5 and the corresponding model parameters are listed in Table 2. We adopt the same RG wind parameters as in Section 3, which are consistent with the results of Iijima (2009). In this picture, the slow component produces the bulk of the GeV emission along with the very steep TeV spectrum at early times, while the fast component produces the hardened TeV emission at later times. This TeV hardening occurs because the fast component achieves both a higher maximum γ-ray energy and, since it sweeps up less mass at early times, a later luminosity peak. Note that the fast component of our multi-shock model has a Sedov-Taylor time of t 3 days, in better agreement with the deceleration timescale implied by X-ray observations (Pandey et al. 2022;Cheung et al. 2022).
The shock velocities we adopt are chosen to produce the best agreement with the observed γ-ray data. While these velocities are broadly consistent with optical data (specifically, v 1 and v 2 denoted in Section 2.5, they differ somewhat from those inferred from X-ray data (Cheung et al. 2022). This discrepancy may be the result of additional shock components, but may also be attributed to uncertainties arising from the conversion of X-ray tem-peratures to shock velocities (see Orio et al. 2022, for a detailed discussion).

Parameter Slow Component Fast Component
Mej 1 × 10 −7 M 1 × 10 −7 M v sh,init 1300 km s −1 4500 km s −1 n0,init 1.2 × 10 10 cm −3 5.0 × 10 7 cm −3 rcrit 1.0 AU 6.0 AU M wind 5 × 10 −7 M yr −1 5 × 10 −7 M yr −1 v wind 30 km s −1 30 km s −1 Table 2. Model parameters used to fit the Fermi and H.E.S.S. data, assuming the observed γ-ray emission arises from two shocks: a slow component that probes a relatively high ambient density, and a fast component that probes a relatively low ambient density. After expanding through these different regions, both components probe the same RG wind. The resulting hydrodynamic evolution can be found in Figure 5. Note that, as with the single shock model, a smoothing function is applied to the density profile, such that the density-related parameters shown here might differ slightly from the ultimate profile that enters into our model.
The resulting light curves and spectra are shown in Figure 6. This two-shock model yields good agreement with both the GeV and TeV observations. In particular, because the fast component takes longer to reach its luminosity peak and has a higher E max at that time, it reproduces the observed TeV delay.
For the sake of simplicity, we take the total γ-ray flux to be the sum of the fluxes from the two components. A different set of model parameters might also fit the data if the relative weights of the two components are modified (i.e., if each component has a angular filling fraction f ≡ ∆Ω/4π not equal to unity, where ∆Ω is the solid angle subtended by the shock), or if the two shocks are launched at different times. Broadly speaking, however, the two components must satisfy several conditions: 1. Shock velocities must differ by a factor of 3. This requirement comes from the need for sufficient stratification in E max (i.e., the fast component must produce substantial TeV emission while the GeV component must not). If the two shock velocities are too similar, we revert to a scenario akin to a single-shock.
2. Assuming the two shocks are launched simultaneously, the extent of the inner homogeneous region probed by the fast component, r crit , is 6 times larger than that of the slow component. As discussed in Section 3, the peak of a component's light curve roughly corresponds to when the shock starts to slow down and to probe the RG wind. Taking the shock age at the TeV peak to be roughly twice that at the GeV peak (consistent with observations), and the fast shock to have a velocity that is 3 times that of the slow shock, we obtain r crit,f /r crit,s 6, where subscripts f and s denote the fast and slow components, respectively.
3. Assuming equal contributions from the two components (i.e., each has the same filling factor, f ), the uniform densities probed by the two components at small radii must differ by a factor of at least a few hundred. The GeV peak luminosity, produced by the slow component, is approximately a thousand times the TeV peak luminosity, produced primarily by the fast component. With the modestly steep photon spectrum dN/dE ∝ E −2.3 returned by our self-consistent model (see Caprioli et al. 2020;Diesing & Caprioli 2021, for the physics behind the obtained spectral slope), we expect the GeV luminosity of the fast component to be roughly 10 times its TeV luminosity. Thus, at GeV energies, L γ,f /L γ,s ∼ 10 −2 , where these L γ are the maximum luminosities of each component. Since, at peak, L γ ∝ ρ 2 0 v 2 sh r 3 crit , our constraints on the ratios of v sh and r crit give ρ 0,f /ρ 0,s 2 × 10 −3 . If we additionally introduce filling factors, f f and f s such that f f = 1 − f s and L γ ∝ f for each component (equivalent to the presence of polar and equatorial components of the shock), we can recast this condition as ρ 0,f /ρ 0,s 2 × 10 −3 (1/f f − 1) 1/2 . Note that the model parameters described in Table 2 differ slightly from this estimate, primarily due to the fact that our model is multi-zone and therefore includes contributions to L γ from multiple epochs of shock evolution.
4. In the case that f f = f s = 0.5, the two components must have comparable ejecta mass (i.e., within a factor of a few). Recall that, to reproduce the shape of the peak, both components must enter the Sedov stage when they reach the extent of their homogeneous region, r crit . Thus, the ejecta mass, M ej must be equal to the mass contained in the homogeneous region, ∝ f ρ 0 r crit . For equal filling factors, the constraints enumerated above yield M ej,f 0.4M ej,s . Otherwise, we have The above arguments make clear that a degeneracy exists between the relative ejecta masses of the fast and slow components and their relative angular filling factors. A similar degeneracy also exists between these filling factors and the relative ambient densities. Nevertheless, the main conclusion of this exercise still stands: while a single shock cannot reproduce the γ-ray observations, scenarios with multiple shocks (or multiple shock components) are capable of doing so.

CONCLUSION
We have modeled the γ-ray emission of RS Ophiuchi's recent outburst using a semi-analytic model of particle acceleration that self-consistently accounts for magnetic field amplification as well as the back-reaction of nonthermal particles. We demonstrate the properties of the observed γ-ray emission is not consistent with a single, external shock. This inconsistency arises from the facts that: a) the finite acceleration time of TeV particles is too short to explain the delay between the GeV and TeV peaks as observed by Fermi-LAT and H.E.S.S.; and (b) the combined GeV-TeV spectra are inconsistent with theoretically-motivated emission models for a single shock, in particular a power-law with an exponential cutoff.
On the other hand, we find that the observed γ-ray emission is naturally reproduced in a scenarios involving multiple shocks. In particular, both the spectra and GeV/TeV light curves are consistent with the combined emission from two shocks: one with a low initial velocity expanding into a dense ambient medium, and one with a fast initial velocity expanding into a comparatively rarefied medium. Different combinations of shocks with non-equal filling factors may also be able to reproduce the observations, as could scenarios with three or more shocks. The key takeaway, then, is that RS-Ophiuchi's recent outburst must be more complicated than the single-shock scenario presented in the literature.
The presence of multiple ejecta components in RS Ophiuchi (and, hence, multiple shocks) is also supported by time-dependent optical spectroscopy, which reveals as many as three distinct velocity components (Sec. 2.5). Multiple shocks are also supported by X-ray observations (Page et al. 2022;Orio et al. 2022); Orio et al. (2022) find evidence for at least two components of shocked plasma, at temperatures kT ≈ 0.75 keV and ≈ 3 keV, respectively, which could in principle be attributed to the "slow" and "fast" shocks in our twoshock scenario. However, the fact that the peak of the hard X-ray light curve is delayed relative to that of either the GeV or TeV γ-ray emission, is challenging to understand and may point to the presence of yet additional shocks. Thermal X-rays from the shocks which dominate the γ-ray emission be absent at early times if they are either absorbed by a large neutral gas column ahead of the shock, or intrinsically suppressed in the multi-phase post-shock region (e.g., Steinberg & Metzger 2018;Nelson et al. 2019).
The ejecta from classical novae are also characterized by multiple velocity components (e.g., Aydi et al. 2020a), with the slower ejecta concentrated in the equatorial plane of the binary and the faster component likely representing a more spherical wind from the white dwarf which then expands more freely into the polar directions (e.g., Chomiuk et al. 2014). If embedded novae like RS Ophiuchi generate ejecta with a similar geometry, then multiple "slow" and "fast" shocks would be generated as these outflows collide with themselves (as in classical novae) and with the surrounding circumbinary environment, which on radial scales interior to the binary separation is also highly aspherical, being significantly denser in the binary equatorial plane than along the polar directions (e.g., Booth et al. 2016). If this interpretation is correct, then one lesson to be gleaned from RS Ophiuchi is that−despite vastly different circumbinary environments−the ejecta properties of embedded/symbiotic novae may be similar to classical novae. This in turn might post a challenge to models in which the binary companion (whose orbital separation from the white dwarf is significantly smaller in classical novae than in symbiotic/embedded novae like RS Ophiuchi) plays a dominant role in shaping the equatorial ejecta in classical novae (e.g., MacDonald 1980;Livio et al. 1990).
Given that symbiotic novae are often treated as rapidly-evolving analogs for supernova remnants (e.g., Martin & Dubus 2013), the complex behavior of RS Ophiuchi revealed in this work is also an important reminder that novae are fundamentally different systems with their own unique properties. That being said, some of this behavior may also be relevant for young supernovae, particularly those expanding into nonuniform media (e.g., Thomas et al. 2022). ACKNOWLEDGMENTS R.D., D.C., and S.G. were partially supported by NASA (grants NNX17AG30G, 80NSSC18K1218, and 80NSSC18K1726) and the NSF (grants AST-1714658, AST-1909778, PHY-1748958, and PHY-2010240