tiny [table]font=tiny

# Time-resolved hadronic particle acceleration in the recurrent nova RS
Ophiuchi

H.E.S.S<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
<EMAIL_ADDRESS>222H.E.S.S. Collaboration authors and
affiliations are listed in the supplementary materials

> Recurrent novae are repeating thermonuclear explosions in the outer layers
> of white dwarfs, due to the accretion of fresh material from a binary
> companion. The shock generated by ejected material slamming into the
> companion star’s wind, accelerates particles to very-high-energies. We
> report very-high-energy (VHE, $\gtrsim 100$ GeV) gamma rays from the
> recurrent nova RS Ophiuchi up to a month after its 2021 outburst, using the
> High Energy Stereoscopic System. The VHE emission has a similar temporal
> profile to lower-energy GeV emission, indicating a common origin, with a
> two-day delay in peak flux. These observations constrain models of time-
> dependent particle energization, favouring a hadronic emission scenario over
> the leptonic alternative. This confirms that shocks in dense winds provide
> favourable environments for efficient cosmic-ray acceleration to very-high-
> energies.

RS Ophiuchi (RS Oph) is a recurrent nova system comprising a white dwarf and a
companion red giant star. Novae are a source of high-energy particles (?, ?),
with non-thermal gamma-ray emission in the range $\sim$100 MeV to $\sim$10 GeV
(?). The RS Oph system is located approximately $1.4$kpc from Earth (?);
analysis of Gaia data suggests larger distances of $2.3$ kpc (?) or $2.7$ kpc
(?), although the reliability of these estimates is questionable due to the
orbital motion in RS Oph. The binary components have a separation of $1.48$
astronomical units (au), close enough for the white dwarf to continually
accrete material from its companion (?). At irregular intervals, enough
material accumulates on the surface of the white dwarf to trigger a
thermonuclear explosion, driving a quasi-spherical shock into the red giant’s
wind (see ?, Figure 2). Eight outbursts were observed between 1898 and 2006,
recurring in intervals of $9\,$to$\,26$ years (?).

On $8^{\rm th}$ August 2021, an outburst of RS Oph was identified from optical
observations (?), with a peak naked-eye visual magnitude of 4.5, over a
thousand times brighter than the quiescent visual magnitude of 12.5. We
observed RS Oph with the High Energy Stereoscopic System (H.E.S.S.), an array
of five atmospheric Cherenkov telescopes. Observations commenced on $9^{\rm
th}$ August 2021 and continued for five nights until $13^{\rm th}$ August
2021. Optical background emission from the Moon prevented good quality
observations for the following ten nights. During each of those five nights,
H.E.S.S. detected point-like gamma-ray emission from the direction of RS Oph,
with a significance of $>6$ sigma on each night (see ?, Table S1). The data
for those five nights combined are shown in Figure 1. Observations recommenced
on $25^{\rm th}$ August 2021, $\sim$17 days after the initial outburst. We
find evidence for a much weaker signal ($\sim 3\sigma$ above the background)
was seen in $\sim$15 hours (after quality cuts) of data accumulated over the
following 14 days.

Figure 1: RS Oph significance maps. Significance maps derived from the
H.E.S.S. $>100$ GeV gamma-ray observations for the early (A) and late (B)
phases of the RS Oph 2021 outburst. T${}_{0}=$ Modified Julian Day (MJD)
59435.25, is the time of peak optical emission. The dashed white circles
indicate the point-spread-function (PSF).

We performed a spectral analysis of the H.E.S.S. data for the first five
observation nights separately. We also separated the data from the array of
four $106\,\textrm{m}^{2}$ mirror area telescopes (designated CT1-4) and the
fifth $612\,\textrm{m}^{2}$ mirror area low-threshold telescope (CT5). We find
that the VHE flux is variable, with a spectral index $>3$ throughout (see ?,
Table S2).

Figure 2: Gamma-ray light curves of RS Oph. Light curves of gamma-ray emission
from RS Oph including data from _Fermi_ -LAT and H.E.S.S. observations. The
H.E.S.S. data (red squares) cover a period of five nights, after which
observations ceased for ten days due to bright moonlight, marked by the shaded
grey band, then recommenced for a period of 14 days. The H.E.S.S. flux is
integrated from 250 GeV to 2.5 TeV, whilst the _Fermi_ -LAT flux is integrated
from 60 MeV to 500 GeV. _Fermi_ -LAT data are shown in 6-hour bins (blue
circles) corresponding to the time windows of the H.E.S.S. observations, and
data outside of these times shown with semi-transparent markers. Error bars
are 1 $\sigma$ statistical uncertainties. A power-law slope model was fitted
to the temporal decay after the time of peak flux for both instruments (red
and blue dashed lines, with uncertainties indicated by the shaded regions).
The vertical dotted black line indicates the peak of the outburst in the
optical waveband, T0.

Figure 2 shows the time evolution of the gamma-ray flux curve for photon
energies between 250 GeV and 2.5 TeV. The VHE gamma-ray flux rises smoothly
from T0, the time of peak optical emission in the V band (?), until a VHE peak
on the third night of observations, after which the VHE gamma-ray energy flux
decays by an order of magnitude over a two-week period. We obtained 60 MeV –
500 GeV data taken by the _Fermi_ -LAT (Large Area Telescope) instrument for
the same time period as the H.E.S.S. observations which are also shown in
Figure 2. The flux varies in the range $\sim$1$\times$10-8 – 2$\times$10-10
erg cm-2s-1, with a peak flux in the _Fermi_ -LAT data on $T_{0}+1$ day. The
VHE gamma-ray emission peak is delayed by a further two days.

After the peak flux, we fitted the decay in time $t$ of the energy flux with a
power-law with exponent $\alpha$, $t^{-\alpha}$ and found best-fitting values
of $\alpha\approx 1.3-1.4$ in both data sets: $\alpha_{\rm HESS}=1.43\pm 0.18$
for H.E.S.S. and $\alpha_{\rm LAT}=1.31\pm 0.07$ for _Fermi_ -LAT, for the
choice of $T_{0}=$1 day. The _Fermi_ -LAT flux and temporal decay are
consistent with that obtained from bins of 24-hour duration, the higher
statistics of the larger bins enabling the detailed spectral analysis shown in
Figure 3.

The combined H.E.S.S. and _Fermi_ -LAT data (?) allow us to measure wide-band
gamma-ray spectra over more than four orders of magnitude in energy and follow
their temporal evolution (Figure 3). The RS Oph spectra are consistent with a
log-parabola model. Comparison of spectra taken on different nights show a
general trend for the flux normalisation to decrease and the parabola to widen
over time (see ?, Table S2).

Figure 3: RS Oph gamma-ray spectra. The H.E.S.S. and _Fermi_ -LAT spectra for
9 (green) and 13 (orange) August fitted with a log-parabola model. The
analysis is applied separately for the H.E.S.S. CT1-4 (squares) and CT5
(triangles) - see text. The _Fermi_ -LAT data (open circles) are integrated
over $24$ h centred at the H.E.S.S. observation times. There is clear spectral
evolution from the 9th to the 13th August, with a noticeable reduction in the
_Fermi_ -LAT flux as well as an increase in the maximum energy of the TeV
spectrum. Error bars are 1 sigma statistical uncertainty, and upper limits are
the $95\%$ confidence level.

The similarity between the spectra of the _Fermi_ -LAT and H.E.S.S. data, and
their similar decay profiles after their respective peaks, indicate a common
origin for the gamma rays from one day to one month after the explosion. We
assume that the particles that generate the gamma rays are accelerated at the
external shock as it propagates into the wind of the red giant (?, Figure S2).
Optical spectroscopic measurements of the 2021 nova indicate shock velocities
in the range $u_{\rm sh}=4000-5000$ km s-1 (?), compatible with measurements
from the previous 2006 outburst of RS Oph (?, ?). High resolution images of
the 2006 event (?) indicated the polar regions of the shock expanded at
$\approx 5\,000$ km s-1 over the first 5 months. We therefore assume that
during the first week following the 2021 outburst the shock velocity did not
fall below several thousand kilometers per second.

The images of the 2006 nova showed a quasi-spherical outflow, pinched at an
equatorial ring (?, ?). This is consistent with a shock expanding into the
wind of the red giant orthogonal to the orbital plane of the binary, but
inhibited close to the plane by the denser gas (?, ?). We expect particles to
undergo diffusive shock acceleration at the external fast moving shocks, above
and below the orbital plane of the binary. We consider two scenarios to
explain the observed spectral and temporal properties of the gamma-ray
emission from RS Oph: gamma-ray production from accelerated protons colliding
with dense gas in the downstream volume (hadronic $\pi^{0}$ decay model), or
gamma-ray emissions from energetic electrons scattering low energy photons
from the nova (inverse Compton model). For both models, the observations place
strong constraints on the physical conditions, particularly the acceleration
efficiencies required to match the measured fluxes and maximum photon energies
(?).

VHE gamma-ray emission requires acceleration of particles to $>$TeV energies.
The maximum energy a particle attains at a shock is determined by either the
time taken before radiative cooling dominates over acceleration, or by when
particles become too energetic and escape upstream of the shock (?). This
confinement limit applies when the accelerating particles are unable to excite
magnetic field fluctuations to a sufficient level ahead of the shock. As
particles spend longer diffusing upstream than downstream, the details of the
downstream magnetic fields can be neglected (?). For upstream magnetic-field
amplification to be effective, a sufficient flux of particles, typically
protons, must escape upstream of the shock. This requires both efficient
transfer of the shock kinetic energy to relativistic protons, and that a
fraction of these protons penetrate far upstream. Escaping particles have
energies concentrated close to the maximum particle energy; less energetic
particles are confined to the shock. The escaping flux per unit area at a
given shock radius can be parameterized as

${j}_{\rm esc}=e\xi_{\rm esc}F_{\varepsilon}/E_{\rm max}~{},$ (1)

where $e$ is the elementary charge, $F_{\varepsilon}=\frac{1}{2}\rho_{\rm
up}u_{\rm sh}^{3}$ is the energy flux density for a high Mach-number non-
relativistic shock, $\rho_{\rm up}$ is the immediate upstream gas density at
that radius, $E_{\rm max}$ is the maximum particle energy which dominates the
escaping flux, and the efficiency parameter $\xi_{\rm esc}$ depends on the
assumed particle spectrum (?). For a wind-like density profile, and neglecting
radiative losses, the confinement limit on the maximum energy for a particle
of charge $|Z|e$ (with atomic number $Z$) is

$E_{\rm max}=1.5|Z|\left(\frac{\xi_{\rm
esc}}{0.01}\right)\left(\frac{\dot{M}/v_{\rm wind}}{10^{11}~{}{\rm
kg~{}m}^{-1}}\right)^{1/2}\left(\frac{u_{\rm sh}}{5000~{}{\rm
km~{}s}^{-1}}\right)^{2}~{}{\rm TeV}\enspace,$ (2)

where $\dot{M}$ and $v_{\rm wind}$ are the mass-loss rate and the wind
velocity of the red giant. $\xi_{\rm esc}$ is predicted to be about $1\%$ for
high Mach-number shocks (?). For RS Oph, $\dot{M}/v_{\rm wind}=6\times
10^{11}~{}{\rm kg~{}m}^{-1}$ (?) which, together with the inferred shock
velocities, indicates a maximum energy $E_{\rm max}\approx 10$ TeV. This is
compatible with the measured maximum photon energies $E_{\rm\gamma,max}\approx
1$ TeV (Figure 3).

