Gamma-ray light curves and spectra of classical novae

ABSTRACT

The nucleosynthesis in classical novae, in particular that of radioactive isotopes, is directly measurable by its γ-ray signature. Despite decades of observations, MeV γ-rays from novae have never been detected – neither individually at the time of the explosion, nor as a result of radioactive decay, nor the diffuse Galactic emission from the nova population. Thanks to recent developments in modelling of instrumental background for MeV telescopes such as INTEGRAL/SPI and Fermi/GBM, the prospects to finally detect these elusive transients are greatly enhanced. This demands for updated and refined models of γ-ray spectra and light curves of classical novae. In this work, we develop numerical models of nova explosions using sub- and near-Chandrasekhar CO white dwarfs as the progenitor. We study the parameter dependence of the explosions, their thermodynamics and energetics, as well as their chemical abundance patterns. We use a Monte Carlo radiative transfer code to compute γ-ray light curves and spectra, with a focus on the early time evolution. We compare our results to previous studies and find that the expected 511-keV-line flash at the time of the explosion is heavily suppressed, showing a maximum flux of only |$10^{-9}\, \mathrm{ph\, cm^{-2}\, s^{-1}}$| and thus making it at least one million times fainter than estimated before. This finding would render it impossible for current MeV instruments to detect novae within the first day after the outburst. Nevertheless, our time-resolved spectra can be used for retrospective analyses of archival data, thereby improving the sensitivity of the instruments.

1 INTRODUCTION

1.1 Nova physics

Classical novae are transient events originating from the thermonuclear runaway (TNR) on the surface of white dwarfs (WDs), resulting in a mass outburst (see e.g. reviews from Starrfield, Iliadis & Hix 2016; Chomiuk, Metzger & Shen 2020). Multiple physical models exist to explain the diversity of the outburst phenomenon, including mass transfer from a companion star, symbiotic star, and dwarf novae (Webbink et al. 1987). The mass transfer model has been studied widely: It includes a binary system in which a compact object, i.e. a WD or neutron star (NS), accretes H- or He-rich matter from its companion main-sequence star through a filled Roche Lobe (Paczyński 1965; Starrfield & Sparks 1987). The later TNR is developed due to the thermally unstable H-burning shell (Rose 1968; Starrfield 1971). After the mass outburst, the luminosity decreases until it re-approaches the initial luminosity so that the mass accretion can resume (MacDonald, Fujimoto & Truran 1985).

The dynamics of novae are non-trivial because it might involve the accretion disc with an aspherical outburst. The ignition and the beginning of the mass outburst begin depend on the mass transfer rate from the companion star, the progenitor mass, and the composition and mixing (Kovetz & Prialnik 1985; Truran & Livio 1986; Townsley & Bildsten 2004; Shen & Bildsten 2009; Denissenkov et al. 2013; Ginzburg & Quataert 2021). The mixing is essential for explaining the observed metal abundances in the ejecta (Prialnik 1986). The rapidly rotating accretion disc creates a rotational instability that triggers element mixing of the H-rich matter in the accretion disc and the CO- (or ONe-) rich matter of the WD (Prialnik & Kovetz 1984; Fujimoto 1988), by turbulent mixing (Fujimoto 1993), convective overshooting (Glasner & Livne 1995), and Kelvin-Helmholtz instability (José, Shore & Casanova 2020). Angular momentum is also transported in the process (Kutter & Sparks 1989). The mixing process driven by the wave breaking is important to the ignition because otherwise, the CO-rich matter cannot diffuse efficiently into the H-envelope, making the TNR weak (Alexakis et al. 2004).

Although novae show an occurrence rate about 100–1000 times higher than supernovae in a galaxy, their contribution of global metal enrichment is small compared to supernovae because they eject a much lower mass (⁠|$\lesssim 10^{-7}\!-\!10^{-5}\, \mathrm{M}_{\odot }$|⁠) during the outburst (Hernanz et al. 1996; José et al. 2020). Nevertheless, novae are ideal laboratories for the production of low-mass elements that eventually mix with the interstellar medium (Gehrz et al. 1998). For example, they produce an enhanced abundance in β-decay isotopes of the CNO cycle, including ¹³C, ¹⁵N, and ¹⁷O (Starrfield et al. 1972). Most reactions rely on the hot-CNO cycle and proton capture. The TNR can also trigger the formation of intermediate (²²Na) and longer (²⁶Al) lifetime radioactive isotopes (Weiss & Truran 1990). The radioactive buildup from the population of novae in the Milky Way is expected to contribute to the diffuse γ-ray emission along the Galactic plane (Diehl et al. 2021). It is also a rich source of ⁷Li (Starrfield et al. 1978, 2020), suggested by the decay signal of ⁷Be (Harris, Leising & Share 1991; Tajitsu et al. 2015). Directly measuring the yields in nova explosions can help to constrain the uncertainties of nuclear reaction rates (José, Coc & Hernanz 2001) and the mystery of the cosmic origin of lithium and fluorine (Spitoni et al. 2018; Ryde et al. 2020).

1.2 Gamma-ray observations

Clayton & Hoyle (1974) pointed out that suited γ-ray telescopes could observe the essential information to prove the intricacies of the nova phenomenon. In particular, the authors suggested the 2312 and 1275 keV line from ¹⁴O and ²²Na, respectively, and the 511 keV line from various β⁺-unstable isotopes including ¹³N, ¹⁴O, and ¹⁵O. Follow-up studies (e.g. Clayton 1981; Leising & Clayton 1987; Hernanz & José 2006; Hernanz 2014) then also suggested that the 478 keV line from ⁷Be should be visible up to distances of ∼2 kpc. To this day, no individual nova nor the diffuse radioactive glow of the nova population have been detected in soft γ-rays.