In the hadronic scenario, the gamma-ray lightcurves are consistent with an
expanding shock in a decreasing density profile. With the adopted distance of
1.4 kpc (?), the measured gamma-ray fluxes require that $>$10% of the post-
shocked medium’s internal energy goes to accelerating protons or other nuclei.
The delay between the peaks in the _Fermi_ -LAT and H.E.S.S. lightcurves would
then reflect the finite acceleration time of the $>$1 TeV protons, or more
specifically, the time taken to populate the high-energy tail of the
distribution (?). A simple calculation of the acceleration time, based on a
comparison of the confinement and Hillas limits, implies that this should
happen on the order of days (?). This is consistent with the spectral
evolution seen in Figure 3: a reduction in the _Fermi_ -LAT flux, accompanied
by a hardening in H.E.S.S. flux and increased $E_{\rm\gamma,max}$ over the
first few days post-outburst. Attenuation of gamma-rays due to the novas
optical and infra-red photon fields is minor below 1 TeV a few hours after the
explosion, and therefore attenuation alone cannot account for the observed
hardening (see ?, Figure S10).

For the alternative, leptonic scenario in which TeV gamma-rays are produced by
VHE electrons, the acceleration needs to overcome the strong radiative losses
due to inverse Compton cooling in the strong photon fields of the nova, as
well as synchrotron cooling in the magnetic field in the shock region. To
achieve this, electrons must accelerate at close to the Bohm rate, i.e. the
scattering rate equal to the rate of gyration in the magnetic field (?). Such
efficient scattering requires strong self-generated magnetic fluctuations
upstream of the shock, which implies the presence of an energetic relativistic
hadronic component. In this scenario, the differences between the spectral
slopes in the _Fermi_ -LAT and H.E.S.S. energy ranges are a consequence of the
energy-dependent cooling rates in time-dependent photon fields. Electrons that
radiate in the VHE band cool on a timescale less than the age of the nova
remnant at the times the observations were taken, while lower-energy un-cooled
electrons accumulate downstream over time. The _Fermi_ -LAT light curve in
this scenario then reflects the evolution of the energy density of soft-photon
targets, while the H.E.S.S. lightcurve traces the full radiative output of
high-energy electrons up to the VHE peak. After the peak, due to the rapidly
decreasing photon energy density, the cooling time increases faster than the
remnant’s age, and the VHE emitting electrons are also slow cooling.

From the time-dependent numerical single-zone model, parameters can be found
to approximately describe the light curves and spectra in both leptonic and
hadronic scenarios, providing quantitative estimates for the acceleration
efficiencies at early times $t$ ($t\leq T_{0}+5$ days). We are unable to
account for the temporal decay at later times since a more complex treatment
of the internal structure of the nova remnant, its non-spherical geometry, and
escape of particles from the emission zone is required. Both the leptonic and
hadronic models at early times are consistent with continuous injection of
particles following a power law spectrum in energy $E^{-2.2}$ with a high-
energy cut-off. To approximately match the observed flux, the hadronic model
requires $>$10% of the shocked gas’ internal energy to be transferred to non-
thermal protons, while the leptonic model requires $>$1% efficiency for non-
thermal electrons.

Such a high fraction of the total energy in non-thermal electrons is
inconsistent with theories of injection at high Mach-number shocks, for which
ion injection efficiency is expected to be much higher than that of electrons
(?). Numerical simulations of high Mach-number shocks find a ratio of
$\ll$10-2 for electron to ion energy densities (?), which is consistent with
multi-wavelength models of supernova remnants (?). A $>$1% efficiency of
conversion to non-thermal electrons cannot be realised in a purely leptonic
model. For this reason, we prefer the hadronic scenario discussed above, for
which both the implied high proton acceleration efficiencies and inferred
maximum energy are in line with theoretical predictions (?). Our findings
support previous hadronic models of gamma-ray novae (?, ?, ?).

The VHE detection of RS Oph demonstrates that particle acceleration to
$\sim$TeV energies can occur within the dense wind environments of recurrent
novae. The total kinetic energy from each nova of RS Oph is estimated to be
$\sim 10^{43}$ erg ($10^{-7}$ solar masses of ejecta at $\sim$4000 km s-1
(?)), with a large fraction of this being converted to relativistic protons
and heavier nuclei, which are the main constituents of Galactic cosmic rays.
Each nova event generates enough cosmic rays to fill a cubic parsec (pc)
volume with an energy density of $\sim$0.1 eV cm-3, similar to the local
Galactic cosmic-ray energy density of $\sim(0.8-1.0)$ eV cm-3 (?) sustained by
supernovae. In the case of RS Oph, the cosmic-ray energy input recurs
approximately every $\Delta t=15-20$ years, leading to an almost continuous
injection of non-thermal particles. For a diffusion coefficient $D$ in the
neighbourhood of RS Oph, the cosmic-ray output from each nova is spread over a
diffusion length $\ell_{\rm diff}=\sqrt{4D~{}\Delta t}$. Using a Galactic
average $D\approx(3-5)\times 10^{28}$ cm2 s-1 (?), we find $\ell_{\rm diff}>1$
pc, and the contribution from each nova is sub-dominant to the average
Galactic cosmic-ray population. If the diffusion coefficient in the
neighbourhood of the nova is much lower than the Galactic average, due to
enhanced turbulence following previous novae for example, such a sustained
source of cosmic rays will raise the local abundance. If efficient
acceleration of particles to TeV energies in recurrent novae is commonplace,
with spectral energy distribution harder than that of the Galactic cosmic-ray
background $\propto E^{-2.7}$, the local contribution from novae would
dominate at TeV energies over volumes $>{\rm pc}^{3}$. The size of the
affected region will depend on the value of the diffusion coefficient, which
can be constrained using measurements of the diffuse gamma-ray emission at
energies $\sim 10-100$ GeV (?).

Our time-resolved gamma-ray emission measurements have implications for the
origin of cosmic rays. Acceleration of cosmic rays to PeV energies in
supernova remnants requires substantial amplification of magnetic fields. Fast
shocks ($\sim 10,000$ km s-1) propagating through the dense winds
($\dot{M}/v_{\rm wind}\sim 10^{13}~{}{\rm kg~{}m}^{-1}$) associated with the
progenitors of supernova remnants from massive ($\gtrsim
8\,\mathrm{M}_{\odot}$) stars, provide the only known environments where the
required conditions can (in theory) be met (?, ?). However, observational
confirmation of this prediction has not been found. The detection of VHE
gamma-rays from RS Oph provides an example of a Galactic accelerator reaching
the theoretical limit for the maximum achievable particle energy via diffusive
shock acceleration (?). If our results can be extrapolated to the most
optimistic supernova conditions, they support the prevailing model of Galactic
PeV cosmic rays originating in supernova remnants from massive stars (?, ?).

## References

  * 1. A. A. Abdo, et al., Science 329, 817 (2010).
  * 2. M. Ackermann, et al., Science 345, 554 (2014).
  * 3. L. Chomiuk, B. D. Metzger, K. J. Shen, Annual Review of Astronomy and Astrophysics 59, 391 (2021).
  * 4. R. K. Barry, et al., RS Ophiuchi (2006) and the Recurrent Nova Phenomenon, A. Evans, M. F. Bode, T. J. O’Brien, M. J. Darnley, eds. (2008), vol. 401 of Astronomical Society of the Pacific Conference Series, p. 52.
  * 5. Gaia Collaboration, et al., A&A 616, A1 (2018).
  * 6. L. Lindegren, et al., A&A 649, A4 (2021).
  * 7. R. A. Booth, S. Mohamed, P. Podsiadlowski, MNRAS 457, 822 (2016).
  * 8. Materials and methods are available as supplementary materials.
  * 9. E. Brandi, C. Quiroga, J. Mikołajewska, O. E. Ferrer, L. G. García, A&A 497, 815 (2009).
  * 10. S. Kafka, Observations from the AAVSO International Database, https://www.aavso.org (2021).
  * 11. J. Mikolajewska, E. Aydi, D. Buckley, C. Galan, M. Orio, The Astronomer’s Telegram 14852, 1 (2021).
  * 12. M. F. Bode, et al., ApJ 652, 629 (2006).
  * 13. J. L. Sokoloski, G. J. M. Luna, K. Mukai, S. J. Kenyon, Nature 442, 276 (2006).
  * 14. M. F. Bode, et al., ApJ 665, L63 (2007).
  * 15. T. J. O’Brien, et al., Nature 442, 279 (2006).
  * 16. R. Walder, D. Folini, S. N. Shore, A&A 484, L9 (2008).
  * 17. A. R. Bell, K. M. Schure, B. Reville, G. Giacinti, MNRAS 431, 415 (2013).
  * 18. P. O. Lagage, C. J. Cesarsky, A&A 125, 249 (1983).
  * 19. L. O. Drury, MNRAS 251, 340 (1991).
  * 20. M. A. Malkov, L. O. Drury, Reports on Progress in Physics 64, 429 (2001).
  * 21. J. Park, D. Caprioli, A. Spitkovsky, Phys. Rev. Lett. 114, 085003 (2015).
  * 22. H. J. Völk, E. G. Berezhko, L. T. Ksenofontov, A&A 483, 529 (2008).
  * 23. B. D. Metzger, et al., MNRAS 457, 1786 (2016).
  * 24. K.-L. Li, et al., Nature Astronomy 1, 697 (2017).
  * 25. E. Aydi, et al., Nature Astronomy 4, 776 (2020).
  * 26. A. C. Cummings, et al., ApJ 831, 18 (2016).
  * 27. A. W. Strong, I. V. Moskalenko, V. S. Ptuskin, Annual Review of Nuclear and Particle Science 57, 285 (2007).
  * 28. M. Ackermann, et al., ApJ 750, 3 (2012).
  * 29. A. Marcowith, V. V. Dwarkadas, M. Renaud, V. Tatischeff, G. Giacinti, MNRAS 479, 4470 (2018).
  * 30. F. Aharonian, et al., Astronomy and Astrophysics 457, 899 (2006).
  * 31. M. Holler, et al., PoS ICRC2015, 847 (2016).
  * 32. T. Ashton, et al., Astroparticle Physics 118, 102425 (2020).
  * 33. B. Bi, et al., Proceedings of 37th International Cosmic Ray Conference — PoS(ICRC2021) (2021).
  * 34. C. C. Cheung, S. Ciprini, T. J. Johnson, The Astronomer’s Telegram 14834, 1 (2021).
  * 35. K. Taguchi, T. Ueta, K. Isogai, The Astronomer’s Telegram 14838, 1 (2021).
  * 36. U. Munari, P. Valisa, The Astronomer’s Telegram 14840, 1 (2021).
  * 37. J. Hahn, et al., Astroparticle Physics 54, 25 (2014).
  * 38. D. Berge, S. Funk, J. Hinton, Astronomy and Astrophysics 466, 1219 (2007).
  * 39. S. Ohm, C. van Eldik, K. Egberts, Astroparticle Physics 31, 383 (2009).
  * 40. T. Murach, M. Gajdus, R. D. Parsons, Proceedings of the 34th International Cosmic Ray Conference (ICRC2015), The Hague, The Netherlands (2015).
  * 41. G. Pühlhofer, et al., Proceedings of 37th International Cosmic Ray Conference — PoS(ICRC2021) (2021).
  * 42. R. D. Parsons, J. A. Hinton, Astroparticle Physics 56, 26 (2014).
  * 43. L. Mohrmann, et al., A&A 632, A72 (2019).
  * 44. C. Nigro, et al., A&A 625, A10 (2019).
  * 45. C. Deil, et al., 35th International Cosmic Ray Conference (ICRC2017) (2017), vol. 301 of International Cosmic Ray Conference, p. 766.
  * 46. M. de Naurois, L. Rolland, Astroparticle Physics 32, 231 (2009).
  * 47. F. Aharonian, et al., A&A 449, 223 (2006).
  * 48. W. B. Atwood, et al., ApJ 774, 76 (2013).
  * 49. M. Wood, et al., 35th International Cosmic Ray Conference (ICRC2017) (2017), vol. 301 of International Cosmic Ray Conference, p. 824.
  * 50. Fermi Science Support Development Team, Fermitools: Fermi Science Tools (2019).
  * 51. S. Abdollahi, et al., ApJS 247, 33 (2020).
  * 52. http://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html.
  * 53. E. Aydi, et al., ApJ 905, 62 (2020).
  * 54. R. Das, D. P. K. Banerjee, N. M. Ashok, ApJ 653, L141 (2006).
  * 55. M. F. Bode, F. D. Kahn, MNRAS 217, 205 (1985).
  * 56. L. I. Sedov, Similarity and Dimensional Methods in Mechanics (Academic Press, New York, 1959).
  * 57. X. Tang, R. A. Chevalier, MNRAS 465, 3793 (2017).
  * 58. L. D. Landau, E. M. Lifshitz, Fluid mechanics (Pergamon Press Ltd., 1959).
  * 59. T. J. O’Brien, M. F. Bode, F. D. Kahn, MNRAS 255, 683 (1992).
  * 60. D. Proga, S. J. Kenyon, J. C. Raymond, ApJ 501, 339 (1998).
  * 61. H. J. G. L. M. Lamers, J. P. Cassinelli, Introduction to Stellar Winds (Cambridge University Press, 1999).
  * 62. U. Munari, P. Valisa, arXiv e-prints p. arXiv:2109.01101 (2021).
  * 63. L. O. Drury, Reports on Progress in Physics 46, 973 (1983).
  * 64. B. Reville, S. O’Sullivan, P. Duffy, J. G. Kirk, MNRAS 386, 509 (2008).
  * 65. E. N. Parker, ApJ 128, 664 (1958).
  * 66. A. J. Kemball, P. J. Diamond, ApJ 481, L111 (1997).
  * 67. A. M. Hillas, ARA&A 22, 425 (1984).
  * 68. A. R. Bell, MNRAS 353, 550 (2004).
  * 69. G. R. Blumenthal, R. J. Gould, Reviews of Modern Physics 42, 237 (1970).
  * 70. A. Evans, et al., ApJ 671, L157 (2007).
  * 71. E. Kafexhiu, F. Aharonian, A. M. Taylor, G. S. Vila, Phys. Rev. D 90, 123014 (2014).

## Acknowledgements

Acknowledgements We appreciate the excellent work of the technical support
staff in Berlin, Zeuthen, Heidelberg, Palaiseau, Paris, Saclay, Tübingen and
in Namibia in the construction and operation of the equipment. This work
benefited from services provided by the H.E.S.S. Virtual Organisation,
supported by the national resource providers of the EGI Federation.

Funding The support of the Namibian authorities and of the University of
Namibia in facilitating the construction and operation of H.E.S.S. is
gratefully acknowledged, as is the support by the German Ministry for
Education and Research (BMBF), the Max Planck Society, the German Research
Foundation (DFG), the Helmholtz Association, the Alexander von Humboldt
Foundation, the French Ministry of Higher Education, Research and Innovation,
the Centre National de la Recherche Scientifique (CNRS/IN2P3 and CNRS/INSU),
the Commissariat à l’énergie atomique et aux énergies alternatives (CEA), the
U.K. Science and Technology Facilities Council (STFC), the Irish Research
Council (IRC) and the Science Foundation Ireland (SFI), the Knut and Alice
Wallenberg Foundation, the Polish Ministry of Education and Science, agreement
no. 2021/WK/06, the South African Department of Science and Technology and
National Research Foundation, the University of Namibia, the National
Commission on Research, Science & Technology of Namibia (NCRST), the Austrian
Federal Ministry of Education, Science and Research and the Austrian Science
Fund (FWF), the Australian Research Council (ARC), the Japan Society for the
Promotion of Science, the University of Amsterdam and the Science Committee of
Armenia grant 21AG-1C085.

Author Contributions: A. Mitchell and S. Ohm led the H.E.S.S. observations of
RS Oph. R. Konno carried out initial on-site analysis, the main H.E.S.S. CT1-4
stereo data analysis and the atmospheric correction. S. Steinmassl performed
the CT5 mono analysis, and J.P. Ernenwein the cross-check analysis. S. Ohm
coordinated the multiple H.E.S.S. analyses and evaluated the systematic
errors. E. de Oña Wilhelmi and T. Unbehaun carried out the _Fermi_ -LAT data
analysis. B. Reville, D. Khangulyan, J.Mackey and E. de Oña Wilhelmi developed
the interpretation and modelling. The manuscript was prepared by A. Mitchell,
B. Reville, S. Ohm, D. Khangulyan, E. de Oña Wilhelmi, R. Konno, and S.
Steinmassl. S. Wagner is the collaboration spokesperson. Other H.E.S.S.
collaboration authors contributed to the design, construction and operation of
H.E.S.S., the development and maintenance of data handling, data reduction or
data analysis software. All authors meet the journal’s authorship criteria and
have reviewed, discussed, and commented on the results and the manuscript.

Competing Interests: The authors declare that they have no competing
interests.

Data and Materials Availability: The H.E.S.S. data are available at:
https://www.mpi-hd.mpg.de/hfm/HESS/pages/publications/auxiliary/2022_RS_Oph/.
This includes the sky maps (cf. Figure 1), the light-curve (cf. Figure 2), the
points of the spectral energy distributions (cf. Figure 3).

Supplementary Material: Authors and affiliations, Materials and Methods,
Figure S1-S10, Tables S1-S4, References (30-71)

Supplementary Material for

Time-resolved hadronic particle acceleration in the recurrent nova RS Ophiuchi

H.E.S.S. Collaboration*

*Correspondence to<EMAIL_ADDRESS>Alison Mitchell (alison.mw.mitchell@fau.de), Stefan Ohm (stefan.ohm@desy.de), Brian Reville<EMAIL_ADDRESS>
This PDF file includes:

Materials and Methods
Supplementary Text
Figs. S1 to S10
Tables S1 to S4

### H.E.S.S. Collaboration authors and affiliations

F. Aharonian1,2,3, F. Ait Benkhali4, E.O. Angüner5, H. Ashkar6, M. Backes7,8,
V. Baghmanyan9, V. Barbosa Martins10, R. Batzofin11, Y. Becherini12,13, D.
Berge10, K. Bernlöhr2, B. Bi14, M. Böttcher8, C. Boisson15, J. Bolmont16, M.
de Bony de Lavergne17, M. Breuhaus2, R. Brose1, F. Brun18, S. Caroff16, S.
Casanova9, M. Cerruti 12, T. Chand8, A. Chen11, G. Cotter19, J. Damascene
Mbarubucyeye10, A. Djannati-Ataï12, A. Dmytriiev15, V. Doroshenko14, C.
Duffy19, K. Egberts20, J.-P. Ernenwein5, S. Fegan6, K. Feijen21, A. Fiasson17,
G. Fichet de Clairfontaine15, G. Fontaine6, M. Füßling10, S. Funk22, S.
Gabici12, Y.A. Gallant23, S. Ghafourizadeh4, G. Giavitto10, L. Giunti12,18, D.
Glawion22, J.F. Glicenstein18, M.-H. Grondin24, G. Hermann2, J.A. Hinton2, M.
Hörbe19, W. Hofmann2, C. Hoischen20, T. L. Holch10, M. Holler25, D. Horns26,
Zhiqiu Huang2, M. Jamrozy27, F. Jankowsky4, I. Jung-Richardt22, E. Kasai7, K.
Katarzyński28, U. Katz22, D. Khangulyan29, B. Khélifi12, S. Klepser10, W.
Kluźniak30, Nu. Komin11, R. Konno10, K. Kosack18, D. Kostunin10, S. Le Stum5,
A. Lemière12, M. Lemoine-Goumard24, J.-P. Lenain16, F. Leuschner14, T.
Lohse31, A. Luashvili15, I. Lypova4, J. Mackey1, D. Malyshev14, D. Malyshev22,
V. Marandon2, P. Marchegiani11, A. Marcowith23, G. Martí-Devesa25, R. Marx4,
G. Maurin17, M. Meyer26, A. Mitchell22,2, R. Moderski30, L. Mohrmann2, A.
Montanari18, E. Moulin18, J. Muller6, T. Murach10, K. Nakashima22, M. de
Naurois6, A. Nayerhoda9, J. Niemiec9, A. Priyana Noel27, P. O’Brien32, S.
Ohm10, L. Olivera-Nieto2, E. de Ona Wilhelmi10, M. Ostrowski27, S. Panny25, M.
Panter2, R.D. Parsons31, G. Peron2, S. Pita12, V. Poireau17, D.A. Prokhorov33,
H. Prokoph10, G. Pühlhofer14, M. Punch12,13, A. Quirrenbach4, P.
Reichherzer18, A. Reimer25, O. Reimer25, M. Renaud23, B. Reville2, F. Rieger2,
G. Rowell21, B. Rudak30, H. Rueda Ricarte18, E. Ruiz-Velasco2, V. Sahakian34,
S. Sailer2, H. Salzmann14, D.A. Sanchez17, A. Santangelo14, M. Sasaki22, J.
Schäfer22, F. Schüssler18, H.M. Schutte8, U. Schwanke31, M. Senniappan13,
J.N.S. Shapopi7, R. Simoni33, A. Sinha23, H. Sol15, A. Specovius22, S.
Spencer19, Ł. Stawarz27, S. Steinmassl2, C. Steppa20, T. Takahashi35, T.
Tanaka36, A.M. Taylor10, R. Terrier12, C. Thorpe-Morgan14, M. Tsirou2, N.
Tsuji37, R. Tuffs2, Y. Uchiyama29, T. Unbehaun22, C. van Eldik22, B. van
Soelen38, J. Veh22, C. Venter8, J. Vink33, S.J. Wagner4, F. Werner2, R.
White2, A. Wierzcholska9, Yu Wun Wong22, A. Yusafzai22, M. Zacharias15,8, D.
Zargaryan1,3, A.A. Zdziarski30, A. Zech15, S.J. Zhu10, S. Zouari12, N.
Żywucka8
1\. Dublin Institute for Advanced Studies, 31 Fitzwilliam Place, Dublin 2,
Ireland
2\. Max-Planck-Institut für Kernphysik, P.O. Box 103980, D 69029 Heidelberg,
Germany
3\. High Energy Astrophysics Laboratory, Russian-Armenian University (RAU),
123 Hovsep Emin St Yerevan 0051, Armenia
4\. Landessternwarte, Universität Heidelberg, Königstuhl, D 69117 Heidelberg,
Germany
5\. Aix Marseille Université, Centre national de la recherche scientifique
(CNRS)/Institut National de Physique Nucléaire et Physique des Particules
(IN2P3), Centre de Physique des Particules de Marseille (CPPM), Marseille,
France
6\. Laboratoire Leprince-Ringuet, École Polytechnique, CNRS, Institut
Polytechnique de Paris, F-91128 Palaiseau, France
7\. University of Namibia, Department of Physics, Private Bag 13301, Windhoek
10005, Namibia
8\. Centre for Space Research, North-West University, Potchefstroom 2520,
South Africa
9\. Instytut Fizyki Ja̧drowej Polskiej Akademii Nauk (PAN), ul. Radzikowskiego
152, 31-342 Kraków, Poland
10\. Deutsches Elektronen-Synchrotron DESY, Platanenallee 6, 15738, Germany
11\. School of Physics, University of the Witwatersrand, 1 Jan Smuts Avenue,
Braamfontein, Johannesburg, 2050 South Africa
12\. Université de Paris, CNRS, Astroparticule et Cosmologie, F-75013 Paris,
France
13\. Department of Physics and Electrical Engineering, Linnaeus University,
351 95 Växjö, Sweden
14\. Institut für Astronomie und Astrophysik, Universität Tübingen, Sand 1, D
72076 Tübingen, Germany
15\. Laboratoire Univers et Théories, Observatoire de Paris, Université PSL,
CNRS, Université de Paris, 92190 Meudon, France
16\. Sorbonne Université, Université Paris Diderot, Sorbonne Paris Cité,
CNRS/IN2P3, Laboratoire de Physique Nucléaire et de Hautes Energies (LPNHE), 4
Place Jussieu, F-75252 Paris, France
17\. Université Savoie Mont Blanc, CNRS, Laboratoire d’Annecy de Physique des
Particules - IN2P3, 74000 Annecy, France
18\. Institute for Research on the Fundamental Laws of the Universe (IRFU),
Commisariat à l’énergie atomique (CEA), Université Paris-Saclay, F-91191 Gif-
sur-Yvette, France
19\. University of Oxford, Department of Physics, Denys Wilkinson Building,
Keble Road, Oxford OX1 3RH, UK
20\. Institut für Physik und Astronomie, Universität Potsdam, Karl-Liebknecht-
Strasse 24/25, D 14476 Potsdam, Germany
21\. School of Physical Sciences, University of Adelaide, Adelaide 5005,
Australia
22\. Friedrich-Alexander-Universität Erlangen-Nürnberg, Erlangen Centre for
Astroparticle Physics, Erwin-Rommel-Str. 1, D 91058 Erlangen, Germany
23\. Laboratoire Univers et Particules de Montpellier, Université Montpellier,
CNRS/IN2P3, CC 72, Place Eugène Bataillon, F-34095 Montpellier Cedex 5, France
24\. Université Bordeaux, CNRS, Laboratoire de Physique des Deux Infinis
(LP2i), Bordeaux, Joint Research Unit (UMR 5797), F-33170 Gradignan, France
25\. Institut für Astro- und Teilchenphysik, Leopold-Franzens-Universität
Innsbruck, A-6020 Innsbruck, Austria
26\. Universität Hamburg, Institut für Experimentalphysik, Luruper Chaussee
149, D 22761 Hamburg, Germany
27\. Obserwatorium Astronomiczne, Uniwersytet Jagielloński, ul. Orla 171,
30-244 Kraków, Poland
28\. Institute of Astronomy, Faculty of Physics, Astronomy and Informatics,
Nicolaus Copernicus University, Grudziadzka 5, 87-100 Torun, Poland
29\. Department of Physics, Rikkyo University, 3-34-1 Nishi-Ikebukuro,
Toshima-ku, Tokyo 171-8501, Japan
30\. Nicolaus Copernicus Astronomical Center, Polish Academy of Sciences, ul.
Bartycka 18, 00-716 Warsaw, Poland
31\. Institut für Physik, Humboldt-Universität zu Berlin, Newtonstr. 15, D
12489 Berlin, Germany
32\. Department of Physics and Astronomy, The University of Leicester,
University Road, Leicester, LE1 7RH, United Kingdom
33\. Gravitation and Astroparticle Physics at the University of Amsterdam
(GRAPPA), Anton Pannekoek Institute for Astronomy, University of Amsterdam,
Science Park 904, 1098 XH Amsterdam, The Netherlands
34\. Yerevan Physics Institute, 2 Alikhanian Brothers St., 375036 Yerevan,
Armenia
35\. Kavli Institute for the Physics and Mathematics of the Universe (World
Premier International Research Center Initiative (WPI)), The University of
Tokyo Institutes for Advanced Study (UTIAS), The University of Tokyo, 5-1-5
Kashiwa-no-Ha, Kashiwa, Chiba, 277-8583, Japan
36\. Department of Physics, Konan University, 8-9-1 Okamoto, Higashinada,
Kobe, Hyogo 658-8501, Japan
37\. Institute of Physical and Chemical Research (RIKEN), 2-1 Hirosawa, Wako,
Saitama 351-0198, Japan
38\. Department of Physics, University of the Free State, PO Box 339,
Bloemfontein 9300, South Africa