The three major objectives, the 1275 keV line from ²²Na, the 478 keV line from ⁷Be, and the 511 keV line flash from various short-lived isotopes, had been studied in several cases: Using CGRO/COMPTEL, Iyudin et al. (1999) found a limiting ejecta mass of ²²Na of |$\lt 2.1 \times 10^{-8}\, \mathrm{M_{\odot }}$| from Nova Cygni 1992. Trying to observe the Galactic ridge, Jean et al. (2001) constrained the ²²Na ejecta mass to |$\lt 3 \times 10^{-7}\, \mathrm{M_{\odot }}$|⁠. By applying a Bayesian hierarchical model to INTEGRAL/SPI data, Siegert et al. (2021) could limit the ²²Na ejecta mass in the Milky Way to |$\lt 2 \times 10^{-7}\, \mathrm{M_{\odot }}$| per nova. The best target for detecting the 478 keV line within the last 20 yr was V5668 Sgr, because of its proximity of 1–2 kpc. Molaro et al. (2016) detected ⁷Be II lines in the UV which suggest a ⁷Be ejecta mass of |$\sim 0.7 \times 10^{-8}\, \mathrm{M_{\odot }}$|⁠. Unfortunately, also, this object was too far away, and Siegert et al. (2018) determined an upper limit of |$\sim 1.2 \times 10^{-8}\, \mathrm{M_{\odot }}$|⁠. A ⁷Be ejecta mass of the order |$10^{-8}\, \mathrm{M_{\odot }}$| is about two orders of magnitude larger than what is theoretically expected, suggesting either a gap in the theoretical understanding or an observational bias.

Finally, the UVOIR observations of novae typically happen at the bolometric maximum luminosity that is expected to be days to weeks after the explosion (Gomez-Gomar et al. 1998). Short-lived isotopes, such as ¹³N (⁠|$\tau _{13} = 10\, \mathrm{min}$|⁠) and ¹⁸F (⁠|$\tau _{18} = 110\, \mathrm{min}$|⁠), are produced during the explosive nucleosynthesis and decay by positron emission. This suggests a strong 511 keV line and low-energy continuum within the first few hours after the explosion – but days to weeks before the nova had been discovered. Therefore, only a retrospective analysis of archival data could find such a signal. The latest of such searches has been conducted by Skinner et al. (2008) using Swift/BAT, but resulted in no detection. Other works that tried to discover the annihilation flash in archival data include Harris et al. (1999, 2000) using TGRS/WIND, Hernanz et al. (2000) using CGRO/BATSE, and Smith (2004) using RHESSI.

1.3 Motivation and structure

Searches in light curves of individual energy bands or lines limit the sensitivity of the instruments to only these particular search windows. In addition, using extracted light curves, i.e. assuming a position in the sky but being agnostic about the temporal and spectral behaviour, will bias the result and consequently the search for transients. The most efficient use of the capabilities of γ-ray instruments is forward modelling of all available information in a complete model of the source to be analysed. This means instead of searching for high count rates in the 511 keV line, for example, a complete spectro-temporal model is convolved with the instrument response, and fitted together with an appropriate description of the instrumental background. One example to implement such models is the Multi-Mission Maximum Likelihood (3ML) framework (Vianello et al. 2015), designed to perform not only multiwavelength analyses but also multi-instrument inferences.

The other notorious problem in soft γ-ray astrophysics is the instrumental background from cosmic ray bombardment of the spacecraft, and also the astrophysical background in all-sky monitoring instruments. Only if these backgrounds are sufficiently well-understood and modelled, weak signals, such as expected from nova explosions, can be disentangled from the raw photon count data. In recent years, significant progress has been made in MeV astrophysics to understand the often erratic backgrounds, and utilize the instruments such as INTEGRARL/SPI (cf. Siegert et al. 2019, 2021) and Fermi/GBM in ways they have not been designed for. In particular, Biltzinger et al. (2020) developed a physical background model for the the all-sky monitor GBM onboard Fermi. This also allows to study diffuse emission and sets the basis for the analysis of hour-scale transients such as nova events.

What is missing are parametrized and public time-resolved spectra of novae that can be used with these recently developed tools to maximize the scientific return. Therefore, in this article, we will first describe our numerical schemes in the stellar evolution and Monte Carlo radiative transfer code of γ-ray lines in Section 2. In Section 3, we present the results of our grid of nova models and discuss their dependence on model parameters. We discuss the corresponding nucleosynthesis with a focus on the radioactive isotopes in Section 4. Section 5 presents our numerical results for the γ-ray light curves and spectra derived from our models. In Section 6, we explore how the generation of the 511 keV line, ortho-Positronium (ortho-Ps) continuum, and Compton scattered photons from novae can be utilized in a retrospective analysis of archival γ-ray data, and how our findings impact previous analyses. Finally, in Section 7, we compare our nucleosynthesis and radiative transfer results with representative works in the literature, discuss numerical issues and caveats of this work, and present our conclusion.

2 METHODS

2.1 Stellar evolution and nova simulation

We use the stellar evolution code mesa version 8118 (Paxton et al. 2011, 2013, 2015, 2018, 2019). The code solves the structure, nuclear reactions, and radiative transfer inside a star with spherical symmetry. Similar to our previous works (Leung et al. 2020; Leung, Fuller & Nomoto 2021b; Leung, Wu & Fuller 2021a), we carry out the stellar evolution and the radiative transfer calculation into two separate steps. The calculation is developed based on the make_co_wd and nova test cases taken from the code. We modify the configuration files and add subroutines for this work.

We first use the package make_co_wd to construct a CO WD with a mass from 0.8 to 1 M_⊙. The code first constructs a 3–7 M_⊙ star until it develops the CO core, and then increases the optical depth of the star to enhance mass loss of the H-envelope. Then we pass the model to the nova package and make the WD accrete H-rich matter. We assume the matter originates from a single degenerate companion star. The matter is mixed with ¹²C and ¹⁶O as impurities. The fraction of CO-rich matter is taken as a model parameter f_CO from 10 per cent to 50 per cent. This mimics different levels of dredge-up of the CO-rich core into the H-rich envelope. Once the H-burning runaway takes place, the code allows excited matter in super-Eddington luminosity to be ejected as wind. We keep track of the synthesized radioactive isotopes of interest (⁷Be, ¹³N, ¹⁴O, ¹⁵O, ¹⁸F, ²²Na, and ²⁶Al) and the ejecta kinematics. In Table 1, we list the isotopes included for the nucleosynthesis calculations. We limit our nuclear network here because we want to focus on the radioactive isotope production described. The chosen network is large enough to include potential nuclear reactions for their formation while maintaining a feasible computational time. In general, the nuclear reaction may proceed to low mass iron-group elements such as Ca (José & Hernanz 1998; José et al. 2001).