## Materials and methods

### H.E.S.S. Observations and Data Analysis

H.E.S.S. is an array of five Imaging Atmospheric Cherenkov Telescopes located
in the Khomas Highland of Namibia, sensitive to VHE gamma rays in the energy
range between a few tens of GeV and a few tens of TeV (?, ?). H.E.S.S.
consists of four 12 m-diameter Davies-Cotton telescopes (CT1-4), equipped with
Cherenkov cameras (?) and have a field-of-view (FoV) of 5∘. The largest
telescope (CT5) is located in the middle of the CT1-4 array, has a 28
m-diameter mirror and is equipped with a Cherenkov camera with 3.4∘ FoV (?).
Observations of classical novae by H.E.S.S. are triggered by external alerts
and follow-up observations are decided on a case-by-case basis taking into
account if a transient signal is reported by _Fermi_ -LAT, if the ejecta
velocity is high ($\gtrsim 1500\,\mathrm{km\,s}^{-1}$), and if the peak
optical magnitude is bright (mV $\lesssim 9$). H.E.S.S. follows up on
$\sim$two nova candidates on average per year, subject to them satisfying at
least one of the aforementioned conditions. Initial reports of an outburst of
RS Oph on $8^{\rm th}-9^{\rm th}$ August 2021 satisfied all three conditions,
with a 6 $\sigma$ detection by _Fermi_ -LAT (?), ejecta velocity of $\gtrsim
2600\,\mathrm{km\,s}^{-1}$ (?, ?) and optical magnitude of mV $\sim 5.0$,
prompting H.E.S.S. observations with high priority. RS Oph was the first nova
visible to H.E.S.S. to meet all three conditions at the same time.

Follow-up observations of RS Oph with H.E.S.S. started on August 9th, 18:17:40
(Coordinated Universal Time, UTC) and concluded on September 7th, 19:47:31
(UTC). H.E.S.S. accumulated a total of 24.3 hours of observations over the
first five nights after the explosion. After the full moon period, during
which H.E.S.S. observations are paused, H.E.S.S. continued RS Oph observations
for a further 32.9 hours. The observations were conducted in the standard
“wobble mode” where telescopes are pointed at alternate offsets of 0.5∘ from
the position of RS Oph (right ascension 17h50min13.16s, declination
-06∘42′28.5′′ (J2000 equinox)). To achieve a low energy threshold, we only
selected data for analysis with zenith angle $\lesssim 35^{\circ}$, with a
resulting average zenith angle of $22^{\circ}$.

Throughout the first five nights, the data set is subject to strongly variable
observing conditions. Due to increasing levels of moonlight, a spectral
analysis of CT5 data is performed for the first three nights only. On the
second and third nights, atmospheric conditions were poor due to a higher-
than-usual aerosol content in the atmosphere, whilst atmospheric conditions
were good on the first, fourth and fifth night. We applied corrections to
account for the night-to-night variations in atmospheric conditions (see
below). Furthermore, we conducted observations under low to moderate moonlight
during nights 2 to 5, which increased the noise from night-sky-background
(NSB) photons by a factor of 1.5 in night 4 and 2.5 in night 5 compared to
dark sky conditions (c.f. Table S1).

Night | $T_{\mathrm{obs}}$ | Livetime | Significance | Atmospheric transparency | NSB noise level | Telescopes
---|---|---|---|---|---|---
| (UTC) | (hours) | ($\sigma$) |  |  |
09 Aug. 2021 | 18:17:40 | 3.2 | 5.8 (6.4) | 0.90 | 1.0 | CT1-5
10 Aug. 2021 | 17:53:46 | 3.7 (2.8) | 9.0 (7.1) | 0.80 | 1.0 | CT1-5
11 Aug. 2021 | 17:44:08 | 3.7 | 9.8 (9.6) | 0.65 | 1.0 | CT1-5
12 Aug. 2021 | 18:17:12 | 2.3 | 13.6 | 1.00 | 1.5 | CT1-4
13 Aug. 2021 | 17:44:43 | 2.8 | 10.5 (9.4) | 1.10 | 2.5 | CT1-5
25 Aug. – 07 Sep. 2021 | 17:48:03; 19:47:31 | 14.6 (13.4) | 3.3 (2.3) | 0.96 | 1.0 | CT1-5
Table S1: Summary of H.E.S.S. observations and data set used in the analysis.
The start time $T_{\mathrm{obs}}$ and duration of observations is given for
each night. Values are for the CT1-4 stereo analysis, except numbers in
brackets which denote CT5 mono analysis results. The atmospheric transparency
is quantified using the Cherenkov transparency coefficient (?). The NSB noise
level is given relative to the nominal Galactic level experienced by the
telescopes. Only observations taken under good atmospheric conditions and low
zenith angle $\lesssim 35^{\circ}$ were used in the data set from 25 Aug. 2021
to 07 Sep. 2021, leaving $\sim$15 hours for analysis out of a total of 32.9
hours of observations.

We performed the analysis of the H.E.S.S. data using both events recorded by
CT1-4 in stereoscopic (stereo) mode, and in a monoscopic (mono) analysis of
CT5. The ring-background method (?) was employed for the extraction of skymaps
for both analyses. The resulting significance maps are shown in Figure 1 for
CT1-4. Background events were classified and rejected using multivariate
analysis techniques (?, ?, ?). The reconstruction of shower properties was
performed using a pixel-based maximum likelihood technique (CT1-4 stereo, (?))
and a multivariate-based reconstruction (CT5 mono, (?, ?)). The spectral
reconstruction and model fitting were performed with a likelihood-based
procedure, using the Gammapy software package (?, ?, ?) version 0.18.2. In the
case of the CT5 mono analysis, the reflected-region (?) method was used to
define background control regions and to derive spectral results.

Significant VHE gamma-ray emission was detected by H.E.S.S. from the direction
of RS Oph on each of the five nights from $9^{\rm th}-13^{\rm th}$ August 2021
(Table S1) in an online analysis and later confirmed offline by two
independent analysis chains (?, ?, ?). The cross-check analysis employs an
independent calibration, reconstruction and background suppression (?). On
night 4, CT5 did not participate in observations for technical reasons, but it
detected RS Oph in VHE gamma rays in the other four nights.

RS Oph was observable again with H.E.S.S. under good low mooonlight conditions
on August $25^{\rm{th}}$ and observations continued as long as the source was
visible to H.E.S.S. until September 7th. The analysis of these observations
revealed a VHE gamma-ray signal from RS Oph at the $3-4\sigma$-level in the
main CT1-4 analysis as well as in the cross-check analysis. No significant
emission was detected with the CT5 mono analysis in this late phase,
consistent with its somewhat lower sensitivity. The corresponding flux upper
limit is $F_{\mathrm{UL}}(E>250\,\mathrm{GeV})=3.7\times
10^{-12}\,\mathrm{cm}^{-2}\,\mathrm{s}^{-1}$. Table S1 summarizes the H.E.S.S.
observations considered in our analysis below.

The quality of the atmospheric conditions is quantified using the Cherenkov
transparency coefficient (?) with lower values corresponding to lower
transmission of Cherenkov light through the atmosphere, as reported in Table
S1 for the RS Oph observations. To correct for the different atmospheric
conditions and resulting light-yield of the telescopes, the energy scale of
the energy migration matrix as well as the effective areas have been scaled
for individual observations with the atmospheric transparency coefficient for
the CT1-4 stereo analysis. To validate this correction, the reconstructed Crab
Nebula spectrum (after scaling the energy scale of the instrument response
functions (IRFs) with the transparency coefficient) is compared to the
published H.E.S.S. spectrum, using the same IRFs as employed for the RS Oph
observations. A good agreement between the spectral parameters of the
corrected Crab Nebula spectrum and a previously published spectrum (?) is
found for a range of atmospheric conditions, with the flux matching to within
10% (see Figure S1.)

Figure S1: Atmospheric transparency correction of Crab Nebula observations.
Illustration of the effect of atmospheric transparency correction to Crab
Nebula observations taken during the same observing season as RS Oph for a
range of atmospheric transparency conditions. The top panel illustrates the
impact of the correction on the normalization $\phi_{0}$ of the reconstructed
flux, $\phi=\phi_{0}\cdot(E/1\,{\rm TeV})^{-\Gamma}$. The bottom panel shows
the reconstructed index $\Gamma$ without and with correction. Vertical blue
lines indicate the atmospheric transparency conditions corresponding to the
five nights of RS Oph observations.

For the CT5 mono analysis, two different sets of IRFs tailored to the
different atmospheric conditions have been used. A comparison study to the
same Crab Nebula data set as described above was performed for the CT5 mono
analysis, confirming that also here the Crab Nebula flux and gamma-ray
spectral index is reconstructed correctly for the range of atmospheric
transparencies experienced during the RS Oph observations (Figure S1).

To test for spectral variation during the outburst, the H.E.S.S. data were
analysed in nightly bins for the five nights with significant detection (see
Table S1) and fitted with a power-law model. Best-fitting parameters are given
in Table S2 for the first three nights using the CT5 mono analysis and for all
five nights using the CT1-4 stereo analysis. Nightly spectral results are
consistent between the two configurations.

| Data set | $\phi_{0}$ | $E_{0}$ | Index $\Gamma$
---|---|---|---|---
|  | [$10^{-11}$ TeV-1 cm-2 s-1] | [TeV] |
mono | 09 Aug. 2021 | $14.9\pm(2.7)_{\text{stat.}}\pm(3.0)_{\text{syst.}}$ | 0.18 | $3.22\pm(0.38)_{\text{stat.}}\pm(0.20)_{\text{syst.}}$
| 10 Aug. 2021 | $25.2\pm(4.7)_{\text{stat.}}\pm(5.0)_{\text{syst.}}$ | 0.18 | $4.01\pm(0.48)_{\text{stat.}}\pm(0.20)_{\text{syst.}}$
| 11 Aug. 2021 | $28.5\pm(3.3)_{\text{stat.}}\pm(5.7)_{\text{syst.}}$ | 0.18 | $3.15\pm(0.23)_{\text{stat.}}\pm(0.20)_{\text{syst.}}$
stereo | 09 Aug. 2021 | $0.91\pm(0.28)_{\text{stat.}}\pm(0.14)_{\text{syst.}}$ | 0.35 | $4.24\pm(0.75)_{\text{stat.}}\pm(0.15)_{\text{syst.}}$
| 10 Aug. 2021 | $1.90\pm(0.32)_{\text{stat.}}\pm(0.38)_{\text{syst.}}$ | 0.35 | $3.32\pm(0.30)_{\text{stat.}}\pm(0.15)_{\text{syst.}}$
| 11 Aug. 2021 | $3.57\pm(0.54)_{\text{stat.}}\pm(0.54)_{\text{syst.}}$ | 0.35 | $4.08\pm(0.42)_{\text{stat.}}\pm(0.20)_{\text{syst.}}$
| 12 Aug. 2021 | $3.00\pm(0.33)_{\text{stat.}}\pm(0.45)_{\text{syst.}}$ | 0.35 | $3.27\pm(0.21)_{\text{stat.}}\pm(0.15)_{\text{syst.}}$
| 13 Aug. 2021 | $1.77\pm(0.25)_{\text{stat.}}\pm(0.35)_{\text{syst.}}$ | 0.35 | $3.24\pm(0.24)_{\text{stat.}}\pm(0.15)_{\text{syst.}}$
stereo | 25 Aug. 2021 - 07 Sep. 2021 | $0.238\pm(0.080)_{\text{stat.}}\pm(0.036)_{\text{syst.}}$ | 0.35 | $3.33\pm(0.45)_{\text{stat.}}\pm(0.15)_{\text{syst.}}$
Table S2:  Nightly spectral measurements of RS Oph. Best-fitting nightly
spectral parameters from H.E.S.S. assuming a power-law model of the form
$\phi_{0}\left(\frac{E}{E_{0}}\right)^{-\Gamma}$, with derived systematic
uncertainties. $E_{0}$ is the reference energy, $\phi_{0}$ is the amplitude at
the reference energy and $\Gamma$ is the spectral index.

#### Systematic Uncertainties

Multiple sources of systematic uncertainties contribute to affect the spectral
measurements of the CT1-4 stereo and the CT5 mono analyses. We summarise the
various sources and their estimated influence on the reconstructed flux and
gamma-ray spectral index. Systematic uncertainties stemming from Monte-Carlo
extensive air shower hadronic interaction models, broken pixels of the
Cherenkov cameras, and the live time of the data set are sub-dominant for both
analyses, but have been taken into account (?).

For the CT1-4 analysis, two different sets of cuts to select gamma-ray-like
events have been applied to the night 1-5 data set, showing a systematic
difference in the reconstructed flux normalisation and spectral index of
$\sim$5% and 0.1, respectively. The tests conducted using Crab Nebula
observations under varying atmospheric conditions capture run-by-run
differences in atmospheric conditions after correction, differences in the
assumed and actual optical telescopes efficiency, and the scaling to correct
for low atmospheric transparency itself. This results in a systematic error of
the flux normalisation of 10%. The impact on the reconstructed spectral index
varies from 0.05 (good, medium) to 0.15 (poor) atmospheric transparency. An
additional 10% systematic uncertainty of the flux normalisation is assumed for
nights 2 and 5, when the Cherenkov cameras were operated at higher camera
trigger settings and during moonlight.

The uncertainty in the CT5 mono analysis is derived from the following
aspects. Following (?), an imperfect description of the background acceptance
leads to a systematic uncertainty of $\sim$15% on the flux normalisation and
0.15 on the spectral index for the CT5 mono analysis. The Crab Nebula spectral
comparison under comparable atmospheric conditions yields a systematic error
on the flux normalisation of 5%, 10% and 15% for good, medium and poor
atmospheric conditions, respectively. The uncertainty on the spectral index is
at a level of 0.15 for the different transparencies.

A comparison to the intrinsic spectrum of the Crab Nebula (?), shows that the
overall energy scale uncertainty of the H.E.S.S. observations is at the 15%
level. Table S2 includes the resulting systematic errors of the RS Oph
spectral measurements for nights 1 to 5 and the CT1-4 stereo as well as the
CT5 mono analysis. The systematic error for the late follow-up observations is
derived in the same way to that described above.

###  _Fermi_ -LAT Data Analysis

We analyse _Fermi_ -LAT data to investigate spectral evolution in bins of 24
hours, centred on 20:00 UTC over the duration of the H.E.S.S. observations. We
use _Fermi_ -LAT (P8R3, pass 8 release 3 (?)) data spanning from 9th August
2021 to 7th September 2021, in the energy range from 60 MeV to 500 GeV. We
retrieve the data from a region of interest (ROI) defined by a radius of
$15^{\circ}$ around the position of RS Oph. The analysis of the LAT data
described above is performed using the fermipy python package (version
0.18.0), based on the Fermi Science Tools (version 2.0.8) (?, ?). We analyzed
SOURCE class events with a maximum zenith angle of $90^{\circ}$ to eliminate
Earth limb events. To derive the daily spectrum, the data are binned in 8
energy bins per decade and spatial bins of $0.1^{\circ}$ size. (The re-binned
points in Figure 3 in the main article were determined only for plotting
purposes, after the best-fitting parameters had been found.) The response of
the LAT instrument is evaluated with the IRFs, version P8R3_SOURCE_V2) and we
include in the model of the region all the LAT sources listed in the _Fermi_
-LAT Fourth Source Catalog (4FGL, (?)) in a radius of $20^{\circ}$ around the
position of the nova. The contributions from Galactic and extra-galactic
diffuse gamma-rays are described using the Galactic (gll_iem_v07) and
isotropic (iso_P8R3_SOURCE_V2_v1) diffuse emission models. The models are
available from the Fermi Science Support Center (FSSC) (?). The energy
dispersion correction is applied to all individual models describing the
diffuse and discrete gamma-ray emission in the ROI, except for the isotropic
diffuse emission model.

Spectral results obtained using 6 hour bins contemporaneous to the H.E.S.S.
observations are found to be compatible with the results of the 24 hour bins,
albeit with higher uncertainty. Additionally, we evaluate the systematic
errors for the daily observations by applying the same analysis to a pre-flare
24 hour observation window. No significant gamma-ray signal is found at the
position of the RS Oph nova in the pre-flare time interval, with a signal
first detected during the optical rise time. This demonstrates that neither
systematic effects in the Galactic diffuse emission model nor source confusion
impact the _Fermi_ -LAT flux measurement.

To compute the light curve, the time range is split into a series of bins of
varying duration. For the first five nights of the H.E.S.S. observations,
contemporaneous 6 h bins are used (and 18 h bins in between these times).
Outside of these five nights, we use bins of 24 h duration (or multiples
thereof) as appropriate, given the statistics. To calculate the energy flux
for a given time window, we model the sources’ gamma-ray emission with a log-
parabola spectral function of the form
$\phi(E)=\phi_{0}\left(\frac{E}{E_{0}}\right)^{-\alpha-\beta\log\left(E/E_{0}\right)}$,
with parameters set to a spectral index $\alpha=2.0$, curvature$\beta=0.13$
and reference energy $E_{0}=1$ GeV, respectively. The flux normalisation
$\phi_{0}$ is left free to vary to obtain the integrated flux in the 60 MeV to
500 GeV energy range.

We use the derived _Fermi_ -LAT spectral points for each observation period to
perform joint spectral fits with the H.E.S.S. mono and stereo data, using the
Gammapy software package (version 0.18.2) (?). The different data sets were
fitted simultaneously using a log-parabola spectral model; the best fitting
parameters for each night are quoted in Table S3. The general trend is for the
flux normalisation to decrease and the parabola to widen, over the course of
the five nights.

Data set | $\phi_{0}$ | $E_{0}$ | $\alpha$ | $\beta$
---|---|---|---|---
| [$10^{-4}$ TeV-1 cm-2 s-1] | [GeV] |  |
09 Aug. 2021 | $7.06\pm 0.58$ | 1.0 | $1.98\pm 0.04$ | $0.19\pm 0.01$
10 Aug. 2021 | $4.27\pm 0.43$ | 1.0 | $2.12\pm 0.05$ | $0.13\pm 0.01$
11 Aug. 2021 | $3.69\pm 0.42$ | 1.0 | $2.01\pm 0.06$ | $0.13\pm 0.01$
12 Aug. 2021 | $1.79\pm 0.32$ | 1.0 | $2.02\pm 0.09$ | $0.11\pm 0.02$
13 Aug. 2021 | $1.94\pm 0.32$ | 1.0 | $2.05\pm 0.07$ | $0.12\pm 0.02$
Table S3: Nightly log-parabola fits for RS Oph. Best-fitting spectral
parameters from a joint fitting of the _Fermi_ -LAT and H.E.S.S. data with a
log-parabola model. $E_{0}$ is the reference energy, $\phi_{0}$ is the
amplitude at the reference energy, $\alpha$ is the spectral index and $\beta$
is the curvature.

## Supplementary Text

### Modelling of gamma-ray emission