2.2 Monte carlo radiative transfer

We use the recorded data and pass it to our Monte Carlo radiative transfer code. The code solves the propagation, scattering, and interaction of photon packets produced by the decay of radioactive nuclei. At each time interval, we estimate the optical depth of the matter for γ-rays using a grey opacity κ_γ = 0.06Y_e (Swartz, Sutherland & Harkness 1995). We choose the mass shells with an optical depth |$\tau = \int _r^R \rho (r^{\prime }) \kappa _{\gamma } \mathrm{ d}r^{\prime } = 5$|⁠. This corresponds to about ∼99 per cent reduction of the luminosity from the innermost layer considered when it arrives the surface.

To generate the photon packet, we calculate the radioactive power per mass q_decay in each mass shell by including the decay of all related isotopes. We locate the mass shells Δm that are ejected and integrate these shells by L_total = ∑q_decayΔm to find the instantaneous γ-ray emission. They correspond to a collection of photons undergoing the same process during their propagation inside the star (Ambwani & Sutherland 1988). Instead of splitting into multiple packets when scattering occurs, they can lose energy. Each packet therefore corresponds to the ‘instantaneous photons emitted per unit time’ instead of a real photon. Thus, a packet carries a luminosity L_packet that corresponds to L_total/N_packet. In general, we find that the time delay from photon emission to escape is short (∼1 ms or lower) so that the static approximation is appropriate.

2.2.1 Photon generation

The code assumes that the γ-rays originate solely from the decay of radioactive isotopes and we calculate the spectra from given time snapshots of the nova model during its mass ejection. We put particular emphasis on the first day of the mass ejection as we are interested in the decays of ¹³N and ¹⁸F, which have short half-life times compared with the expansion time-scale of the nova (∼10 d for our models). These isotopes are distinctive from the others because their main decay channel is β⁺-decay, which emits a positron or in addition to (in the case of ¹⁵O) an energetic photon. In Table 2, we tabulate the principle parameters for the isotopes of interest.

If an isotope decays and subsequently emits a photon, a photon packet is created by assigning it the corresponding γ-ray line energy and an arbitrary direction. The γ-ray energies are taken from the co-moving frame with the moving nuclei. In our calculation, the thermal fluctuation and collective motion of the nuclei is small enough that the Doppler effect is negligible.

When an isotope undergoes β⁺-decay, the positron quickly loses its energy through Coulomb interaction and captures an electron to form Ps. Ps has two spin states: para-Ps and ortho-Ps. The former emits two photons and the latter three due to charge and spin conservation (see Ore & Powell 1949; Berko & Pendleton 1980, for fundamental features of Ps). Quantum statistics limits their ratio to a maximum of para:ortho=1:4.5. The exact ratio depends strongly on the matter density and temperature (e.g. Leising & Clayton 1987). In our case, where the matter is opaque and dense, we expect the ratio to approach the quantum limit. When two-photon emission occurs, we assign two photon packets of energy 511 keV. The first one has an arbitrary direction, with the direction of the second packet chosen by momentum conservation. When three-photon emission occurs, the individual photon energy is set by Monte Carlo process that reproduce the observed probability distribution, and the sum of the three photons conserves the total energy 1022 keV and the initial momentum. Then we assign an arbitrary direction for one photon packet and obtain the directions for the other two photon packets by conservation of momentum. The time delay from the formation of Ps to its decay is 10⁻⁹–10⁻⁶ s (Czarnecki & Karshenboim 1999), which is much shorter than its escape time.

2.2.2 Photon interactions

We consider the following three types of interaction processes that change the energy of a photon packet (Pozdnyakov, Sobol & Syunyaev 1983):

• Compton scattering: the photon transfers energy to an electron and loses energy. The energy before E and after E′ the interaction is related by:

  $$\begin{eqnarray} E^{\prime } = \frac{E}{1 + (E/m_\mathrm{ e} c^2) \cos \theta }, \end{eqnarray}$$

  (1)

  where m_e is the electron mass and θ is the scattered angle in the centre-of-mass frame.
• Photopair production: the photon packet is assumed to lose all its energy to an electron and the electron later emits an e⁻-e⁺-pair similar to the two-photon emission described above. The cross-section of this process is calculated according to the photon energy:

  $$\begin{eqnarray} \sigma &=& A (E - 1.022) Z^2 \times 10^{-27} {\rm cm^2}, ~ 1.022 \lt E \lt 1.5 {\rm MeV}; \\ \sigma &=& [B_1 + B_2 (E - 1.5)] Z^2 \times 10^{-27} {\rm cm^2}, ~ E \ge 1.5 {\rm MeV}, \end{eqnarray}$$

  (2)

  with A, B₁, and B₂ being 1.0063, 0.481, and 0.301, respectively. E is the photon packet energy in units of MeV.
• Photoelectric absorption: the photon is assumed to be absorbed by a nucleus. The cross-section takes the equation:

  $$\begin{eqnarray} \sigma = C (E / {\rm 100~keV})^{-3.2} \times 10^{-24} {\rm cm^2}, \end{eqnarray}$$

  (3)

  with C = 0.0448 for elements with Z ∼ 7. Although the exact value of C depends on the atomic number, σ ∼ E⁻³ holds true for a wide range of elements, which means photoelectric absorption dominates the destruction of photon packets as the photon energy reaches below ∼100 keV.

We assume that the photon has at least one interaction when it travels across δτ = 1. To determine which process takes place, in each step, we assign a random number Δτ ∈ (0, 1) that corresponds to the distance travelled by the photon packet before the next interaction. The mean free path Δr_i = Δτ/(ρκ_i) for κ_i is taken from the three processes described above. By finding the minimum Δr_i, we assign the corresponding probability for each process to take place. Then we generate another random number to decide which process is chosen. In the case where the travelling distance crosses the mass shell (defined by the stellar evolution model), we also update the local thermodynamical properties experienced by the photon packet.