Recurrent novae are complex hydrodynamic processes accompanied by the
formation of multiple shock waves. The shocks are either internal, i.e.,
formed by collision of different parts of ejecta, or external, i.e., formed in
the circumbinary medium. Both internal and external shocks might be
responsible for the production of gamma-ray emission. While internal shocks
remain a hypothesis (see, however, (?)), external shock formation is
unavoidable and the emission produced at the external shock is expected to be
long-lasting. This was evident from the multi-wavelength observations of the
2006 RS Oph nova event (?, ?, ?, ?, ?). We therefore focus on the external
shock component. Similarly to supernova explosions, the external shocks of
recurrent novae have different phases of evolution (?). Initially, when the
circumbinary medium cannot substantially influence the explosion dynamics, the
blast wave expands with approximately constant velocity, the so-called ejecta-
dominated phase. Once a sufficient amount of material is swept up, the shock
decelerates, transitioning to the Sedov-Taylor phase. The shock deceleration
rate further increases in the radiative phase. For a point-like explosion in a
radial wind with a mass density ($\rho$) profile $\rho(r)\propto r^{-2}$, the
shock radius ($R_{\rm sh}$) satisfies $R_{\rm sh}\propto t^{2/3}$ during the
Sedov-Taylor phase (?). When cooling dominates the shock enters the radiative
phase, with $R_{\rm sh}\propto t^{1/2}$. Observations of the previous 2006
outburst of RS Oph indicate that after day $\approx 6$ the spectroscopic data
is consistent with a shock speed $u_{sh}\propto t^{-0.6}$ (?, ?). While
$u_{sh}\propto t^{-0.3}$ over the first 3 weeks, it is inconsistent with the
expected X-ray flux (?), indicating a self-similar solution does not apply.
Follow-up imaging of the expanding shock from the 2006 outburst (?) showed
that the polar expansion continued at about $5\,500$ km s-1 for at least 155
days, apparently with minimal deceleration. This is consistent with 3D
hydrodynamical simulations showing late-time expansion at $>\,4\,000$ km s-1
in the polar directions (?). This points to the potential complexity of the
source geometry and difficulty in interpreting the measurements in the context
of a model with a single shock velocity. The mid-latitude regions with slower
shocks, denser gas and stronger radiative cooling, could dominate the spectral
measurements. Assuming a similar evolutionary profile for the 2021 outburst,
the source was detected with H.E.S.S. over the first five days primarily
during the ejecta-dominated phase, during which the velocity is in the range
$(2-4)\times 10^{3}$ km s-1. The velocity may have decreased by the time
H.E.S.S. observations recommenced after the bright moonlight break. Optical
spectroscopy observations suggest that on day 2 the shock was faster than
found in earlier measurements made on day 1 (?, ?). This we attribute to the
shock propagation in an inhomogeneous stellar wind, as found in previous 3D
simulations (?).

Figure S2: Schematic of external shock model. The explosion originates near
the surface of the white dwarf (WD). Within one day, the shock is expanding as
a bipolar blast wave moving orthogonal to the accretion disk, into the wind of
the red giant (RG). The colour gradient of the shock indicates the expansion
velocity, whilst the red gradient of the surrounding medium indicates the
density of the RG wind. Material internal to the shock is shown in blue.
Figure is not to scale.

##### Shock conditions:

We assume that the explosion origin is located at the white dwarf star, which
moves around a red giant star on an orbit of radius $r_{\rm orb}\approx
1.48\rm\,au$ (?). The accretion onto the white dwarf proceeds through a dense
disk in the orbital plane of the system. This results in the quasi-spherical
blast wave initially moving faster in directions perpendicular to the orbital
plane, as two sections of a spherical shock, as illustrated in Figure S2. This
picture is consistent with observations of the previous nova event of RS Oph
(?). The period relevant for the H.E.S.S. observations likely covers the
transition of the outflow expansion from the bipolar to quasi-spherical
regimes (and also from the ejecta-dominated to the Sedov-Taylor phases), thus
self-similar solutions (see, e.g. (?, ?)) cannot describe the entirety of this
period consistently.

We employ a simple numerical model that describes the motion of the ejecta
(namely distance covered by the ejecta $R_{\mathrm{ej}}$, and its speed
$v_{\mathrm{ej}}$) and the dynamics of the explosion’s forward shock (its
radius $R_{\mathrm{sh}}$, and speed $u_{\mathrm{sh}}$). These quantities are
obtained under the following assumptions: (i) the density and pressure across
the shocked region are constant, and determined by the Rankine–Hugoniot
conditions (?); (ii) the upstream medium density is determined by the distance
covered by the shock as

$\rho_{\mathrm{up}}=\frac{\dot{M}}{4\pi v_{\rm w}\left(r_{\rm
orb}^{2}+R_{\mathrm{sh}}^{2}\right)}\,.$ (S3)

Here $\dot{M}$ is the mass-loss rate of the red giant star and $v_{\rm w}$ is
its wind speed, which depends on the distance as

$v_{\rm w}=v_{\infty}\left(1-\frac{R_{*}}{\sqrt{r_{\rm
orb}^{2}+R_{\mathrm{sh}}^{2}}}\right)^{\beta_{\rm w}}\,,$ (S4)

with $v_{\infty}$ the terminal wind speed. We take the value $R_{*}\approx
0.35\rm\,au$ for the radius of the red giant (?) and adopt for the exponent of
the wind acceleration $\beta_{\rm w}=3$ (?). For the red giant star in RS Oph
the wind parameter $\dot{M}/v_{\infty}$ is estimated to be in the range
$(5.7~{}-~{}7.4)\times 10^{11}~{}{\rm kg~{}m}^{-1}$ (?). We therefore adopt a
fiducial value of $\dot{M}/v_{\infty}=6\times 10^{11}~{}{\rm kg~{}m}^{-1}$.
For a typical terminal wind velocity of $v_{\infty}=(10-30)\rm\,km\,s^{-1}$
for red giants (?), this corresponds to a mass-loss rate of
$\dot{M}\approx(1-3)\times 10^{-7}M_{\odot}\rm\,yr^{-1}$ (where $M_{\odot}$ is
the solar mass).

For a strong shock in single-atom polytropic gas, one can estimate the shocked
gas ($\rho_{\mathrm{ds}}$) mass-density as

$\rho_{\mathrm{ds}}\approx 10^{-14}\left(\frac{\sqrt{r_{\rm
orb}^{2}+R_{\mathrm{sh}}^{2}}}{1\rm\,au}\right)^{-2}~{}{\rm g~{}cm}^{-3}\,,$
(S5)

where we set $v_{\rm w}=20\,\rm km\,s^{-1}$. The numerical value outside the
parentheses in the shocked gas density expression equates to a proton number
density of $n_{\mathrm{ds}}\approx 5\times 10^{9}{\rm cm}^{-3}$. The
corresponding pressure $p_{\mathrm{ds}}$ in the shocked region

$p_{\mathrm{ds}}\approx\frac{3}{16}\rho_{\mathrm{ds}}u_{\mathrm{sh}}^{2}\,,$
(S6)

decelerates the ejecta according to

$m_{\mathrm{ej}}\frac{dv_{\mathrm{ej}}}{dt}=-4\pi
R_{\mathrm{ej}}^{2}p_{\mathrm{ds}}\,.$ (S7)

Here $m_{\mathrm{ej}}$ is the ejecta mass, $v_{\mathrm{ej}}$ is ejecta speed
(initially the ejecta is assumed to move with $v_{\mathrm{ej},0}$), and
$R_{\mathrm{ej}}$ is the distance traveled by the ejecta:

$R_{\mathrm{ej}}(t)=\int\limits_{0}^{t}dt^{\prime}v_{\mathrm{ej}}(t^{\prime})\,.$
(S8)

The radius of the forward shock $R_{\mathrm{sh}}$ is determined by the
condition that the shocked region density is equal to the value given by Eq.
(S5):

$\rho_{\mathrm{ds}}=\frac{M_{\mathrm{ds}}}{4\pi/3\left(R_{\mathrm{sh}}^{3}-R_{\mathrm{ej}}^{3}\right)}\,,$
(S9)

where the total mass of the shocked region is given by

$M_{\mathrm{ds}}(t)=4\pi\int\limits_{0}^{R_{\mathrm{sh}}(t)}dr^{\prime}r^{\prime}{}^{2}\rho_{\mathrm{up}}(r^{\prime})\,.$
(S10)

We solve the above system of equations using an iterative method and obtain
the time dependence of the forward shock speed, shocked region density, and
ejecta speed. The kinematics of the forward shock and ejecta are shown in Fig.
S3, including the shock/ejecta speed (Fig. S3A) and the path covered by the
shock (Fig. S3B). For the considered geometry the forward shock propagates in
a medium with a constant density during the first $\sim 12$ hours, Fig. S3B.
Afterwards, the upstream medium is characterized by a decreasing density,
which results in a period of acceleration of the forward shock. The forward
shock speed achieves its maximum value after approximately 3 days. The shock
subsequently starts to decelerate, as expected from self-similar solutions in
a medium with the density decreasing as $\rho_{\mathrm{up}}\propto r^{-2}$
(?).

Figure S3: Temporal shock evolution. A: Model shock speed as a function of
time since the explosion. Markers correspond to the expansion velocity
estimated from the absorption center and the blue absorption edge of the
$H\alpha$ and $H\beta$ P Cyg profiles obtained with Varese (?, ?) and SALT (?)
spectroscopy of RS Oph. The shock speed is obtained in the framework of a
single-zone model for initial ejecta velocity
$v_{\mathrm{ej},0}=3\,000\rm\,km\,s^{-1}$. B: distance covered by the forward
shock and distance between the center of the RG star and forward shock.

For the shock speeds in Figure S3, the post-shock temperature is $>10^{8}$ K
and the only effective radiative cooling process of the thermal plasma is
bremsstrahlung. Fig. S4 shows that the cooling time of the post-shock gas
behind the forward shock (assuming the red giant wind is compressed by a
factor of 4) is much longer than the expansion time of the nova for the
duration of the H.E.S.S. observations, and so we expect the forward shock to
be adiabatic, in contrast to models of slower internal shocks in classical
novae (?).

Figure S4: Cooling times. Bremsstrahlung cooling time ($t_{\rm cool}$) of the
thermal plasma behind the forward shock as a function of the time since nova
eruption ($t_{\rm e}$) in our model (blue solid line), with the time elapsed
since eruption also plotted for reference (orange dashed line). The pre-shock
gas density, $n_{\rm H}$, is also shown (blue dotted line, right-hand
$y$-axis).

In the simplest treatment of diffusive shock acceleration at the external
shock, assuming Bohm scaling (?), the acceleration rate for relativistic
particles of charge $q$ is determined by the shock speed and the post-shock
mean magnetic field strength, $B_{\mathrm{sh}}$ as

$\left.\frac{dE}{dt}\right|_{\mathrm{acc}}=\frac{qB_{\mathrm{sh}}u_{\mathrm{sh}}^{2}}{\eta
c}\,,$ (S11)

where $\eta\geq 2\pi$ is a dimensionless parameter corresponding to the ratio
of the relativistic particle-gyrofrequency to the scattering rate. Its inverse
can be taken to characterize the acceleration efficiency and is a measure of
the energy in the turbulent magnetic field. If non-linear magnetic field
amplification occurs, $\eta$ may fall well below $2\pi$. This should however
be done with care, since magnetic field amplification is unlikely to occur on
all scales simultaneously, which is the essential principle of Bohm scaling
(?). Particle acceleration to higher energies can be saturated either by
escaping the acceleration region (confinement limited) or, in the case of
electrons, when inverse Compton (IC) or synchrotron energy losses overtake the
acceleration (cooling limited). The mean magnetic field strength at the
forward shock is determined by the magnetic field of the red giant wind. We
follow the Parker (?) description of a stellar wind expanding radially
outwards from a rotating star, for which the magnetic field near the poles
approaches a split-monopole ($B\propto r^{-2}$) and near the equator is swept
into a Parker spiral with a toroidal-dominated field ($B\propto r^{-1}$).
Since the nova expansion is understood to be bipolar, and the rotation axis of
the red giant should be similar to the orbital axis of the binary system, we
adopt the split-monopole approximation for the magnetic field strength at
distance $R$ from the red giant:

$B_{\rm w}(R)=B_{*}\left(\frac{R}{R_{*}}\right)^{-2}\,,$ (S12)

where $B_{*}$ is the red giant surface magnetic field, typically of order
$1\,\mathrm{G}$ (?). We consider two cases, adopting $B_{*}=1$ and $10\,{\rm
G}$ (baseline and extreme values, respectively). The computed strength of the
magnetic field in the blast-wave downstream is shown in Figure S5A. The wind
should be kinetic-energy dominated at a radius of a few times $R_{*}$, so a
limit on the magnetic field (assuming wind expansion at $v_{\infty}$) is

$B_{\rm w}(R)<\frac{\sqrt{\dot{M}v_{\infty}}}{R}\approx
0.2\,\mathrm{G}\left(\frac{R}{1\,\mathrm{au}}\right)^{-1}\;.$ (S13)

This is consistent with our assumed surface field of order $(1\,\mathrm{G})$
at $R=R_{*}$. Substantially larger surface fields can be ruled out, as they
would imply magnetically dominated winds at large distances.

Figure S5: Magnetic and photon field evolution. A: Model magnetic fields as a
function of time since the explosion. B: Energy density in the optical photon
field at the forward shock and in the magnetic fields as a function of time
since the explosion.

##### Maximum Energy:

The VHE H.E.S.S. detection confirms the acceleration of particles to at least
several TeV during the ejecta-dominated phase. The Hillas limit (?) for a
particle of charge $q=Ze$ in the above radial field (i.e. ignoring strong
magnetic field amplification) is

$E_{\max}=q\frac{u_{\rm sh}}{c}R_{\rm sh}B(R_{\rm sh})\approx
10Z\left(\frac{u_{\rm sh}}{5\,000~{}{\rm
km~{}s^{-1}}}\right)\left(\frac{B_{*}}{1~{}{\rm G}}\right)\left(\frac{R_{\rm
sh}}{\rm au}\right)^{-1}~{}{\rm TeV}$ (S14)

assuming $R_{*}=0.35$ au. In the first week the shock advances approximately
$2-3$ au per day, such that strong magnetic field amplification above that of
the background value $B_{\rm w}$ at the forward shock is essential to account
for ongoing TeV particle acceleration $>3$ days after the explosion. To
facilitate acceleration, the magnetic fields are necessarily amplified
upstream of the shock, which can only be triggered by energetic particle
currents. For $B_{\rm w}\propto R^{-2}$ the magnetic energy density falls
rapidly, suggesting that only a small fraction of cosmic-ray energy needs to
be converted to magnetic energy, for example via the Bell instability (?), to
amplify the fields above background levels. This is an exact analogue of the
core-collapse supernova scenario (?), for which there is a formalism to
predict the maximum energy for a fast shock in a dense wind. We summarize that
calculation (?) below.

Neglecting shock curvature, and assuming isotropy of the particles, the
differential accelerating flux (i.e. the rate at which particles increase
their energy from $E$ to $E+dE$) at a shock is determined by the shock
compression (?)

$\psi(E)=\frac{\Delta u}{3}E~{}N_{E}$ (S15)

where $N_{E}(E)$ is the differential non-thermal particle density, and $\Delta
u=u_{\rm sh}-u_{\rm ds}$, i.e. the difference between the upstream and
downstream flow velocities. In the model (?), particles must generate their
own self-confining waves at the expense of an escaping flux. Setting the
escaping current density to the accelerating flux of the maximum energy
particles when the shock was at a radius $r$ i.e. $j_{\rm esc}(r)=q\psi(E_{\rm
max})$, and for simplicity assuming an $N_{E}\propto E^{-2}$ spectrum, it
follows that

$\frac{j_{\rm esc}E_{\rm max}}{q}=\xi_{\rm esc}\frac{1}{2}\rho_{\rm up}u_{\rm
sh}^{3}$ (S16)

where

$\xi_{\rm esc}=\frac{\Delta u}{u_{\rm sh}}\left(\frac{U_{\rm cr}}{\rho_{\rm
up}u_{\rm sh}^{2}/2}\right)\frac{1}{\log(E_{\rm max}/mc^{2})}$ (S17)

with $U_{\rm cr}$ the energy density of cosmic-rays above the rest mass
energy. The efficiency parameter $\xi_{\rm esc}$ corresponds to the fraction
of energy density flux processed by the shock that is lost to the upstream
escaping energetic particles, and since $E_{\rm max}\approx 10TeV$, we expect
$\log(E_{\rm max}/mc^{2})\approx 10$. For acceleration efficiencies of
${U_{\rm cr}}/(\frac{1}{2}\rho_{\rm up}u_{\rm sh}^{2})\approx 10\%$,
approximately $1\%$ of the energy density flux processed by the shock is lost
to the upstream. To determine if particles are confined, the fluctuations
driven by the non-resonant instability (?) need sufficient time to grow to a
level that can effectively scatter the highest energy particles. This requires
many growth times of the fastest growing instability. The current that passes
through a given fluid element at upstream radius $R$ is the integral of the
escaping current from the shock over the history of its expansion, diluted by
a factor $(r/R)^{2}$ (here we assume radial magnetic field and ignore
deceleration of the shock). Using the definition of the maximum growth rate
for cosmic-ray driven instabilities (?), and demanding at least $5$ growth
times, we find

$\int\Gamma_{\rm max}dt=\int\sqrt{\frac{\pi}{\rho_{\rm up}c^{2}}}j_{\rm
esc}dt>5$ (S18)

Using equation S16 to determine $j_{\rm esc}$, by the time the shock reaches a
fluid cell located at radius $R$, the confinement condition is

$\xi_{\rm esc}q\int_{R_{0}}^{R}\sqrt{\frac{\pi}{\rho_{\rm
up}c^{2}}}\frac{\rho_{\rm up}u_{\rm sh}^{2}}{E_{\rm max}}{r^{2}}dr>5R^{2}$
(S19)

where $R_{0}$ is the birth location of the shock. Differentiating both sides
with respect to $R$, we find

$E_{\rm max}<\frac{qc\xi_{\rm esc}}{20}\sqrt{\frac{\dot{M}}{v_{\rm
w}}}\left(\frac{u_{\rm sh}}{c}\right)^{2}\,,$ (S20)

which is independent of both the magnetic field strength and shock radius. Eq.
(S20) is the so-called confinement limit. This can also be used to estimate
the effective scattering magnetic field $B_{\rm eff}$ felt by the maximum
energy particles found by equating the above expression for the confinement
limit with the first equation in equation S14. This results in a magnetisation
parameter of

$\frac{B_{\rm eff}^{2}}{4\pi\rho_{\rm up}u_{\rm sh}^{2}}>\left(\frac{\xi_{\rm
esc}}{20}\right)^{2}\,.$ (S21)

This is only a representative value of the upstream field. Stronger magnetic
fields concentrated on scales much less than the gyroradius of the maximum
energy particles as measured in $B_{\rm eff}$ will not affect their
acceleration. Further amplification can also occur in the post-shock medium,
though this has only a minor effect on the acceleration rate, as the
acceleration time is dominated by the upstream residence time of the diffusing
particles. Because $E_{\rm max}$ depends only on the properties of the wind,
the acceleration time for particles with energy $E_{\rm max}$ is always
approximately equal to the instantaneous age of the shock.

Using the shock and wind parameters discussed above for RS Oph, the maximum
energy is several TeV. The prediction is consistent with the detection of TeV
gamma-rays from RS Oph. The confinement limit is the dominant constraint for
protons. Electrons must compete also with continuous IC, bremsstrahlung and
synchrotron energy losses, so their maximum energy may be limited further. In
the numerical leptonic model, the maximum electron energy does not exceed the
confinement limit, but this still requires an escaping flux with $\xi_{\rm
esc}\approx 0.01$.

##### Radiative Losses:

High-energy electrons can interact via different processes losing their energy
on time-scales comparable to the acceleration time. The typical cooling times
for the synchrotron $t_{\rm syn}$, IC (in the Thomson regime) $t_{\rm IC,Th}$,
and bremsstrahlung $t_{\rm br}$ processes are

$t_{\rm syn}\approx
400\left(\frac{E}{1\rm\,TeV}\right)^{-1}\left(\frac{B}{1\rm\,G}\right)^{-2}\rm\,s\,;$
(S22) $t_{\rm IC,Th}\approx
15\left(\frac{E}{1\rm\,TeV}\right)^{-1}\left(\frac{w_{\rm
ph}}{1\rm\,erg\,cm^{-3}}\right)^{-1}\rm\,s\,;$ (S23)

and

$t_{\rm br}\approx
10^{6}\left(\frac{n_{\mathrm{ds}}}{10^{9}\rm\,cm^{-3}}\right)^{-1}\rm s\,.$
(S24)

Here $w_{\rm ph}$ is the energy density of target photons and
$n_{\mathrm{ds}}$ is the number density of particles in the background plasma.
These estimates show that for the conditions expected in the production region
(see Fig. S5), the IC mechanism plays the dominant role. For VHE electrons,
Klein-Nishina suppression (?) decreases the rate of IC scattering, and
synchrotron cooling might become dominant in this range. We therefore consider
these two processes below.

As the radiative cooling becomes more important with increasing electron
energy, we are interested in accurately estimating its contribution for VHE
electrons. TeV electrons efficiently up-scatter target photons with energy
below the limit set by the Klein-Nishina cutoff, $<0.25\rm\,eV$. Thus, optical
and infra-red (IR) observations provide information about the intensity of the
soft photon fields, which serve as a target for IC losses and emission.
Assuming these photon fields result from thermal free-free emission, the
expected spectrum in this energy band can be approximated as

$S_{\nu}\propto\nu^{-2}B_{\nu}(T(t))\,,$ (S25)

where $B_{\nu}$ is the Planck function, $T$ is the (time-dependent)
temperature of the gas producing free-free emission, and the factor $\nu^{-2}$
appears due to free-free opacity. For the typical post-shock temperatures
$T>10^{8}\,{\rm K}$, we are interested in the lower energy part of the
spectrum, $h\nu<k_{B}T$, and the blackbody spectrum can be approximated by the
Rayleigh–Jeans law, $B_{\nu}\propto\nu^{2}$. Here $h$ and $k_{B}$ are Planck
and Boltzmann constants, respectively. The relevant part of the spectrum
should have a flat spectral energy distribution (SED). In our model
calculations we directly use the time-dependent information from the optical
B, V, R, and I bands (?), and below the I band we approximate the target
photon spectrum by a flat SED component. This approximation agrees with the IR
spectra measured from RS Oph after its eruption in 2006 (?). The time-
dependence of the photon target energy density is shown in Figure S5B.

Figure S6: Particle energy and luminosity evolution. A: Maximum energy of
accelerated electrons and protons at the forward shock. B: Energy crossing the
forward shock per time unit, $L_{\rm sh}$, and observed optical-IR luminosity
of the nova, $L_{\rm ph}$. The dip occurring around day 1 was in all optical
bands.

The evolution of the non-thermal particles follows the acceleration-losses
equation:

$\frac{dE}{dt}=\left.\frac{dE}{dt}\right|_{\mathrm{acc}}+\left.\frac{dE}{dt}\right|_{\mathrm{losses}}\,,$
(S26)