3 EVOLUTION OF NOVAE

We generate a grid of nova models by varying the WD mass M (0.8–1.2 M_⊙), CO fraction f_CO (10–50 per cent), mass accretion rate |$\dot{M}$|⁠, and convection parameters. In Fig. 1, we show the luminosity and radius as a function of time for |$M = 1.0\, \mathrm{M}_{\odot }$|⁠, f_CO = 0.5 and |$\dot{M} = 10^{-9}\, \mathrm{M}_{\odot }$| yr⁻¹. We plot two pulses in the evolution of a recurrent nova. During the nova outburst, the luminosity can increase by five orders of magnitude and reach |$\sim 10^{4}\, \mathrm{L}_{\odot }$|⁠. The WD radius also increases by one order of magnitude during the nova outburst. The two pulses are separated by about 7000 yr. The exact timing depends on the mass accretion rate, the mass of the WD and mixing of accreted matter. For recurrent novae, the mass accumulated before the TNR occurs is roughly similar. In the right-hand panel of the same plot, we magnify the first nova outburst. Zooming into the first outburst, it becomes clear that the whole outburst lasts about 0.2 yr until most thermally excited matter is ejected. It takes a much longer time (∼10–100 yr) for the nova to relax and to restore its original luminosity. The decline of the luminosity satisfies L = L₀t⁻¹, suggesting the WD cools by blackbody radiation. The decrease of the radius happens because the super-Eddington wind gradually removes the hot matter in the outer parts of the nova. The mass ejection is efficient in removing the deposited matter until it reaches its original mass before mass accretion.

In Fig. 2, we show the stellar profiles for the same nova model as above, focusing on the thermodynamics and kinematics. The snapshot is taken before the wind mass loss removes the thermally excited matter.

The nova consists of two layers, the compact CO core with a radius |$\sim 10^{-2}\, R_{\odot }$| and the sparse accreted H-rich matter that extends to |$\sim 10^{-1}\, R_{\odot }$|⁠. The temperature bump near |$r = 10^{-2}\, R_{\odot }$| also shows the location where explosive H-burning has happened. This is consistent with the luminosity profile, also being increased at the interface and gradually decreasing as the nuclear reactions occur in a mass shell with a small mass |$\sim 10^{-5}\, \mathrm{M}_{\odot }$|⁠. From the velocity profile, only the outer part of the H-envelope is expanding, but with a very low velocity ∼1–|$3\, \mathrm{km\, s^{-1}}$|⁠. The expansion does not lead to a dynamical mass outburst because the typical escape velocity is |$\sim 10^3\, \mathrm{km\, s^{-1}}$|⁠. Instead, the nova loses its matter by radiative winds. During TNR of the H-envelope, the heat creates strong convection, which allows the produced radioactive isotopes to be efficiently transported and some γ-rays to escape easily. Only those near the surface are important for the early γ-ray signal because of the very short lifetime of ¹³N and ¹⁸F.

We surveyed the nova models by spanning the progenitor WD mass and the mixing efficiency. In Fig. 3, we show the recurrence times for novae as depending on the mixing, accretion rate, and WD composition. The recurrence time is defined by the time lapse between the first and second outburst. The models show a tight exponential relation for models with 0.8–1.2 M_⊙ for CO WDs. Beyond 1.2 M_⊙, the degeneracy leads to models deviating from this trend. A similar relation is observed for ONe WDs between 1.1–1.2 M_⊙.

We determine a smooth relation of the recurrence times by fitting the function t_recur = Aexp (− BM), where M is the WD mass, and A and B are fitted constants as listed in Table 3. The fitting suggests that B ∼ −5.5 for a wide range of f_CO. The choice of the mass accretion rate does not change the trend of the models. An exception occurs for the CO WD model with |$\dot{M} = 10^{-10}\, \mathrm{M}_{\odot }$| yr⁻¹ and |$M = 1.3\, \mathrm{M}_{\odot }$|⁠, which has an exceptionally high recurrence time. The ONe WDs also show a similar trend which can be well-described by the same scaling.

To estimate the expansion velocity, we show the velocity and density profiles as a function of time. For this part, we repeat the model but suppressed all super-Eddington mass loss so that we can follow how the ejecta develop towards their asymptotic velocity during expansion. As an example, we show the nova models of a CO WD with |$\dot{M} = 10^{-9}\, \mathrm{M}_{\odot }$| yr⁻¹, f_CO = 0.5, and a mass of |$M = 0.8\, \mathrm{M}_{\odot }$| and |$M = 1.3\, \mathrm{M}_{\odot }$|⁠, respectively, in Figs 4 and 5.

In the first model it takes about 8–10 d for the motion to reach the surface and develop an expansion flow similar to homologous expansion. At early times (≲ 8 d), the ejecta have a low velocity of |$\sim 1\, \mathrm{km\, s^{-1}}$| and are decelerating with time, showing that they remain bounded by gravity. Only around day 8 when the momentum reaches the low-density tail of the ejecta, they gradually develop the rapid expansion with an asymptotic velocity 200–|$300\, \mathrm{km\, s^{-1}}$|⁠. Once such a velocity is reached, the ejecta no longer accelerate. Higher mass models have a stronger deceleration at early times, but the density structure for the low-density tail is similar. Both profiles suggest that the mass ejection process is similar regardless of the progenitor. The slow velocity at early times also suggests that the ejecta with radioactive isotopes remain compact for a time longer than the half lives of, for example, ¹³N and ¹⁸F. Such an opaque environment largely suppresses the propagation of γ-rays.

In the temperature and luminosity profiles, we observe that during the TNR of H, the peak temperature can reach as high as 10^8.2 K. This provides the necessary thermal pressure for the expansion, where we see the long tail of matter with its temperature falling to about ∼10⁵ K around day 10. The high temperature allows most elements in the ejecta to remain ionized. The peak of nuclear reactions can raise the luminosity to about |$10^8\, \mathrm{L}_{\odot }$|⁠. After that the luminosity quickly falls to a constant level of about |$10^{4.5}\, \mathrm{L}_{\odot }$| when the surface matter expands.