where the acceleration rate is given by equation (S11) and for electrons the
loss term accounts for both synchrotron and IC mechanisms,
$\left.\frac{dE}{dt}\right|_{\mathrm{losses}}=\left.\frac{dE}{dt}\right|_{\mathrm{syn}}+\left.\frac{dE}{dt}\right|_{\mathrm{ic}}$.
Relativistic bremsstrahlung is important at early times, but negligible
compared to synchrotron and IC by the time the H.E.S.S. observations
commenced. Similarly, for the expected densities, the energy losses for
protons are negligibly small. The maximum energy determined from the solution
to equation (S26) for the two cases considered is shown in Figure S6. For this
plot we adopted $\eta=2\pi$ both for electrons ($\eta_{e}$) and protons
($\eta_{p}$). In the case of protons, we assume the maximum energy follows the
confinement limit described above with $\xi_{\rm esc}=10^{-2}$. For protons
the maximum energy does not depend on the magnetic field strength (except for
a short time interval after the onset of the acceleration, which is not
distinguishable on the scale of the figure). The maximum energy grows rapidly
at early time due to the strong magnetic fields close to the star, and the
$\eta=2\pi$ example shown in Figure S6 likely over-estimates the acceleration
rate. Following the rapid growth, the maximum energy saturates at the
confinement limit and follows this value for later times. Thus non-linear
magnetic field amplification is implicit in the numerical model. The
transition from fast to slow cooling is also evident in the electron spectrum
.

The obtained maximum energy determines the spectral cut-off of the non-thermal
particles, which is injected downstream at each moment:

$\frac{dN_{\rm inj}}{dEdt}=\left\\{\begin{matrix}A_{e/p}E^{-\alpha_{e/p}}{\rm
exp}{\left[-\left(\frac{E}{E_{\mathrm{max},e/p}}\right)^{\beta_{e/p}}\right]}&E\geq
E_{\mathrm{min},e/p}\\\ 0&E<E_{\mathrm{min},e/p}\,.\end{matrix}\right.$ (S27)

Here $E_{\mathrm{max},e/p}$ is obtained by solving Eq. (S26) subject to the
limit in Eq. (S20), for electrons and protons; $E_{\mathrm{min},e/p}$ is a
free parameter set to $100\rm\,MeV$ for leptons and $2\rm\,GeV$ for protons to
avoid non-relativistic effects. The normalization coefficient $A$ is computed
by the condition that a fixed fraction, $\kappa_{e/p}$, of the energy flux
through the forward shock, $L_{\mathrm{sh}}$, is transferred to non-thermal
particles:

$A_{e/p}=\frac{\kappa_{e/p}L_{\mathrm{sh}}}{\int\limits_{E_{\mathrm{min},e/p}}^{\infty}dEE^{1-\alpha_{e/p}}{\rm
exp}{\left[-\frac{E}{E_{\mathrm{max},e/p}}\right]}}\,.$ (S28)

The energy flux through the shock is determined by the upstream medium density
and the forward shock speed:

$L_{\mathrm{sh}}=2\pi R^{2}\rho_{\mathrm{up}}u_{\mathrm{sh}}^{3}\,.$ (S29)

The model dependence of this parameter is shown in Figure S6B.

The spectrum of non-thermal particles, $n=dN/dE$, is described by the single-
zone evolution equation

$\frac{\partial n}{\partial t}+\frac{\partial}{\partial E}\left[\dot{E}_{\rm
losses}n\right]=\frac{dN_{\rm inj}}{dEdt}\,,$ (S30)

where the injection term is given by Eq. (S27) and the initial condition is
$n\big{|}_{t=0}=0$.

Figure S7: Single-zone particle spectra evolution. A: Electron spectra for
five instants in time post-outburst and two different values for the strength
of the red giant surface magnetic field. B: Spectra of non-thermal protons for
five instants in time (the impact of magnetic field strength is minor, due to
confinement limited saturation).

The non-thermal particle spectra expected in the source at days $1$, $2$, $3$,
$4$, and $5$ after the explosion are shown in Figure S7. In the case of non-
thermal protons the time evolution of the spectrum is simple. The impact of
cooling on the proton spectrum is minor and we expect an accumulation of non-
thermal protons accompanied by minor spectral evolution. In contrast, the
development of the non-thermal electrons is quite complex, and the spectral
properties are affected by the strength of the magnetic field. There are
features caused by the decrease in energy loss rates with radius: increasing
maximum energy and spectral hardening caused by the transition from the fast
cooling to slow cooling regimes. For the chosen parameters, we do not expect
acceleration of electrons to very-high energies during the first day after the
explosion for the case of weak external magnetic field. On the other hand, the
flux level in the VHE band is smaller than that at GeV energies thus a
sufficient VHE flux can be produced by particles with energy exceeding the
cutoff energy, in the tail of the distribution. Therefore for the radiation
calculations we consider the case with $B_{*}=1\,{\rm G}$, which as discussed
above is more consistent with the value expected at the RG surface.

In Figure S8 we compare the model IC spectra produced by the shock-accelerated
electrons to observational data obtained with _Fermi_ -LAT in the GeV band and
H.E.S.S. data obtained in the VHE regime. We set $\kappa_{e}=0.1$ and compute
the IC spectra for five different instants in time post-outburst: $t=1$, $2$,
$3$, $4$, and $5$ days after the explosion. The IC model under-predicts the
flux on night 1; although an increased IC flux for the first night can be
obtained by increasing the density of the wind in the model, this change would
also accelerate the evolution, with the wind deceleration commencing between
days 1 and 2 post explosion. A denser wind is therefore problematic due to the
rapid deceleration of the shock (incompatible with multi-wavelength
observations) and the need to introduce a suppression of the p-p channel to
maintain dominant IC (inconsistent with the theory of shock acceleration).

In Figure S9 we show the result of hadronic emission spectra (?) due to
collisions between the non-thermal protons and shock compressed external
medium for the same times, using $\kappa_{p}=0.5$. Here we assume that the
gamma-ray emission is produced in a region close to the external shock where
the density is twice higher than the averaged value (i.e., the filling factor
is $50\%$ ). The inverse of the acceleration-rate efficiency for electrons was
set to $\eta_{e}=10\pi$ and the maximum energy of protons was calculated for
$\xi_{\rm esc}=10^{-2}$ and $\eta_{p}=30\pi$.

Figure S8: Nightly Inverse Compton spectra. A-E IC emission of shock-
accelerated electrons for nights 1-5 respectively, together with spectral
measurements obtained with H.E.S.S. (filled points) and _Fermi_ -LAT (hollow
points). The model parameters are summarized in Table Radiative Losses:.
Figure S9: Nightly hadronic emission spectra. Same as Fig. S8, but for proton-
proton emission. The model parameters are summarized in Table Radiative
Losses:.

The local intense photon fields, especially during the first days after the
explosion, require consideration of $\gamma-\gamma$ attenuation. We provide
some general estimates, which allow us to define the photon frequencies and
epochs when $\gamma-\gamma$ absorption might be important. For a source with
luminosity $L$ and size $R_{\rm sh}$ we can put the following upper limit for
the optical depth, $\tau$, for photons with energy $\varepsilon$:

$\tau\ll\tau_{\rm max}=\frac{0.2\sigma_{\rm T}\varepsilon L}{4\pi R_{\rm
sh}c(3.5m_{e}^{2}c^{4})}=2\left(\frac{L}{10^{37}{\rm\,erg\,s}^{-1}}\right)\left(\frac{\varepsilon}{1{\rm\,TeV}}\right)\left(\frac{R_{\rm
sh}}{1{\rm\,au}}\right)^{-1}\,.$ (S31)

Here $0.2\sigma_{\rm T}$ is the maximum cross-section for the pair creation
process expressed in the Thomson cross-section value ($\sigma_{\rm T}$). This
maximum is achieved for $\varepsilon h\nu_{\rm target}\approx
3.5m_{e}^{2}c^{4}$, where $h\nu_{\rm target}$ is the target photon energy. The
attenuation of TeV gamma-rays occurs due to the interaction with optical – IR
photons, and X-rays provide the dominant target for GeV photons.

Equation (S31) shows that the $\gamma-\gamma$ attenuation is the most relevant
for VHE photons. For this gamma-ray energy range, optical and IR photons
provides the dominant target, and below we compute the corresponding optical
target photons. During the first several days after the eruption the typical
optical – IR luminosity remains at the level of $(1-5)\times 10^{37}\rm
erg\,s^{-1}$ (see Fig. S6 and (?)), while the typical X-ray luminosity is
fainter. For example, X-ray fluxes at the level of
$1\rm\,keV\,cm^{-2}\,s^{-1}$ were reported at days $3-6$ after the eruption in
2006 (?). For the adopted source distance of $1.4\rm\,kpc$ this corresponds to
the X-ray luminosity of $\sim 3\times 10^{35}\rm erg\,s^{-1}$. This implies
that the $\gamma-\gamma$ attenuation might be important for GeV photons only
during the first hours after the eruption, when the source size is small,
$R_{\rm sh}\ll 1{\rm\,au}$.

To estimate its impact we compute also the attenuation factor, ${\rm
exp}\left[-\tau\right]$, averaged over the forward shock. The resulting
factors, again assuming a distance of 1.4 kpc are shown in Figure S10 for five
different epochs: $t=1$, $2$, $3$, $4$, and $5$ days after the explosion. The
absorption is not sufficient to account for the observed difference between
the HE and VHE bands, and therefore must reflect a feature in the particle
population responsible for the gamma rays. The impact of the $\gamma-\gamma$
attenuation is shown in Fig. S10B where it can be seen that the spectrum
transformation by the absorption is minor. We note that the spectral energy
distributions in Figs. S8 and S9 account for the attenuation factor obtained
for the source distance of $1.4{\rm\,kpc}$ shown in Fig. S10A.

The idealised single zone model presented here does not accurately account for
the inhomogeneous structure of the post-shock medium, in which strong magnetic
field and gas density gradients naturally develop. The purpose of this model
is therefore not to precisely fit the data, but rather to reproduce essential
features of the system’s evolution, characterised by a minimally sufficient
number of free variables.

Figure S10: Effect of gamma-gamma attenuation. A: Attenuation factor due to the gamma-gamma absorption on the thermal photons produced by the explosion assuming a distance of $1.4$ kpc. B: The impact of the gamma-gamma absorption on the emission spectrum from night 1 for the source assumed distance of $1.4{\rm\,kpc}$. The data points are the same as in Fig. S9. Parameter | Symbol, unit | p-p model | IC model
---|---|---|---
Acceleration slope electrons | $\alpha_{e}$ | – | $2.2$
Acceleration slope protons | $\alpha_{p}$ | $2.2$ | –
Cutoff exponent electrons | $\beta_{e}$ | – | $0.5$
Cutoff exponent protons | $\beta_{p}$ | $0.5$ | –
Fraction of energy in electrons | $\kappa_{e}$ | $0$ | $3\%$
Fraction of energy in protons | $\kappa_{p}$ | $50\%$ | $0$
Acceleration efficiency of electrons | $\eta_{e}$ | – | $10\pi$
Acceleration efficiency of protons | $\eta_{p}$ | $30\pi$ | –
Escape efficiency | $\xi_{\rm esc}$ | $10^{-2}$ | –
Electron low energy cutoff | $E_{\mathrm{min}}$ | – | $10^{2}m_{e}c^{2}$
Proton low energy cutoff | $E_{\mathrm{min}}$ | $2.m_{p}c^{2}$ | –
RG surface magnetic field | $B_{*},\rm\,G$ | $1.$
RG radius, au | $R_{*},\rm\,au$ | $0.35$
RG mass-loss rate | $\dot{M}/v_{\rm w},\rm\,{\rm g~{}cm}^{-1}$ | $6.3\times 10^{11}$
WD orbit radius | $r_{\rm orb},\rm\,au$ | $1.48$
Distance from Earth | $\rm kpc$ | $1.4$
Ejecta initial speed | $v_{\mathrm{ej},0},~{}{\rm km\,s}^{-1}$ | $3000$
Ejecta mass | $m_{\mathrm{ej}}$ | $10^{-7}M_{\odot}$
Table S4: Numerical model parameters. Assumed parameters of the model used in
equations (S3) through (S30) above, for the two cases of proton-proton (p-p)
emission (see Figure S9) and a leptonic Inverse Compton (IC) model (see Figure
S8) respectively.