To understand how the matter ejection depends on the WD progenitor, we examine the velocity evolution of the near-Chandrasekhar mass CO WD model with |$M = 1.3\, \mathrm{M}_{\odot }$|⁠. The evolution is almost identical to the lower mass models despite the large difference in mass and initial radius. The asymptotic velocity is slightly higher, about |$300\, \mathrm{km\, s^{-1}}$|⁠, but is marginally below the escape velocity. The ejecta maintain a high velocity on day 10 when it reaches a radius of |$10\, R_{\odot }$|⁠, suggesting that they will become unbounded at later times. Again, the development of the velocity flow takes about 9 d before the velocity reaches its asymptotic value. Thus, the short lifetime radioactive isotopes, such as ¹³N and ¹⁸F – the drivers for the annihilation flash – have decayed before the matter becomes transparent to γ-rays.

To further explore the trend of the mass deposition and ejection of the nova events, we plot the total ejected mass in Fig. 6. The models show a clear falling trend with increasing progenitor WD mass. The typical ejected mass decreases from ∼2–|$3 \times 10^{-5}\, \mathrm{M}_{\odot }$| for WDs with masses of 0.8 M_⊙ to ∼5–|$10 \times 10^{-7}\, \mathrm{M}_{\odot }$| for WDs near the Chandrasekhar mass. The models show a small variation among different mixing efficiencies. A similar falling trend is observed for ONe WDs with a 10-fold larger ejecta mass compared to CO WDs with the same mass.

We show the global maximum temperature of the models during the outburst in Fig. 7. Contrasted to the ejected mass, the maximum temperature experienced in the star increases with the progenitor mass approximately exponentially. The variation among models for the same mass is small, within ∼0.05 in logarithmic scale. The maximum temperature always occurs at the bottom of the H-burning zones. A higher maximum temperature is expected because the corresponding density for the H-burning shells is higher for a higher mass WD. As a result, the nuclear reactions can proceed more readily once the temperature is sufficiently high. ONe WDs also show an increasing trend with a similar slope, however, at higher temperatures around 220 MK. Such a high temperature may allow more radioactive isotopes to be synthesized.

4 NUCLEOSYNTHESIS IN NOVAE

The stellar evolution models with the extended nuclear reaction network allows us to compute explicitly the detailed chemical composition, and especially those of radioactive isotopes, in the ejected matter. In each model, we capture the chemical abundance before the super-Eddington wind mass loss begins. This ensures that the all synthesized metal is properly captured.

4.1 Productions of radioactive isotopes

In Fig. 8, we show the nucleosynthesis results of our nova models. The models are classified based on the composition (CO or ONe) and the mass accretion rate from its companion stars.

It is evident that there is no clear trend for the lighter isotopes, including ⁷Be, ¹³N, and ¹⁸F. The scatter of the mass can range from one to three orders of magnitude and without a consistent mass dependence. ⁷Be has a narrower range from 10⁻¹²–|$10^{-11}\, \mathrm{M}_{\odot }$|⁠, while ¹³N (¹⁸F) can range between 10⁻⁹ and |$10^{-6}\, \mathrm{M}_{\odot }$| (10⁻⁸–|$10^{-5}\, \mathrm{M}_{\odot }$|⁠). There is also no specific trend with the WD composition. On the other hand, there is a clear decreasing trend for ²²Na as a function of WD mass. ONe WDs have a similar trend but with a typical mass about 10–100 times higher than CO WDs (10⁻¹²–|$10^{-10}\, \mathrm{M}_{\odot }$|⁠).

4.2 Production of stable isotopes

We also examine the chemical abundance pattern of our nova models with a focus on the lower mass elements C, N, and O. Most of them are produced in massive stars, but with a metallicity comparable to solar values. Novae are often associated with super-solar production. Here, we examine how the production of these elements compare to the observational data. We use measurements reported in Schmidt et al. (2018) who derived the chemical abundances from the planetary nebula K4-47. The low ¹²C/¹³C, ¹⁴N/¹⁵N, and ¹⁶O/¹⁷O ratios highly deviate from the solar value, implying that standard massive star explosions are unlikely to be the origin. Their ratio represents how the proton capture competes with the CNO cycle. Another tight measurement from Nova Oph 2017 is given in Joshi et al. (2017), which also shows a very low ¹²C/¹³C ratio of 1.6 ± 0.3. In Fig. 9, we show the three ratios of our models as a function of WD progenitor mass in comparison to these measurements.

There are subtle trends in the models albeit the large fluctuations among models. For ¹²C/¹³C, the mass ratio increases with the progenitor mass. ONe WD models have a higher mass ratio than the CO WD models for high WD masses. For ¹⁴N/¹⁵N, the mass ratio clearly decreases with the progenitor mass. The ONe WDs show distinctively lower mass ratios. There is no clear trend in the ¹⁶O/¹⁷O mass ratio against the progenitor mass and ONe WDs again have the lowest value among all models. The clustering of the data points suggests two trends to classify the CO and ONe WDs. ONe WDs tend to have lower ¹⁴N/¹⁵N (<100) and lower ¹⁶O/¹⁷O (<50) ratios. These thresholds could serve as an indicator of the WD composition.

For the observational data from Schmidt et al. (2018), each observed isotope ratio suggests different mass ranges or WD types from our models. Given the ¹²C/¹³C-values from our models, progenitor masses of 1.0–1.1 M_⊙ (or |$\sim 0.8\, \mathrm{M}_{\odot }$|⁠) are suggested. From the ratio ¹⁴N/¹⁵N as well as ¹⁶O/¹⁷O, WD masses with 1.1–1.3 M_⊙ are favoured. However, the overlap region around |$1.1\, \mathrm{M}_{\odot }$| would either suggest a low accretion rate on a CO WD (¹²C/¹³C), high accretion rate on a CO WD (¹⁴N/¹⁵N), or an ONe WD (¹⁶O/¹⁷O). The displacement of the models with the data suggests that further input physics is required in matching all constraints at the same time, and it will be an interesting future project for an extensive exploration.

For the observational data from Joshi et al. (2017), our models show that CO WDs with a mass 1.0–1.1 M_⊙ or a |$1.1\, \mathrm{M}_{\odot }$| ONe WDs can approach the narrow ¹²C/¹³C mass ratio observed in this nova. The results are consistent with the spectrographic data in the UV band that also suggested that the observed WD has a mass of 1.0–1.1 M_⊙ (Mikołajewska & Shara 2017). In their work, the ONe WD is excluded. Our CO WD models with a low f_CO can fit the mass and the mass ratio simultaneously.

5 GAMMA-RAY RADIATIVE TRANSFER

After obtaining the kinematics of the ejecta from Section 3 and the chemical abundance from Section 4, we pass the information to our Monte Carlo radiative transfer code and compute the corresponding γ-ray spectrum. From the spectrum, we derive the γ-ray luminosity by aggregating the escaped photon packets. To construct the spectra, we use logarithmic bins for photon energies below 508 keV and above 514 keV, up to 3 MeV.

5.1 Nova Gamma-ray spectra

In Fig. 10, we show the γ-ray spectra snapshots for different models on day 0.25, 0.50, 1.00, and 2.00. In all spectra, we observe that the low energy part 30–200 keV is dominated by the back-scattering of higher energy photons. Up to day 0.50, the 511 keV line from the β⁺-ecays of ¹³N and ¹⁸F is dominating the spectra. Starting on day 1.00, the 478 keV line from ⁷Be also emerges, and becomes important for later time as well due to its longer lifetime. After day 1.00, the 511 keV line almost vanishes due to the exhaustion of ¹³N and ¹⁸F. The spectral shape is almost unchanged for later times except that the absolute magnitude is decreasing with time. There is an inversion in the CO WD models where the absolute magnitude of the spectra is greater on day 1.00 than 2.00. This is because ⁷Be and ²²Na, which support the later γ-ray luminosity, locate at the inner ejecta so that it takes a longer time for the photosphere to recess and expose these isotopes.

We find that the γ-ray luminosity is larger for ONe WD models compared to CO WD models. ONe WD models have a higher accretion and outburst mass because ONe-rich matter requires higher masses to incinerate before the TNR to take place. As a result, it generates more radioactive isotopes and hence stronger γ-ray emission.

A higher progenitor mass is beneficial for the emission of γ-rays. This is because the maximum temperature reached during TNR is higher when the mass is higher. The higher temperature allows more rapid nuclear reactions as well as higher thermal pressure, which helps to expel the matter. Even though the progenitor WD is more compact for a higher mass WD, the overall burning is still stronger.

Finally, the γ-ray signal is sensitive to the mixing ratio f_CO of the CO (and without loss of generality ONe) matter. Because the presence of CO matter facilitates the outburst process not only by its interval, but also its tendency to burn. A high f_CO allows the outburst of the envelope with a smaller ejecta mass and a higher velocity that facilitates the transport of γ-rays.

5.2 Nova Gamma-ray light curves

By collecting the corresponding γ-ray luminosity in the energy band from 10 keV to 3 MeV, we show the light curve of representative models in Fig. 11. Unlike for type Ia supernovae, where the luminosity is very low at early times, the γ-ray luminosity for novae is the highest at the beginning: Because of the ‘exposed’ ¹³N and ¹⁸F on the surface, the direct formation and decay of Ps leads the γ-rays to readily escape from the star. The luminosity reaches as high as |$\sim 10^{28}\, \mathrm{erg\, s^{-1}}$|⁠. However, the luminosity quickly drops by five to six orders of magnitude soon after. The recession of the photosphere is slow so that the escape of photon packets from the inner layers is very inefficient. After day 1.00 after the outburst, the γ-ray luminosity gradually increases and reaches an asymptotic value of |$\sim 10^{25}\, \mathrm{erg\, s^{-1}}$|⁠. We find that the early emission is dominated by the 511 keV line as has been found in previous studies.

We examine the dependence of the γ-ray light curve on the WD progenitor. A higher mass WD generates a higher γ-ray luminosity but it converges quickly for models with |$M \gt 1.1\, \mathrm{M}_{\odot }$|⁠. The changes due to different |$\dot{M}$| are small, as the light curve is primarily depending on the TNR strength and the synthesized radioactive isotopes (cf. previous sections). The WD composition has a significant effect on the light-curve shape. CO WDs show a rapidly falling luminosity, which is different from the parabolic shape of light curves appearing in ONe WD models. The parabolic shape implies that the lines related to early times are stronger in the ONe WD models, compared to CO WD models.

To further extract the γ-ray signals, we plot in Fig. 12 the components of our reference models assuming |$M = 1.20\, \mathrm{M}_{\odot }$|⁠, |$\dot{M}=10^{-9}\, \mathrm{M}_{\odot }$| yr⁻¹ and f_CO = 0.5 for the CO WD (left-hand panel) and the ONe WD (right-hand panel). Most of the early-time light-curve power in the CO model comes from the scattering background below 80 keV and from the 511 keV line (here we define 511 ± 6 keV as the line), albeit weaker than the scattering background. On day 1 and beyond, the decay of ⁷Be becomes dominant in the light curve. In contrast, 511 keV line becomes insignificant because most short-life radioactive isotopes ¹³N and ¹⁸F have decayed. The ONe model shows a similar hierarchy in its luminosity components. The difference in the chemical composition makes the early time (<1 d) scattering background and 511 keV line stronger for a longer time.

The 487 keV line from ⁷Be decay surpasses the 511 keV line after day 1. Beyond day 10, as the ejecta become transparent in regions, where they have been accelerated by thermal expansion, the escaped γ-ray luminosity approaches the total luminosity by radioactivity decay. For example, the reference model takes about 30 d until the escaped γ-ray luminosity is 99 per cent of the total luminosity.

6 APPLICATIONS TO OBSERVATIONS

In general, our new models can be used to improve γ-ray searches from individual novae as well as to study the cumulative effect of the whole nova population in the Milky Way. First, we note that our flux estimates are at least four to eight orders of magnitude smaller within the first ∼30 d after the explosion compared to previous calculations (e.g. Gomez-Gomar et al. 1998; José & Hernanz 1998; Hernanz 2014, the first paper is called JH98 hereafter). This would make any detection prior to a month after the explosion impossible for a nova happening at any astrophysical distance for current instrumentation. In Table 4, we provide a list of integrated fluxes for characteristic times and energy bands, comparing the brightest model variants of a CO and an ONe WD with |$1.2\, \mathrm{M}_{\odot }$| each, as well as the CO WD simulation setup from JH98. In particular, the supposed ‘511 keV-flash’ h after the outburst that had been suggested to show fluxes on the order of 10⁻³–|$10^{-1}\, \mathrm{ph\, cm^{-2}\, s^{-1}}$| would be reduced to 10⁻¹⁰–|$10^{-9}\, \mathrm{ph\, cm^{-2}\, s^{-1}}$| in our new model calculations.

Such an extremely reduced flux in the entire 0.02–2 MeV energy band has also an impact on the interpretability of diffuse emission measurements of the 478 and 1275 keV lines. For example, Siegert et al. (2021) performed the latest analysis of diffuse emission from novae in the Milky Way. Because the early-time light curves from ⁷Be and ²²Na have typically been assumed to follow an exponential radioactive decay with maximum luminosity around days 4–6 after the explosion (Hernanz 2014), the ejecta masses are wrongly calibrated in such observations. Extrapolating from years after the explosion with pure radioactive decay will lead to higher ejecta mass estimates than what might actually be the case, especially if the peak γ-ray luminosity occurs only several weeks after the explosion. Therefore, with our new models, the cumulative effect of novae in the Milky Way can be estimated with higher accuracy in a population synthesis model, and then fitted via a Bayesian Hierarchical Model similar to one performed by Siegert et al. (2021).

Another application concerns the positron annihilation flash shortly after the explosion. Instead of using only the 511 keV line or individual energy bands, as has been done in previous retrospective studies (e.g. Harris et al. 1999, 2000; Hernanz et al. 2000; Smith 2004; Skinner et al. 2008), MeV data analysis is severely more sensitive to expected signals if the entire time-variable spectrum is used. Soft γ-ray analysis has been recently progressing into the regime of full forward-modelling, i.e. taking into account the full response and position of the source in a single step – without assuming generic spectral shapes such as power laws. Because our models are time resolved, and in addition would include several other parameters such as the mixing fraction, accretion rate, or WD mass, it is possible to perform both, a directed retrospective search for known novae, and a blind search for unknown objects. Because the optical maximum from classical novae occurs 2–10 d after the explosion (Gomez-Gomar et al. 1998), the expected 511 keV flash has already passed. However, there is a chance that individual γ-ray instruments observed the nova by chance, or that an all-sky monitor showed enhanced fluxes but which could not be identified. In a directed retrospective search for MeV emission from past nova outbursts, the expected time frame is revisited in data archives and investigated for coincident emission. With our new time-resolved models and the ability to perform simultaneous fits in time and energy (for example with the 3ML framework, Vianello et al. 2015), a directed retrospective search for γ-ray flashes from novae can be conducted. Ideal instruments currently in space would be Swift/BAT, Fermi/GBM, and ISGRI and SPI on INTEGRAL. With their large fields of views, the chances are high that individual novae happened close enough for them to be detectable. In fact, the number of novae that should have been detected by SPI within the last 19 yr, based on previous flux estimates from JH98 of the 511 keV line alone, is between 1 and 40. Since no nova or 511 keV flash has been reported from SPI, our results of an extremely reduced flux seem to agree with these non-detections.

A general retrospective search for MeV emission of unknown novae in the Milky Way is now also possible: The nova rate in the Milky Way is |$50 \pm 25\, \mathrm{yr^{-1}}$| (Shafter 2017), however, only 10–30  per cent are detected at UVOIR wavelengths each year. Such a search would be very time consuming because in addition to the time and energy domain the position also would need to be determined. Because the highest fluxes are still expected within the first 2–3 h after the explosion, several overlapping time frames of observations will be required to obtain reasonable baselines to search for unexplained emission in addition to known sources. Thanks to developments in Fermi/GBM background modelling (Biltzinger et al. 2020) that also allows to search for longer signals, as well as SPI (Biltzinger et al., in preparation) and Swift/BAT (DeLaunay & Tohuvavohu 2021) data analysis, such a blind search might be worthwhile.

Given the enhanced line and continuum sensitivity of the new NASA SMEX mission COSI (Tomsick et al. 2019), on the order of a few novae should be detectable within its nominal mission time of 2 yr. Our models suggest that the brightest MeV γ-ray emission only occurs several weeks after the nova explosion and not days. Because COSI is observing the full sky within 24 h, an observation strategy is not required here. However, this shifted maximum in the 478 keV line, for example, should be considered for targeting observatories such as INTEGRAL whenever a nearby nova will trigger dedicated observations.

7 DISCUSSION

7.1 Comparison with literature models

7.1.1 Gomez-Gomar et al. (1998)

In JH98, radiative transfer of selected CO and ONe nova models are presented. They used a 1D hydrodynamics model to compute the deposition of matter and their outburst with a large 100-isotope networks up to ⁴⁰Ca. The used accretion rate was |$\dot{M} = 2\times 10^{-10}\, \mathrm{M}_{\odot }$| yr⁻¹ with f_CO = 0.5. To check that our code can reproduce qualitatively similar results, we calculate an artificial model by setting the expansion velocity to 3000 km s⁻¹ with the radioactive isotope masses in the same order of magnitude as theirs. In Fig. 13, we show the spectra at selected times and the light curve from 0.01 to 10 d.

By comparing the energetics and thermodynamics of this and our models, our model shows a comparable production of ⁷Be, ¹³N, and ²²Na. The ¹⁸F abundance is lower by one to two orders of magnitude. Our ejecta velocity is also smaller in general. Their models show a rapid expansion of about 3000 km s⁻¹, while our models show a very slow expansion right after the TNR. The ejecta reach about 300 km s⁻¹ only at ∼day 10. The peak temperature reached in their models is higher than our models by 0.2 in log₁₀ scale.

Our spectra show agreeing results with theirs. At the beginning, a flux of |$\sim 10^{-2}\, \mathrm{ph\, cm^{-2}\, s^{-1}}$| is observed in the 511 keV line. On day 1 and beyond, the 511 keV line becomes subtle while the 478 keV line becomes prominent. Meanwhile, the weak but observable 1275 keV line emerges as well. In Fig. 14, we show the spectra of our brightest models as well as the resulting spectrum from JH98 for three selected times, 15 min, 3 h, and 1 d after the explosion, and compare them to instrument sensitivities. Within a distance of 1 kpc, none of the new models would result in measurable emission.

Our light curves also show similar results compared to theirs: the early light curve is rapidly falling and then reaches an asymptotic value of about 10³⁴ erg s⁻¹ after day 1. The artificially high velocity makes the ejecta transparent at much earlier times, so that the escaping γ-rays approach the total radiative decay power much faster than in our models. A factor of two is observed between total and escaped γ-rays. Also shown in Fig. 14 is a comparison of the light curves between our models and previous calculations. Again, it is evident that current instrumentation would not be able to observe the short-duration nova flash.

7.2 Caveats

The mixing process is known to be complex for novae as the mixing changes the composition of the material (Fujimoto 1988), which affects the recurrence time, the ignition temperature, and the corresponding nuclear reactions (Shen & Bildsten 2009; Denissenkov et al. 2013). The stellar evolution model used here assumes spherical symmetry, while the mixing involving Kelvin-Helmholtz and baroclinic instabilities naturally involve modelling with 2D or 3D simulations (e.g. Casanova et al. 2016; Casanova, José & Shore 2018). Thus it is unclear whether the constant mixing and the mixing rate used in this work is completely agreeing with the actual value. The combination of using multidimensional simulations to extract the mixing parameters with application to spherical symmetric models (José et al. 2020) can provide a better estimate how the mixing takes place.

The production of ⁷Be has been in tension as theoretical models produce about 10⁻⁵ in mass fraction (Chugai & Kudryashov 2020), which is an order of magnitude lower than the averaged observational value ∼10⁻⁴ (Molaro, Cescutti & Fu 2020). This might, however, be an observational bias towards extremely large values. Our models, where ⁷Be shows a mass fraction of about 10⁻⁵ agree with the consensus of theoretical models but require further input physics to match the observational data for these extreme measurements. Recent works searching for optimized parameters suggest that initial composition, including ⁴He, still plays a role to the exact abundance pattern (Chugai & Kudryashov 2020; Denissenkov et al. 2021), but cannot completely resolve the tension. Since the ⁷Be decay has a half life ∼54 d, the exact amount will largely contribute to the later γ-ray spectra when the ejecta become transparent. The tension also suggests that more careful treatment in the pre-outburst evolution is necessary. Similar to type Ia supernovae, a single detection in the MeV γ-ray band will not necessarily clarify this conundrum because there is a variety of parameters that will influence the resulting spectrum and in turn our understanding of classical novae.

Our models predict in general the ejecta with a velocity much lower than previous works (e.g. José et al. 2001). Differences can be seen from Figs 6 and 7, where their models show a higher T_max and M_ej. A major difference comes from the choice of radiative opacity. They have used the table from Iben (1975) while ours use (Iglesias & Rogers 1993, I93). As shown in the two figures, their models show a large reduction of T_max and M_ej when the I93 opacity table is used. The more recent opacity table predicts softer matter, which can be compressed and heated more easily, leading to a TNR faster and more frequent. These effects impact to the γ-ray transport, so that it takes more time for the ejecta to become optically thin and consequently for the γ-rays to escape.

7.3 Conclusion

In this work, we have used the stellar evolution code mesa for computing an array of CO and ONe nova models for a given WD mass M_WD from 0.8 to 1.3 M_⊙, mass accretion rate |$\dot{M} = 10^{-10}$|–|$10^{-9}\, \mathrm{M}_{\odot }$| yr⁻¹, and mixing ratio f_CO = 0.1–0.5. We study the parameter dependence of the nova recurrence, thermodynamics, and chemical abundance pattern. We find that some of the nova models can fit to the extremely low values of ¹²C/¹³C, ¹⁴N/¹⁵N, and ¹⁵O/¹⁶O observed in planetary nebulae and novae.

We then carry out Monte Carlo radiative transfer simulations to study the soft γ-ray spectra and corresponding light curves. A major result of this work is that novae take a longer time (∼5–10 d) before the ejecta reach a high velocity (∼100–300 km s⁻¹) for the matter to become transparent to γ-rays. As a result, most of the 30–500 keV emission and in particular the 511 keV line coming from the short-lived radioactive isotopes such as ¹³N and ¹⁸F is heavily suppressed.

In particular, we find that the previously expected 511 keV line flash with fluxes on the order of 10⁻³–|$10^{-1}\, \mathrm{ph\, cm^{-2}\, s^{-1}}$| within the first few hours after the explosion is reduced to 10⁻¹⁶–|$10^{-9}\, \mathrm{ph\, cm^{-2}\, s^{-1}}$|⁠, depending on the model configuration. The high end of these flux expectations are only reached for high mixing values f_CO ∼ 0.5. With current instrument capabilities from transient monitors like Swift/BAT or Fermi/GBM, as well as from serendipitous observations by INTEGRAL, this early emission would be invisible for astronomical distances. Nevertheless, because there are still considerable uncertainties in these model parameters because classical novae have never been observed in the MeV γ-ray band, we propose directed and blind retrospective searches for nova outbursts in data archives. Thanks to our parametrized spectro-temporal models, the sensitivity for such searches is greatly increased, and might reveal nearby outbursts that escaped UVOIR detection.

ACKNOWLEDGEMENTS

SCL acknowledges support by NASA grants HST-AR-15021.001-A and 80NSSC18K1017. Thomas Siegert is supported by the German Research Foundation (DFG-Forschungsstipendium SI 2502/3-1) and acknowledges support by the Bundesministerium für Wirtschaft und Energie via the Deutsches Zentrum für Luft- und Raumfahrt (DLR) under contract number 50 OX 2201.

DATA AVAILABILITY

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

This project is done with the use of software mesa (Paxton et al. 2011, 2013, 2015, 2018, 2019) version 8118 and python libraries: matplotlib (Hunter 2007), pandas (The pandas development team 2020), and numpy (Harris et al. 2020).

REFERENCES