© ESO 2021

1. Introduction

Classical novae (CNe) are thermonuclear explosions on the surface of white dwarfs (WDs) and belong to the most frequent transient phenomena within a galaxy. Despite their high occurrence rate in the Milky Way of 50 ± 25 yr⁻¹ (Shafter 2017), only about 20% are seen in UVOIR wavelengths in one year. At X-ray (Ness et al. 2007) and GeV (Franckowiak et al. 2018) energies, only a few dozen objects in total could be analysed in detail to date. At MeV energies, CNe have never been observed, leaving a large gap in the understanding of nucleosynthesis in these events, as well as in the transition region between thermal and non-thermal processes.

The absolute nucleosynthesis yields from CNe are very uncertain and have only been inferred indirectly from abundance ratios in the expanding nova clouds, the latter of which generally agree very well with theoretical expectations (e.g., Livio & Truran 1994). The interiors of these objects can be studied using γ-rays, and absolute abundances can be determined (e.g., Clayton & Hoyle 1974; Clayton 1981; Leising & Clayton 1987; Hernanz & José 2006; Hernanz et al. 2014). The most promising messengers of ongoing nucleosynthesis in CNe include ⁷Be and ²²Na, which emit soft γ-rays at 478 keV and 1275 keV, respectively, as by-products of the decays to their stable daughter-nuclei ⁷Li and ²²Ne.

Individual objects are too dim for the best current γ-ray line spectrometer telescope, INTEGRAL/SPI (Winkler et al. 2003; Vedrenne et al. 2003), to be detected beyond a few hundred parsec (pc). Because many objects contribute to the total Galactic luminosity, however, we have attempted to constrain the ejecta masses of ⁷Be and ²²Na from the whole population of known and unknown CNe. The cumulative effect of this population of sub-threshold sources is expected to result in a radioactive glow along the Galactic plane in both isotopes. This effect is well measured in the case of ²⁶Al from the nucleosynthesis within massive stars (e.g., Diehl et al. 2006; Bouchet et al. 2015; Pleintinger et al. 2019). In the case of CNe, earlier measurements with COMPTEL focused on the diffuse emission of the 1275 keV line. In the inner Galaxy (|l|< 30°), Jean et al. (2001) found in an upper limit of the ²²Na ejecta mass from the diffuse emission of ONe novae (see Sect. 2.1) of < 3 × 10⁻⁷ M_⊙. The strongest constraints so far for an individual object were found by Iyudin et al. (1999), with a limit of < 2.1 × 10⁻⁸ M_⊙ from Nova Cygni 1992. A similar but rougher extended emission can also be expected from the 478 keV line, which has not been considered yet, however. The most stringent upper bound on the ejecta mass of ⁷Be was found by Siegert et al. (2018) for the CO nova V5668 Sgr, with < 1.2 × 10⁻⁸ M_⊙, assuming a distance of 1.6 kpc. All these limits are close to the highest ejecta mass estimates, for example from UV measurements (e.g., Molaro et al. 2016, finding < 0.7 × 10⁻⁸ M_⊙ for V5668 Sgr), which rely on a canonical value for the total ejecta mass, however, which itself is uncertain by one order of magnitude.

In addition to the diffuse part, more than 100 individual CNe are known to have occurred during the time of the INTEGRAL observations since 2003. Their contribution must therefore also be taken into account. This is possible in the framework of a Bayesian hierarchical model (Gelman et al. 2013) in which we attempt to link the physics of individual CNe and thus determine ejecta mass estimates from the whole CN population for both CO and ONe novae.

This paper is structured as follows: in Sect. 2 we introduce the nucleosynthesis that is expected to occur in CNe and predict individual and cumulative γ-ray signals from these expectations in general. Sect. 3 describes the INTEGRAL/SPI data set, followed by Sect. 4, which includes the sample of known CNe during the time of the observations. The general data analysis for extracting fluxes from the raw count data is shown in Sect. 5, together with the refined approach of the Bayesian hierarchical model. We present our results from individual sources, diffuse emission, and hierarchical modelling in Sect. 6. In Sect. 7 we discuss these results in terms of CN nucleosynthesis constraints and the contribution of ONe novae to the Galactic positron puzzle. We summarise and conclude in Sect. 8.

2. Nuclear astrophysics expectations for novae

2.1. Explosive burning

Nova explosions are the result of mass accretion onto a WD in a close binary system. The explosion itself is described as a thermonuclear runaway reaction of accreted material that gradually becomes degenerate, heats up under the additional pressure from still-accreting matter, and finally ignites when specific conditions are met (e.g., José & Hernanz 1998). Depending on the composition and consequently mass of the WD, different isotopes up to mass number ∼40 are produced and ejected into its surroundings (e.g., José & Hernanz 1998; José et al. 2001a; Starrfield et al. 2009, 2020). Additionally, the composition of the accreted material has to be considered as well as its fraction in the final mixture on the WD surface, resulting in a broad range of possible abundance ratios from this explosive nucleosynthesis. In particular, two types of generic WD compositions are typically considered: CO WDs in the mass range 0.6–1.15 M_⊙, and ONe WDs in the mass range above 1 M_⊙. From an independent measurement of a WD mass, the composition and CN type can thus be determined approximately (e.g., Gil-Pons et al. 2003).

In hydrodynamical models, typical accretion rates are of the order of 10⁻¹⁰ M_⊙ yr⁻¹, which results in an envelope mass of about M_env ∼ 10⁻⁵ M_⊙ during an accretion time of 10⁵ yr. The explosion sets in at a temperature around 2 × 10⁷ K and reaches a peak of several 10⁸ K about a minute after the ignition. The complete thermonuclear runaway results in the canonical total ejecta mass of $M_{\mathrm{tot}}^{\mathrm{ej}}\sim {10}^{-5}{M}_{\odot }$. In a scenario in which WDs gradually gain mass and eventually explode as type Ia supernovae, it would be expected that $M^{\mathrm{env}}>M_{\mathrm{tot}}^{\mathrm{ej}}$ (Starrfield et al. 2020). Expansion velocities are typically expected in the range of 500–3000 km s⁻¹.

UVOIR observations of CNe typically happen around maximum light when the objects are detected days to weeks after the initial explosion, as the expanding nova cloud is opaque to low-energy photons (Gomez-Gomar et al. 1998). Insight into the onset of the explosion is therefore linked to strong assumptions from theoretical modelling. During the explosive nucleosynthesis, short-lived isotopes such as ¹³N (τ₁₃ = 10 min) and ¹⁸F (τ₁₈ = 110 min) are also produced. They decay by positron emission. The positrons would quickly find electrons and lead to a strong 511 keV line during the first hour after the explosion. This has never been observed, however, mainly because this strong annihilation flash occurs days to weeks before the CN is detected. Retrospective searches for individual known objects or complete archives have not found any of these counterparts either (Skinner et al. 2008).

An alternative to complete the information in CN explosions is provided by longer-lived nucleosynthesis products. In particular, the electron-capture decay of ⁷Be (τ₇ = 76.78 d) results in an excited state of ⁷Li, which de-excites quickly by the emission of a γ-ray photon at 478 keV, so that enough ⁷Be is still present for follow-up observations after the detection of a CN. The synthesis of ⁷Be is believed to occur via ³He(α, γ)⁷Be in both CN types. The CO nova abundances of ⁷Be, however, are expected to be about one order of magnitude larger than for ONe novae (José & Hernanz 1998). In contrast, ONe novae are expected to produce great amounts of ²²Na and ²⁶Al, whereas those abundances in CO are expected to be several orders of magnitude smaller. The production of ²²Na mainly occurs through proton captures on seed nuclei, with the main reactions ²⁰Ne(p, γ)²¹Na(p, γ)²²Mg(β⁺)²²Na and ²⁰Ne(p, γ)²¹Na(β⁺)²¹Ne(p, γ)²²Na as part of the nuclear reaction network in the NeNaMgAl region. With a lifetime of τ₂₂ = 3.75 yr, ²²Na decays (90β⁺; 10% electron capture) to ²²Ne∗, which emits a γ-ray photon at 1275 keV. This allows us to search for ²²Na even years after the explosion and from multiple sources at the same time (see Sect. 2.3). The long lifetime of ²⁶Al (τ₂₆ = 1.05 Myr) results in a mixture with other more dominant ²⁶Al sources in the Galaxy, so that individual CN explosions can no longer be traced back, and only the cumulative diffuse emission remains.

We therefore restrict this study to the intermediate-lifetime isotopes ⁷Be and ²²Na. Models predict ⁷Be yields between 10⁻¹¹ and 10⁻⁸ M_⊙ for a single CO nova event (e.g., José et al. 2001a,b, 2006; Hernanz & José 2006; Starrfield et al. 2009, 2020; Hernanz et al. 2014), obtaining an upper limit for the detectability of 500 pc with INTEGRAL/SPI in the most optimistic case. The yields for ²²Na in ONe novae might reach up to 10⁻⁶–10⁻⁷ M_⊙, depending on the WD mass, composition, and mixture (Starrfield et al. 2009). However, older estimates also suggest lower ²²Na ejecta masses in ONe novae with up to 10⁻⁸ M_⊙ (e.g., José & Hernanz 1998).

During the now 18 mission years of INTEGRAL, about 100 CNe are known to have occurred inside the Milky Way (see Sect. 4), whereas about 900 are expected owing to the CN rate. Because the occurrence rate is so high compared to the decay rates of ⁷Be and ²²Na, an equilibrium mass of radioactive material inside the whole Milky Way can be expected, which supersedes that of a single CN event. In the following, we detail from predictions from first principles what to expect from γ-ray measurements.

2.2. Gamma-rays from individual objects

From the ejected mass $M_a^{\text{ej}}$ of each isotope a, the maximum γ-ray flux F_0, a at the time of the explosion T₀ from a CN at distance d can be estimated by

$$\begin{array}{c}\hfill F_{0,a}=\frac{1}{4\pi d^2}\frac{p_a^{\gamma }{M}_a^{\mathrm{ej}}}{m_a{\tau }_a},\end{array}$$(1)

where m_a is the atomic mass of isotope a, τ_a its characteristic lifetime, and $p_a^{\gamma }$ the probability of the daughter nucleus to emit a photon due to nuclear de-excitation (Thielemann et al. 2018). From the radioactive decay law, the flux as a function of time t follows an exponential decay,

$$\begin{array}{c}\hfill F_a(t)=F_{0,a}exp(-(t-T_0)/{\tau }_a)\mathrm{\Theta }(t-T_0),\end{array}$$(2)

where Θ(t − T₀) is the Heaviside function, setting the starting time (explosion date) to T₀. Given the angular resolution of SPI of ∼2.7°, individual CNe appear as point sources. The full spatio-temporal model thus reads

$$\begin{array}{c}\hfill F_a(l,b,t)=F_a(t)\delta (l-l_0)\delta (b-b_0),\end{array}$$(3)

where (l₀, b₀) are the coordinates of a CN in Galactic longitude and latitude, respectively. For the two considered isotopes in this study, ⁷Be and ²²Na, the atomic masses are m₇ = 7.017 u and m₂₂ = 21.994 u, their decay times are τ₇ = 76.8 d and τ₂₂ = 3.75 yr, and the probabilities of emitting a 478 keV and 1275 keV photon are $p_7^{\gamma }=0.1044$ and $p_{22}^{\gamma }=0.999$, respectively. The maximum flux at 478 keV and 1275 keV is therefore

$$\begin{array}{c}\hfill F_{0,7}=22887\left(\frac{M_7^{\mathrm{ej}}}{M_{\odot }}\right){\left(\frac{d}{\mathrm{kpc}}\right)}^{-2}\mathrm{ph}{\mathrm{cm}}^{-2}{\mathrm{s}}^{-1}\end{array}$$(4)

and

$$\begin{array}{c}\hfill F_{0,22}=3838\left(\frac{M_{22}^{\mathrm{ej}}}{M_{\odot }}\right){\left(\frac{d}{\mathrm{kpc}}\right)}^{-2}\mathrm{ph}{\mathrm{cm}}^{-2}{\mathrm{s}}^{-1}.\end{array}$$(5)

Clearly, individual objects can only be seen by SPI if the distance is about few hundred pc, given its nominal 3σ narrow line sensitivity¹ of 7 × 10⁻⁵ ph cm⁻² s⁻¹ (478 keV) and 5 × 10⁻⁵ ph cm⁻² s⁻¹ (1275 keV). We note that especially for ⁷Be, a sensitivity estimate or distance threshold that is based on the maximum flux becomes flawed rather quickly because the exponential decay law decreases the number of received photons faster than the significance will increase by a typical square-root of exposure-time scaling. Likewise, for ²²Na the Doppler broadening of several 1000 km s⁻¹ adds significantly to the instrumental resolution, so that the actual CN line sensitivities are worse (cf. Sect. 6). Nevertheless, the additional information of the exponential decay can be used when applying the SPI response (Sect. 5.1) to this model to predict the number of expected counts, and thus provide a more robust estimate of the ejected mass.

2.3. Gamma-rays from diffuse emission

On average, only about 20% of the expected CNe per year in the Milky Way are detected in UVOIR wavelengths (see Sect. 4). These and all undetected sources still contribute to the Galactic-wide γ-ray emission. Because the decay times of both ⁷Be and ²²Na are longer than the average waiting time between two CNe, τ_N := 1/R_N ≈ 7 d, the Galaxy in 478 and 1275 keV can be described by a diffuse glow of unseen CNe due to the radioactive build-up. The quasi-persistent γ-ray luminosity of the population of CNe in a galaxy can be estimated by the sum over all unknown individual objects n, hence

$$\begin{array}{c}\hfill L_a^{\mathrm{diff}}={\displaystyle \sum _n{L}_{a,n}^{\mathrm{PS}}=\sum _n\frac{M_a{p}_a^{\gamma }}{{\tau }_a{m}_a}{e}^{-\frac{t-T_{0,n}}{{\tau }_a}}\mathrm{\Theta }(t-T_{0,n})=p_a^{\gamma }\frac{M_a^{\mathrm{ej}}}{m_a}{R}_N.}\end{array}$$(6)

In Eq. (6), $L_{a,n}^{\mathrm{PS}}$ is the luminosity of a point source n, emitting photons from isotope a, and T_0, n = n/R_N = nτ_N is the average explosion time of each object. We note that the final expression is independent of the decay time τ_a for any isotope. This is reasonable because Eq. (6) only describes the average luminosity of an entire galaxy that produces an average mass of $M_a^{\text{ej}}$ at a rate R_N. If instead the diffuse (population) luminosity is considered as the sum of individual CN luminosities with maximum $L_{0,a}^{\mathrm{PS}}$, Eq. (6) becomes

$$\begin{array}{c}\hfill L_a^{\mathrm{diff}}=L_{0,a}{R}_N{\tau }_a,\end{array}$$(7)

which defines whether the galaxy-wide luminosity is dominated by a single source (τ_N ≫ τ_a) or the population (τ_N ≪ τ_a). Examples would be ²⁶Al mainly from massive stars and their supernovae with a decay time of τ₂₆ = 1.05 Myr compared to the core-collapse supernova rate in the Milky Way of R_CCSN = 0.02 yr⁻¹ (e.g., Diehl et al. 2006), showing that about 10⁴ supernovae contribute to the diffuse 1.8 MeV emission. Conversely, ⁴⁴Ti, which is also produced mainly in core-collapse supernovae, has a decay time of only 86 yr, so that the Milky Way in ⁴⁴Ti decay photons is dominated by one (or a few) supernova remnant(s) at each time (The et al. 2006).

The diffuse γ-ray flux can be estimated similarly, resulting in a direct conversion between a measured flux and the luminosity or the ejected mass of each object, given a known CN rate,

$$\begin{array}{c}\hfill F_a^{\mathrm{diff}}={\displaystyle \sum _n\frac{L_a^{\mathrm{diff}}}{4\pi d_n^2}=\omega L_a^{\mathrm{diff}},}\end{array}$$(8)

where $\omega ={\sum }_n{(4\pi d_n^2)}^{-1}$ is related to the effective distance of the diffuse emission, taking the probability into account that a CN occurs at a distance d_n from the Sun. The conversion factor ω can be determined in two equivalent ways by assuming a 3D density distribution of the population of CNe, which we describe in the following. Shafter (2017) showed that CNe in the Milky Way can be described as a linear combination of a De Vaucouleurs profile (ρ₁(x, y, z)) for the bulge and a doubly exponential disc (ρ₂(x, y, z)), where the weights of the components are taken as the relative CN rates in bulge and disc, f₁ = 0.1 and f₂ = 0.9. The normalised density profiles, such that ∫ dV ρ(x, y, z) = 1, are given by

$$\begin{array}{c}\hfill {\rho }_1(x,y,z)=\frac{16a^{17/2}}{2027025{\pi }^{3/2}{R}_e^3}exp(-a{\left(\frac{R}{R_e}\right)}^{1/4}){\left(\frac{R}{R_e}\right)}^{-7/8}\end{array}$$(9)

and

$$\begin{array}{c}\hfill {\rho }_2(x,y,z)=\frac{1}{4\pi r_e^2{z}_e}exp(-\frac{r}{r_e})exp(-\frac{|z|}{z_e}),\end{array}$$(10)

with $R=\sqrt{x^2+{\mathit{y}}^2+z^2}$, a = −7.669, R_e = 2.7 kpc, $r=\sqrt{x^2+{\mathit{y}}^2}$, r_e = 3.0 kpc, and z_e = 0.25 kpc. The distribution of distances can now be directly sampled from the density profiles, for example via 3D rejection sampling of (x, y, z)-coordinates, from which the distances are calculated as $d=\sqrt{{(x-x_s)}^2+{(\mathit{y}-{\mathit{y}}_s)}^2+{(z-z_s)}^2}$, where (x_s, y_s, z_s) = (8.179, 0, 0.020) kpc is the position of the Sun. We show 10⁴ samples of distances calculated from the distribution ρ_tot = f₁ρ₁ + f₂ρ₂ in Fig. 1. We find that this distribution can be adequately approximated by a Γ-distribution with α_d = 4.25 and β_d = 1/2.45, with the expectation value α_d/β_a = 10.4 kpc. We use this distribution of CN distances later in Sect. 4.2 to define a prior for objects with unknown distances. With a large sample size, the infinite sum in Eq. (8) converges to ω = 0.00291 kpc^-2, or to an effective distance of d_eff = (L^diff/(4πF^diff))^1/2 = 5.23 kpc.

Fig. 1. Distance distribution of CNe according to Eqs. (9) and (10) through rejection sampling. Shown are 104 samples (grey histogram, corresponding to a timescale of 200 yr, left axis), and an approximation with the Γ-distribution in red. The cumulative distributions are shown in black for the sample and in cyan for the Γ-distribution (right axis).

As an alternative to estimating the infinite sum, we calculated the diffuse emission map exactly by line-of-sight integration,

$$\begin{array}{c}\hfill F^{\mathrm{diff}}(l,b)=\frac{1}{4\pi }{\displaystyle {\int }_0^{+\infty }\mathrm{d}s\rho (x^{\prime }(s),y^{\prime }(s),z^{\prime }(s)),}\end{array}$$(11)

where (x′(s),y′(s),z′(s)) = (x_s − scos(l)cos(b),y_s − ssin(l)cos(b)), z_s − ssin(b)) is the line-of-sight vector starting from the Sun in all directions (l, b) (cf. also Erwin 2015). The total integrated flux of this map is consequently

$$\begin{array}{c}\hfill F^{\mathrm{diff}}={\displaystyle \int }\mathrm{d}\mathrm{\Omega }{F}^{\mathrm{diff}}(l,b).\end{array}$$(12)

The intrinsic luminosity of a given density distribution is given by

$$\begin{array}{c}\hfill L^{\mathrm{diff}}={\displaystyle \int \mathrm{d}\mathrm{\Omega }{\int }_0^{+\infty }\mathrm{d}s{s}^2\rho (x^{\prime }(s),y^{\prime }(s),z^{\prime }(s)),}\end{array}$$(13)

so that the conversion factor reads

$$\begin{array}{c}\hfill \omega =\frac{1}{4\pi }\frac{\int \mathrm{d}\mathrm{\Omega }{\int }_0^{+\infty }\mathrm{d}s\rho (x^{\prime }(s),y^{\prime }(s),z^{\prime }(s))}{\int \mathrm{d}\mathrm{\Omega }{\int }_0^{+\infty }\mathrm{d}s{s}^2\rho (x^{\prime }(s),y^{\prime }(s),z^{\prime }(s))}.\end{array}$$(14)

Diffuse quasi-persistent γ-ray emission from the population of CNe in the Milky Way is therefore

$$\begin{array}{c}\hfill F_a^{\mathrm{diff}}(l,b)=L_a^{\mathrm{diff}}\omega F^{\mathrm{diff}}(l,b).\end{array}$$(15)

In Fig. 2 we show the diffuse emission template and overlay all individual CNe that we considered for this study (Sect. 4). For the two isotopes, the diffuse flux can finally be related to the ejected mass by

$$\begin{array}{c}\hfill F_7^{\mathrm{diff}}=178\left(\frac{M_7^{\mathrm{ej}}}{M_{\odot }}\right)\left(\frac{R_N}{{\mathrm{yr}}^{-1}}\right)\mathrm{ph}{\mathrm{cm}}^{-2}{\mathrm{s}}^{-1}\end{array}$$(16)

Fig. 2. Diffuse emission template from line-of-sight-integrated nova density distribution, Eq. (15), together with nova sample (cyan points). The contours indicate the De Vaucouleurs profile for the bulge and the exponential disc.

and

$$\begin{array}{c}\hfill F_{22}^{\mathrm{diff}}=528\left(\frac{M_{22}^{\mathrm{ej}}}{M_{\odot }}\right)\left(\frac{R_{\mathrm{ONe}}}{{\mathrm{yr}}^{-1}}\right)\mathrm{ph}{\mathrm{cm}}^{-2}{\mathrm{s}}^{-1},\end{array}$$(17)

where $R_{\mathrm{ONe}}\approx \frac{1}{3}{R}_N$ is the ONe nova rate (Gil-Pons et al. 2003). Using Eqs. (16) and (17) and the theoretical expectations, Sect. 2, we expect a diffuse flux for the 478 keV line of about 10⁻⁸–10⁻⁵ ph cm⁻² s⁻¹. In the energy region of the ⁷Be decay line, there is strong Galactic background emission, mainly from the ortho-positronium continuum and inverse Compton scattering (e.g., Churazov et al. 2005, 2011; Jean et al. 2006; Bouchet et al. 2010; Siegert et al. 2016, 2019a), which we take into account when the ejected ⁷Be is calculated in a later step (Sect. 6.2). Even though the expected fluxes for individual objects as well as the diffuse emission of the 478 keV line barely scratch the sensitivity of SPI, a combined fit that takes all objects into account that share similar physics can provide stronger limits on the ejected masses. Likewise, the 1275 keV line from ²²Na can be expected to show a diffuse flux of about 10⁻⁵–10⁻⁴ ph cm⁻² s⁻¹. This is well within the sensitivity threshold of SPI. However, the instrumental background line at 1275 keV does not follow the variation of cosmic-ray intensity defined by the solar cycle strictly, but builds up as a function of mission time and thus increases the background at these energies (Diehl et al. 2018).

2.4. Expectations from cumulative signals

From these considerations and theoretical expectations, we compiled a first-order population synthesis model for both isotopes and compared this to the expected background of our chosen data set (see Sect. 3) to obtain a signal-to-noise ratio (S/N). In both cases, we sampled 3D positions according to the combined 3D density distributions, Eqs. (9) and (10), from which distances and Galactic coordinates were calculated. The occurrence of CNe in the Milky Way is a Poisson process (Ross 2008) with a rate of R_N ≈ 50 yr⁻¹ for CO novae and $\approx \frac{1}{3}{R}_N$ for ONe novae. This means that the waiting time for each CNe after the last one is exponentially distributed. We therefore sampled the waiting times as ΔT ∼ Exp(R_N). This provided the explosion date of event n by $T_{0,n}=\sum _i^n\mathrm{\Delta }{T}_i$. For the logarithm of ejected masses of ONe novae, we assumed a normal distribution $lg{M}_{22}^{\mathrm{ej}}\sim \mathcal{N}(-8.25,0.5^2)$ and for CO novae $lg{M}_7^{\mathrm{ej}}\sim \mathcal{N}(-10,1^2)$, presenting conservative estimates from theoretical expectations (Sect. 2.1). The widths of the distributions were chosen to cover the plausible range of theoretical ejecta masses and to not produce unreasonably high γ-ray fluxes that would have already been seen by previous instruments. In a more realistic version of this population synthesis model, the ejecta masses would be distributed according to the mass distribution of WDs on which CNe occur. Because both distributions are fairly uncertain, we considered the general case and only included the range of theoretical models from Starrfield et al. (2009, 2020). Finally, fluxes at the sampled positions were calculated via Eqs. (4) and (5).

We show the characteristic time profile for one realisation in both lines integrated across the whole sky in Fig. 3. Clearly, the flux of the 478 keV line appears more erratic because the lifetime of ⁷Be is shorter than that of ²²Na. Because the lifetime of ⁷Be is still longer than the average waiting time of about one week, however, the total integrated flux averages to 2–11 × 10⁻⁵ ph cm⁻² s⁻¹ at any given time. Large peaks, such as around year 18 and 116, are due to individual objects that occurred close to the observer. The diffuse part of the 478 keV emission can be estimated from this to be about 10⁻⁵ ph cm⁻² s⁻¹, whereas at each given time, the total flux is dominated by about ten individual objects (cf. R_Nτ₇ ≈ 10.5). For ²²Na, the decreased ONe nova rate yields a quasi-persistent flux between 4 and 9 × 10⁻⁵ ph cm⁻² s⁻¹. In the case of the 1275 keV line, the diffuse flux is dominant, around 5 × 10⁻⁵ ph cm⁻² s⁻¹, and individual objects only contribute significantly if they are close to the observer, such as the event in year 116.

Fig. 3. Population synthesis of nova events in a Milky Way-like galaxy. Left: year 116.3 of the modelled Poisson process, where a nearby CO nova at (l, b) = (+123° , + 8° ) outshines the remaining galaxy in 478 keV emission (top). At the same time, 1275 keV emission is not enhanced (bottom). The angular resolution is chosen according to the SPI 2.7° (FWHM). Right: time profiles (light curves) of the 478 keV and 1275 keV since the beginning of the synthesis for a duration of ∼150 yr. This timescale is enough to reach convergence and to show characteristic features. The vertical purple line indicates the time shown in the plots on the left.

We estimated the range of possible S/N, given these assumptions, using a generic background model (Sect. 5) and 1000 realisations of our CN sample (Sect. 4) and the diffuse emission (Eq. (15)) in Fig. 4 for the case of ²²Na. In both cases, individual objects are most of the time below the detection threshold of SPI. However, the cumulative signal of 97 targets that are known to emit at these wavelengths result in an average S/N of ∼5σ for ⁷Be and ∼3σ for ²²Na. In the diffuse-emission case, the 478 keV line would be seen on average with a S/N of ∼1σ and the 1275 keV line with ∼3σ. Therefore both emission lines would show a S/N of ∼5σ when this INTEGRAL data set is considered. However, the allowed range of theoretical ejecta masses leads to considerable uncertainties, being consistent with undetectable signals for SPI even after 16 years. Before the launch of INTEGRAL, Jean et al. (2000) performed a similar study to estimate the chances for SPI to detect the 1275 keV line and suggested that measurements of the diffuse emission could probe ejecta masses ≲10⁻⁸ M_⊙. Here we show that these pre-launch simulations were too promising.

Fig. 4. Estimated S/N given the SPI data set (Sect. 3) and the range of plausible 22Na ejecta masses for our sample of 97 novae. Top: the S/N of individual objects does not increase with $\sqrt{T_{\exp }}$ because the exponential decay is faster than the accumulation of 1275 keV photons (cyan lines). The cumulative signal of point sources saturates on a timescale of a few years (dark blue). If all objects are assumed to be ONe novae, the S/N would be about three times higher (pale red). Bottom: comparison between cumulative point-source emission significance and diffuse-emission significance during the INTEGRAL mission timescale. The bands contain the 68th percentile from 1000 realisations of the population synthesis model (boundaries marked by dashed lines, median by solid lines; Sect. 2.4), to distribute ejecta masses among sources and diffuse emission.

3. SPI data sets

We used publicly available SPI data between 2003 and 2018, including INTEGRAL revolutions 43–1951. Except for a few exceptions, such as V5668 Sgr (Siegert et al. 2018), INTEGRAL did not purposely target CN outbursts, so that the exposure times of our selected sample (see Sect. 4) often contain large gaps or objects are only observed in the partially coded field of view. In a first step, we applied selection criteria based on orbital parameters and instrumental count rates and sensors: Until 2015, we selected orbital phases 0.1–0.9 to avoid the Van Allen radiation belts. Afterwards we only selected orbital phases 0.15–0.85 because of an orbit manoeuvre. We excluded times when the count rate of the SPI anticoincidence shield was increased significantly, for example due to solar flares. We furthermore excluded times when the cooling plate temperatures showed a difference larger than 1 K. Revolutions 1554–1558 were removed from the data set because of the outburst of the microquasar V404 Cygni, which was the brightest soft γ-ray source in the sky during this time. Revolutions in which a detector failure occurred were completely excluded as well.

In a second step, we fit this data set with our instrumental background model and the diffuse emission map (Sect. 5) and removed all pointed observations that showed residuals larger than 10σ. This removed about 0.1% of all selected pointings. We show the resulting exposure map in Fig. 5.

Fig. 5. Exposure map of the selected data set. The effective areas around 478 and 1275 keV are ∼80 and ∼55 cm^2, respectively.

3.1. ⁷Be

The electron capture of ⁷Be results in a photon at 477.62 keV from an excited state of ⁷Li ($p_7^{\gamma }=0.1044$). We used all SPI single events, that is, photons that only interact once with the detector, resulting in 100 105 pointed observations, with a typical exposure time between 1800 and 3600 s. The total dead-time-corrected exposure time, taking failed detectors into account, is 180.1 Ms.

The spectral resolution at 478 keV is about 2.1 keV (Diehl et al. 2018). Additional astrophysical Doppler broadening of the line can be expected from the expansion velocity of the nova ejecta, which typically is in the range of 500–3000 km s⁻¹ (see Sect. 2.1). To include more than 99% of the γ-ray line photons in a single bin for an ejecta velocity of 2000 km s⁻¹, we used an 8 keV wide energy bin between 474 and 482 keV. This includes 7.0 keV broadening due to homologous ejecta, and was added in quadrature to the instrumental resolution.

3.2. ²²Na

The decay of ²²Na to ²²Ne, either through β⁺ decay ($p_{22}^{+}=0.904$) or electron capture, results in a γ-ray photon of 1274.58 keV with $p_{22}^{\gamma }=0.999$. At this energy range, electronic noise features in the SPI spectrum can emerge, which can be filtered through pulse shape discrimination (PSD). The selection for PSD events reduces the size of the data set to 99 880 pointings and the exposure time to 179.6 Ms. We note that the fluxes resulting from PSD events have to be re-scaled by an efficiency correction factor ≈1/0.85. The SPI spectral resolution around 1275 keV is about 2.8 keV, so that we performed our analysis in a single 20 keV energy bin between 1265 and 1285 keV to account for possible astrophysical Doppler broadening (as above, Sect. 3.1).

4. Nova sample

4.1. Selecting from databases and catalogues

In order to infer information about the ejecta masses from flux measurements, the distance to the objects and their unambiguous classification as CNe is important. While many astrophysical transients are called ‘novae’, we clearly wished to include only the thermonuclear explosions on the surface of WDs, which are called ‘classical novae’. For current soft γ-ray instrumentation, it is hopeless to detect any CN outside the Milky Way, so that we restricted our selections to our own Galaxy, that is, to objects with galactocentric distances of ≲25 kpc (cf. Fig. 1). This is justified considering the expected quasi-persistent 478 and 1275 keV line fluxes from M31 or the Large Magellanic Cloud (LMC): Using the same formalism from Sect. 2.4 with adjusted CN rates (R_N(M31) = 100 yr⁻¹, R_N(LMC) = 2 yr⁻¹) and distances (d(M31) = 780 kpc, d(LMC) = 50 kpc), we expect at most 10⁻⁷ ph cm⁻² s⁻¹ (1275 keV) and 9 × 10⁻⁹ ph cm⁻² s⁻¹ (478 keV) from the LMC and 2 × 10⁻⁸ ph cm⁻² s⁻¹ (1275 keV) and 2 × 10⁻⁹ ph cm⁻² s⁻¹ (478 keV) from M31.

Several comprehensive and well-maintained CN catalogues are available online, but they are scattered in multiple publications and have different information. We first selected CNe from recent peer-reviewed articles, in particular Özdönmez et al. (2018), Shafter (2017), Hounsell et al. (2016, SMEI), Saito et al. (2013, VVV), and Walter et al. (2012, SMARTS). We complemented the information from individual websites that host and update CN lists. These include Koji Mukai’s ‘List of Recent Galactic Novae’ (2008–2020)², the CBAT ‘List of Novae in the Milky Way’ (seventeenth century until 2010)³, Bill Gray’s ‘List of Galactic Novae’ (‘Project Pluto’; seventeenth century until 2019)⁴, Christian Buil’s ‘Nova Corner’ (1999–2015)⁵, and the ARAS Spectral Data Base of Novae (2012–2020)⁶.

INTEGRAL launched on October 17, 2002 (Winkler et al. 2003), and publicly accessible data are available starting with INTEGRAL revolution 43 (MJD = 52683; February 13, 2003). Owing to the lifetime of ²²Na of 3.75 yr, it is also reasonable to include CNe before the launch and first measurements of INTEGRAL. Given the shorter lifetime of ⁷Be of 0.21 yr and to make a coherent sample, we included all objects from January 1, 2002 (MJD = 52275), to June 30, 2018 (MJD = 58299), that is, 16.5 yr, and those that were at least in the partially coded field of view of SPI (∼30° corner to corner) for 100 ks.

We included all objects whose type was a CN called either ‘fast’, ‘moderately fast’, ‘slow’, or other. This also excluded nova-named objects, such as V838 Mon (Bond et al. 2003) and similarly classified events. The absolute magnitude in UVOIR wavelengths is no selection criterion for this study. We did not include known recurrent CNe in our sample as they would require additional special treatment considering their long-term light curves and ejecta distribution. This further excludes in particular U Sco (the last two known outbursts occurred in 1999 and 2010), RS Oph (1985 and 2006), and IM Nor (1920 and 2002).

In total, this provides a sample of known objects within the considered time frame of 97 CNe. Given the CN rate of ∼50 yr⁻¹, about ∼12% of all expected objects are included as point sources with known positions in this study (Sect. 2.2). The remaining ∼88% of unknown sources plus hundreds of CNe from before the start of the considered time frame then make up the diffuse component of the γ-ray emission (Sect. 2.3). The full list of objects can be found in Table B.1.

4.2. Handling unknown distances

For 44 CNe in this sample – mainly older objects –, a distance estimate is available from the different input catalogues and websites. When two or more estimates were available for one object, we used the most recent value. Most of the distance estimates for CNe come from the maximum magnitude relation with decline time (MMRD, Zwicky 1936; McLaughlin 1940; Buscombe & de Vaucouleurs 1955). For objects located inside a certain cluster or galaxy, this method provides reliable estimates for the population distance. However, individual CNe inside the Milky Way still carry uncertainties on their distances of about 30–50%. For example, the distance estimate of 3 ± 1 kpc for V5115 Sgr (Table B.1) would suggest an ejecta mass uncertainty of at least 66% (≈0.2 dex), ignoring the flux uncertainties.

The problem of converting estimated fluxes into ejecta masses for objects with unknown distance is even more severe: When we consider the distance distribution of CNe inside the Milky Way from the point of Earth, Fig. 1, it is clear that most sources are expected around the Galactic centre. However, an object found toward the Galactic anticentre with unknown distance will most certainly not be at such a large distance; these distances are still possible, but only with a probability P(7 < d < 9) ≈ 17%. Very UVOIR-bright CNe could be located close to Earth, but with an unavailable distance estimate: The probability that any object is located inside a sphere of radius 0.5 kpc (1.0 kpc) is P(d ≤ 0.5) ≈ 3 × 10⁻⁵ (5 × 10⁻⁴). Given the CN rate, less than one object is expected within these distances during the INTEGRAL observations, consistent with no γ-ray detection so far. These considerations make the distribution in Fig. 1 and its approximation with d_unknown ∼ Γ(α_d = 4.25, 1/β_d = 1/2.45) a reasonable first-order distance estimator for any CNe inside the Milky Way, independent of its direction as seen from Earth. We use the approximated distribution as prior information for inference in Sect. 5.2.1 to construct a full posterior for the ejecta masses that properly takes all uncertainties into account.

We estimated the effect of using this generic distribution by naive scaling factors as above: As described in Sect. 2.3, the expectation value of distances is ⟨d_unknown⟩ = 10.4 kpc, with a standard deviation σ_dunknown = 5.1 kpc. Half of all CNe observed from Earth are found within a sphere of 9.6 kpc, thus including the Galactic bulge, nearby high-latitude objects, and the Galactic anticentre, and thus also covering most known CN distances (Table B.1). When we consider the symmetric 90% interval around the expectation value, a canonical distance uncertainty of 7.6 kpc can be given. This makes the ejecta masses of objects with unknown distances uncertain by at least 150% (≈0.4 dex). Given the large total volume of the Milky Way, this appears as a reasonable estimate.

5. SPI data analysis

5.1. General method

The SPI data d_jpe are detector (j) triggers per unit time (pointing p) and energy (e), and consequently follow the Poisson distribution,

$$\begin{array}{c}\hfill d_{\mathit{jpe}}\sim \mathrm{Poisson}(m_{\mathit{jpe}})=\frac{m_{\mathit{jpe}}^{d_{\mathit{jpe}}}{e}^{-m_{\mathit{jpe}}}}{d_{\mathit{jpe}}!},\end{array}$$(18)

where m_jpe is the (modelled) rate parameter. In our analysis, the energy is fixed and we omit the index e in the following. The likelihood of measuring d_jp counts in the selected data set D, given a model M with expectation m_jp, is

$$\begin{array}{c}\hfill \mathcal{L}(D|M)={\displaystyle \prod _{\mathit{jp}}\mathrm{Poisson}(d_{\mathit{jp}}|m_{\mathit{jp}}).}\end{array}$$(19)

We modelled the SPI counts as a linear combination of sky components and background models,

$$\begin{array}{c}\hfill m_{\mathit{jp}}={\displaystyle \sum _{i=1}^{N_N+1}{F}_{i,p}{R}_{\mathit{jp}}^{\mathit{lb}}{S}_{i,lb}+\sum _{k=L,C}{\beta }_k{B}_{k,jp}.}\end{array}$$(20)

In Eq. (20), N_N is the number of CNe in our sample, plus one model component to account for diffuse emission. F_i, p is the flux of sky model S_i, lb to which the imaging response function $R_{\mathit{jp}}^{\mathit{lb}}$ was applied (mask coding). The background model B_k, jp was split into two components, k = L and C, accounting for instrumental line and continuum background, respectively. Temporal variations in the background were determined by the parameters β_k, scaling the amplitudes of the two components (cf. Siegert et al. 2019b). The expected flux of point-like sources (individual CNe) was determined according to Eqs. (4) and (5). Because the exponential decay times are known, the pointing-to-pointing variation can be fixed in Eq. (20), and only one parameter was fitted for each source and isotope. The expected diffuse flux was modelled according to Eqs. (16) and (17). In the individual CN case, the spatial model is one point at the coordinates of the CN. For diffuse emission, the flux is constant in time and distributed according to the line-of-sight-integrated density structure Eq. (15). Finally, the fit parameters of interest are the maximum flux F at the time of explosion T₀, which can be split into a function of ejecta mass M^ej and distance d for point sources, and mass and CN rate R_N for diffuse emission.

Clearly, this would result in an ill-defined likelihood function if the distances and the CN rate are not constrained. In the case of a known distance for CNe in our sample, we incorporated a normal prior on the distance. When the distance to an object was not known, we used the information that the CN occurred inside the Milky Way, which translates into a distance prior according to a Γ-distribution d_unknown ∼ Γ(α_d = 4.25, β_d = 1/2.45) (see Sects. 2.3 and 4.2). This sets any object for which the distance is unknown on average to α_d/β_d = 10.4 kpc with a standard deviation ${({\alpha }_d/{\beta }_d^2)}^{1/2}\approx 5.1$ kpc (cf. Fig. 1). For the CN rate, we set a prior according to the most recent literature value in units of yr⁻¹ of R_N ∼ 𝒩₊(50, 25²) (index +: truncated at zero) for all CN types and $R_{\mathrm{ONe}}=\frac{1}{3}{R}_N$ for ONe novae (Shafter 2017). Especially the prior for the CN rate is particularly broad and can, in principle, be consistent with zero. A zero rate, however, is unphysical because CNe constantly occur throughout the year. We nevertheless included this extreme to show a broad range for the remaining parameter space. We set a uniform prior on the logarithm of the ejecta mass, lgM^ej ∼ 𝒰(−11, −4), so that each decade considered in the parameter space obtained the same prior probability, and a wide range of theoretical and otherwise plausible values were sampled. The lower bound on the ejecta mass prior of 10⁻¹¹ M_⊙ is data driven, as above considerations (Sect. 2.4) show that SPI is incapable of probing these values except for distances smaller than 100 pc.

Details on the background modelling procedure can be found in Diehl et al. (2018) and Siegert et al. (2019b). In short, the background is modelled by long-term monitoring of the complete SPI spectrum and the spectral response changes. This results in a background and response database per INTEGRAL orbit and detector for several 100 background lines, and the underlying continuum. Because the background dominates the measured count rate at any time, and in addition, the coded mask pattern from celestial objects smears out on a timescale of tens of pointings (typically 50–100 in one INTEGRAL orbit), a background response can be constructed from this database. The short-term pointing-to-pointing variation is fixed by an onboard counting rate, in our case, the saturating germanium detector events. Finally, any variation that is not captured by this procedure is handled by the introduction of the background re-scaling parameters β_k. In both cases, we find that one parameter per INTEGRAL orbit per background component is enough to provide an adequate fit. The total number of fitted parameters for one object finally is twice the number of INTEGRAL orbits, plus two for the flux, that is, n_par = 2 ⋅ 1674 + 2 = 3350. For the 478 keV line (1275 keV), the number of data points excluding dead detectors is 1598964 (1595430).

5.2. Hierarchical modelling approach

In order to combine the extracted posterior distributions of the flux, we applied a hierarchical model (Gelman et al. 2013) in three different steps. In Fig. 6 we show the full graphical model as separated into interesting parameters, nuisance parameters, fixed parameters, resulting spatio-temporal models, and their interdependencies. When we consider the sky emission from top to bottom, the hierarchy flows from superordinate hyper-parameters μ and τ, which describe the distribution of ejecta masses M_i of each CN i, which then determines the fluxes, given additional parameters: Depending on the distances d_i, or in the diffuse case, the CN rate R_N, absolute fluxes F_i and F_diff result from Eqs. (4), (5), (16), and (17). These are either point-like at source positions (l, b)_i or distributed according to the line-of-sight integration of ρ (cf. Eqs. (9), (10), and (15)), resembling the sky models $F_i^{\mathit{lb}}$ and $F_{\mathrm{diff}}^{\mathit{lb}}$, respectively. For individual objects, the discovery times T_0, i are known⁷, after which the flux is exponentially decaying according to isotope-specific parameters (Eq. 1). The diffuse emission component is constant in time. In particular, this scheme is included in our Poisson process, Sect. 2.4. We now try to infer the values of M_i and/or μ and τ.

Fig. 6. Complete graph of different model variants used in this study. When the source models are used individually, that is, regardless of a common underlying physical process, the most conservative estimates are found in a no pooling setting (Sect. 5.2.1). The opposite extreme assumes that all objects eject the same amount of matter, so that instead of ∼100 different masses, only one is fitted. This resembles complete pooling and defines tight constraints but is subject to bias if the true variance is large (Sect. 5.2.2). A compromise between no and complete pooling is achieved by invoking the hyper-parameters μ and τ that determine the mean and the spread of ejecta masses, if any (partial pooling, dash-dotted box, Sect. 5.2.3).

Ideally, all sky models are convolved with the imaging response to convert them into a SPI-compatible format and are then fitted simultaneously to the raw SPI count data, together with the two-component background model. Given the large number of data points, fitted parameters, and especially their interdependence in the full hierarchy (see Sects. 5.2.1–5.2.3), this is computationally very expensive and not feasible without very many high-performance computing hours. Instead, we applied the hierarchical model in a two-step process:

First, we extracted the fluxes of each of the N_N = 97 CNe and the diffuse emission individually (Eq. (20)). This resulted in posterior distributions for the fluxes of all sky models considered. To first order, these posteriors follow a Γ-distribution⁸, so that instead of the previous 1.5 million data points, we condensed the information into the two shape parameters, α and 1/β of the Γ-distribution for each sky model. Thus, in a second step, we applied the hierarchical model to a reduced data set consisting of (97 + 1) ⋅ 2 shape parameters.

Because the actual measurement information (Poisson counts) is lost in this procedure, the resulting combined values are approximations to the true distribution. We already note here that none of the sources considered show a significance above background of more than 2.5σ, so that most flux posteriors are sharply peaked at zero flux, with a long tail according to the Γ-distribution.

5.2.1. No pooling

The most agnostic view about the physics at play is provided by treating each object individually (no pooling). This means that we let the data speak for themselves if uninformative or weakly informative priors are used, as is the case for our general method: The flux extraction step, Sect. 5.1, is equivalent to a no pooling analysis. Because we are interested in the ejecta masses, we fit the two shape parameter of the flux posteriors, in particular their expectation values ⟨F_i⟩=α_i/β_i, according to the likelihood

$$\begin{array}{c}\hfill P(⟨F_i⟩|{\beta }_i,M_i^{\mathrm{ej}},d_i,R_N)={\displaystyle \prod _{i=0}^{N_N+1}\mathrm{\Gamma }(⟨F_i⟩;\stackrel{\sim }{F_i}{\beta }_i,{\beta }_i),}\end{array}$$(21)

where $\stackrel{\sim }{F_i}$ are the latent fluxes according to Eqs. (4), (5), (16), and (17). By definition, this results in a distribution for the flux of each object that is consistent with the input distribution and an expectation value $⟨F_i⟩=\stackrel{\sim }{F_i}$. By adding the distance (or CN rate) information in a prior, we obtain the posterior distribution of ejecta masses for each individual object (and diffuse emission) independently, according to the measured flux and its asymmetric uncertainties. This determines 98 values for $M_i^{\text{ej}}$ from 98 data points and their 98 uncertainties. For completeness, we repeat the prior probabilities that are used here and elsewhere in the paper for the same parameters:

$$\begin{array}{cc}& P(d_{\mathrm{known}})\sim {\mathcal{N}}_{+}({\mathrm{\mu }}_d,{\sigma }_d^2)\hfill \\ \hfill & P(d_{\mathrm{unknown}})\sim \mathrm{\Gamma }({\alpha }_d,{\beta }_b)\hfill \\ \hfill & P(lg{M}^{\mathrm{ej}})\sim \mathcal{U}(lg{M}_{min}^{\mathrm{ej}},lg{M}_{max}^{\mathrm{ej}})\hfill \\ \hfill & P(R_N)\sim {\mathcal{N}}_{+}({\mathrm{\mu }}_{R_N},{\sigma }_{R_N}^2)\hfill \\ \hfill & P(R_{\mathrm{ONe}})\sim {\mathcal{N}}_{+}({\mathrm{\mu }}_{R_{\mathrm{ONe}}},{\sigma }_{R_{\mathrm{ONe}}}^2),\hfill \end{array}$$(22)

with μ_d and σ_d from Table B.1 in units of kpc, α_d = 4.25 and β_d = 1/2.45 such that α_d/β_d = 10.4 kpc is the expectation value of CNe with unknown distances, $M_{min}^{\mathrm{ej}}={10}^{-11}{M}_{\odot }$ and $M_{max}^{\mathrm{ej}}={10}^{-4}{M}_{\odot }$, μ_RN = 50 yr⁻¹ and σ_RN = 25 yr⁻¹, and ${\mathrm{\mu }}_{R_{\mathrm{ONe}}}=16.7{\mathrm{yr}}^{-1}$ and ${\sigma }_{R_{\mathrm{ONe}}}=8.3{\mathrm{yr}}^{-1}$. Instead of using truncated normal distributions for the event rates, we also tested log-normal distributions with the same means and variances for the hierarchical model fits. This avoids the unphysical boundary of a (nearly) zero rate. The posterior distribution of the ejecta masses in each object, marginalised over the nuisance parameters distance and rate, finally reads

$$\begin{array}{c}\hfill P(M_i^{\mathrm{ej}}|⟨F_i⟩,{\beta }_i)\propto {\displaystyle \int \int }d{d}_{i}dRP(⟨F_i⟩|{\beta }_i,M_i^{\mathrm{ej}},d_i,R)P(M^{\mathrm{ej}},d_i,R).\end{array}$$(23)

In Eq. (23), P(M^ej, d, R) is the joint prior distribution, built as multiplicative distribution from Eq. (22).

We note that neither the true Poisson count data nor the condensed data set can constrain the distance (or rate) information, which consequently is only a nuisance parameter and included to properly convert the flux values into ejecta masses, given the uncertainties on the distances (or rate). Of our 97 objects, certainly not all are ONe novae. Nevertheless, we determined the upper bounds for the 1275 keV line flux and the resulting ²²Na ejecta mass bounds from these objects. A more detailed discussion of the results and limits from individual objects is presented in Sect. 6.1. The complete results table can be found in Table B.1.

5.2.2. Complete pooling

One way to improve upon estimates from individual objects – which in our case are rarely constraining (Sect. 6.1) – is to include the whole population and assume that each object (and the cumulative diffuse emission) stems from the same ejected mass. When we consider Fig. 6, this means that $M_i=\stackrel{\sim }{M}$, and the resulting fluxes ⟨F⟩_i (and ⟨F⟩_diff) share a common dependence in the complete pooling setting. This is sometimes called ‘stacking’ (e.g., Malz 2021), which is not equivalent in detail, however, because the extracted spectra (or fluxes) are not stacked to obtain an average spectrum, but the model parameters are shared among all objects. For a completely pooled (combined) posterior for the ejected mass, we thus sampled the fluxes according to the likelihood

$$\begin{array}{c}\hfill P(⟨F_i⟩|{\beta }_i,{\stackrel{\sim }{M}}^{\mathrm{ej}},d_i,R_N)={\displaystyle \prod _{i=0}^{N_N+1}\mathrm{\Gamma }(⟨F_i⟩;\stackrel{\sim }{F_i}{\beta }_i,{\beta }_i).}\end{array}$$(24)

The only difference between Eqs. (24) and (21) is that there is only one parameter of interest, ${\stackrel{\sim }{M}}^{\mathrm{ej}}$, which is then representative for the whole population. This determines one value for M^ej from 98 data points and their 98 uncertainties, which are dependent on each other, and it still includes all uncertainties from distances, rates, and individual flux measurements.

Complete pooling results in the most optimistic but also most biased estimate of the ejecta masses. Because the known effects of different ejecta masses for different WD masses and compositions are ignored in this setting, the resulting estimate is equivalent to the mean of the population of CNe. This mean is not necessarily a useful number because it is dominated by the most abundant objects, most frequent CN types, or observations with the smallest uncertainties.

5.2.3. Partial pooling

Because of the different astrophysical and measurement-related effects, the true ejecta masses of CNe in the Milky Way may follow a distribution that can be described by a mean value μ and a width τ, such that

$$\begin{array}{c}\hfill lg{M}_i^{\mathrm{ej}}\sim \mathcal{N}(\mathrm{\mu },{\tau }^2).\end{array}$$(25)

This means that we first sampled ejecta masses according to the distribution, Eq. (25), and then again sampled the flux distributions as before, but now the distributions depended on the hyper-parameters μ and τ rather than $M_i^{\mathrm{ej}}$ individually. The likelihood becomes

$$\begin{array}{c}\hfill P(⟨F_i⟩|{\beta }_i,\mathrm{\mu },\tau ,d_i,R_N)={\displaystyle \prod _{i=0}^{N_N+1}\mathrm{\Gamma }(⟨F_i⟩;\stackrel{\sim }{F_i}{\beta }_i,{\beta }_i).}\end{array}$$(26)

This partial pooling setting is a compromise between no pooling and complete pooling: depending on the value of τ, the samples of ⟨F_i⟩ either converge to one common value (complete pooling, τ → 0) or spread out towards their independent values (no pooling, τ → ∞). The prior for the mean of this hyper-distribution was chosen similar to that of the individual masses, μ ∼ 𝒰(−11, −4). The prior for the width τ is discussed in the literature (Gelman et al. 2013, and references therein), depending on the purpose of the analysis. We tested two priors for τ: first, a uniform distribution between 0 and 10, which is agnostic about the fact that each object should indeed rather be similar in ejecta mass and should not show a difference of several orders of magnitude. However, especially in the case of the 1275 keV emission, such a prior would be useful because the ejecta masses would indeed be either close to zero (CO novae) or around a certain mean (ONe novae). Second, we tested a half-Cauchy prior for τ with width 1. This assumes that most objects follow a similar trend, but also allows for a broader distribution that might include several outliers.

We note that our data are mostly dominated by the measurement uncertainties so that the choice of the prior for τ has no effect on the result, and both converge around the same point (Sect. 6.3). Partial pooling determines the two shape parameters of the ejecta mass distribution, Eq. (25). The upper bounds of μ can then be interpreted as conservative upper bounds on the ejecta masses of ⁷Be and ²²Na in individual CN events in the Milky Way, taking into account the complete population and its uncertainties from the SPI measurements.

6. Results

6.1. Individual objects

In Fig. 7 we illustrate the bounds on the ²²Na ejecta masses for each object at its (un)known distance, and how the estimates relate to upper bounds on the flux. All our upper bounds are quoted on the 99.85th percentile level if not mentioned otherwise. Apparently bright objects for which the flux is not consistent with zero within 2σ, the full posterior distributions are shown (here only V1535 Sco). Clearly, the mass estimates are inherently bound by the observation time towards the target and the exponential decay law winning over the square-root of time increase in sensitivity. Even though several objects have been frequently re-visited since the beginning of the INTEGRAL mission, all become saturated in terms of significance above background. This leads to an apparently common upper bound of the fluxes between 3 × 10⁻⁵ and 3 × 10⁻⁴ ph cm⁻² s⁻¹ for ²²Na and 10⁻⁴ to 3 × 10⁻³ ph cm⁻² s⁻¹ for ⁷Be (see Fig. A.1), which cannot be improved on the individual object basis. Only a nearby ONe nova (≲1.0 kpc) or one with exceptionally high ejecta mass would be seen with INTEGRAL/SPI. The remaining variance in flux bounds and hence ejecta mass bounds stems from the varying background level in this 16-year data set, as well as on exposure differences.

Fig. 7. Summary of the no pooling analysis. Shown are the upper bounds (99.85th perc.) on the ejecta mass of 22Na for each object with known (orange) and unknown (black) distance. The dependence on distance lets the bounds on mass appear along the lines of constant flux (dashed lines). Full posteriors are given for objects with the best flux bound (V5115 Sgr), the best mass bound (V5113 Sgr), and whenever the significance is higher than 2σ (V1535 Sco). Theoretical expectations (Starrfield et al. 2009) are indicated by the green bands for 1.25 and 1.35 M⊙ ONe novae. Other theoretical expectations are bound to values below ≈ 10−8 M⊙ (e.g., José & Hernanz 1998). See text for a discussion of individual objects.

Nova V1535 Sco shows the highest significance above instrumental background with 2.5σ. With a distance of 5.3 ± 1.7 kpc, this would be an unexpectedly high ²²Na ejecta mass. Furthermore, it is not known whether V1535 Sco is an ONe or CO nova. Hard X-ray emission (> 1 keV) has been found in V1535 Sco within a few weeks after its discovery on February 11, 2015, with a peculiar behaviour that is thought to be an indication of the CN being embedded in a red giant star (Linford et al. 2017). V1535 shows no to weak indications of GeV emission (Franckowiak et al. 2018), probably because of its large distance. We find this enhanced S/N only for the 1275 keV line, which would be unexpected if the emission was due to a CN outburst.

In our sample, we obtain the lowest upper bound on the 1275 keV line flux for V5115 Sgr. The large exposure toward the Galactic centre (cf. Fig. 5) and its outburst time on March 28, 2005 (Nakano et al. 2005), leads to an upper bound on the 1275 keV line flux of < 3.6 × 10⁻⁵ ph cm⁻² s⁻¹. V5115 Sgr is probably an ONe nova with a WD mass of 1.2 ± 0.1 M_⊙ (Hachisu & Kato 2006; Shara et al. 2018). At a distance of 2.8 ± 1.0 kpc, the upper bound on the flux converts into an upper bound on the ²²Na ejecta mass of < 1.3 × 10⁻⁷ M_⊙. This excludes WD masses around 1.3 M_⊙ (see Fig. 7), fully consistent with theoretical expectations (Starrfield et al. 2009) and independent measurements (Hachisu & Kato 2006). Using this bound, we can set a lower bound on the distance of V5115 Sgr of ≳2.0 kpc because otherwise, SPI would have seen the 1275 keV line. We note again that longer exposures for this object will not result in any improvements on the flux or ejecta mass because most ²²Na has already decayed.

We obtain the tightest limit on the ²²Na ejecta mass for the object closest to us, as expected, here V5113 Sgr with a distance estimate of 0.9  ±  0.2 kpc. At this distance, we infer a bound on the mass up to < 2 × 10⁻⁸ M_⊙. If V5113 Sgr was an ONe nova, this would set tight constraints on the nucleosynthesis on such objects, as it would fall below the lowest theoretical estimates (e.g., Starrfield et al. 2009). There are no indications that V5113 Sgr is indeed an ONe nova (Tanaka et al. 2011; Mróz et al. 2015), nor is the mass of the WD constrained.

We note that the mass estimates for ⁷Be and 478 keV are one to two orders of magnitude above theoretical expectations in most cases. We refer to Appendix A for these results.

6.2. Diffuse emission

In Fig. 8 we show the posterior distribution of the diffuse 1275 keV line flux in the Galaxy, accounting for possible continuum emission, and separated into ONe nova rate and ²²Na ejecta mass. The illustrated flux posterior takes into account that a weak diffuse Galactic continuum emission is present in the band 1265–1285 keV. We determined the total flux in this band according to $F_{\mathrm{tot}}^{22}=F_{\mathrm{Na}}^{22}+F_{\mathrm{conti}}^{22}$, where $F_{\mathrm{conti}}^{22}\sim \mathcal{N}(7.6,2.5^2)$ in units of 10⁻⁵ ph cm⁻² s⁻¹ (Wang et al. 2020). From this, we find an upper bound on the diffuse 1275 keV line flux of < 3.9 × 10⁻⁴ ph cm⁻² s⁻¹ in the entire Galaxy, which converts into an upper bound on the ²²Na ejecta mass of < 2.7 × 10⁻⁷ M_⊙. We wish to point out that this value is similar to the one found by Jean et al. (2001) with COMPTEL (< 3 × 10⁻⁷ M_⊙), which only focused on the inner Galactic ridge (|l|< 30°). Our analysis takes the uncertainties on the ONe nova rate, the asymmetric flux estimate from the SPI count data, as well as possible continuum contribution into account, whereas the analysis by Jean et al. (2001) considered different spatial distributions for CNe.

Fig. 8. Posterior distributions for diffuse emission of the 1275 keV line. The flux (top right) is separated into its two determining components, the 22Na ejecta mass (top left) and the ONe nova rate (bottom right, treated as a nuisance parameter). Theoretical predictions are indicated in green from Starrfield et al. (2009). The joint posterior (68.3th, 95.4th, and 99.7th percentile) of ejecta mass and nova rate (bottom left) shows the expected anticorrelated behaviour according to Eq. (17).

The energy band 474–482 keV is detected with more than 11σ above instrumental background with a flux of $F_{\text{tot}}^7=(1.2\pm 0.1)\times {10}^{-3}\text{ph}{\text{cm}}^{-2}{\text{s}}^{-1}$. Most of this flux originates in the three-photon decay of ortho-Positronium in the interstellar medium, plus small contributions from inverse Compton scattering and unresolved high-energy sources for an integrated continuum flux of $F_{\text{conti}}^7\sim \mathcal{N}(99.2,{11.7}^2)$ (e.g., Strong et al. 2005; Bouchet et al. 2011; Siegert et al. 2016). When we take this into account, we obtain an upper bound on the diffuse 478 keV line flux of < 5.9 × 10⁻⁴ ph cm⁻² s⁻¹, which converts into a ⁷Be ejecta mass of < 4.1 × 10⁻⁷ M_⊙. Again, the ⁷Be values are not constraining, and we refer to Appendix A for the results of the 478 keV line.

6.3. Combined analysis

Because none of the objects are individually detected, the hierarchical analysis mainly takes into account the asymmetric uncertainties on these flux measurements. This allows us to define common (population) upper bounds on the ejecta mass, given the individual uncertainties from fits to the SPI raw count data.

In Fig. 9 we show the summary of our hierarchical model for the ²²Na case. As τ goes to zero in the case of complete pooling, all values (here: upper bounds on the ejecta mass) converge to one single point, $M_{\mathrm{ej},\mathrm{CP}}^{22}<0.7\times {10}^{-7}{M}_{\odot }$. This assumes that the ejecta mass of all individual point sources in our CN sample and the diffuse emission to originate from the same ²²Na ejecta mass. We know that this is not the case (see Sect. 6.1), and we describe two ways of lifting these restrictions below.

Fig. 9. Summary of the hierarchical model. Shown are the upper bounds on the 22Na ejecta masses (left axis) for individual objects (black) as a function of the population width τ, together with the complete pooling estimate (τ → 0, dashed orange), and the population mean μ (solid orange). Individual objects of interest (best (blue) worst (aqua) bound, diffuse emission (red)) are indicated. As τ → ∞, each upper bound converges to its no pooling value (Table B.1). Given the SPI flux posteriors, a population width of $\tau =0.{52}_{0.23}^{+0.18}$ (right axis, grey shaded) is found. This value is dominated by the flux uncertainties.

In an attempt to account for the fraction of 1/3 of all CNe to be ONe novae, we used a finite mixture model with exactly this defined ratio R_CO : R_ONe = 2 : 1 and repeated the CP analysis. This mixture model selected a sub-sample of one-third of our sample and fitted it according to Eq. (24), whereas the remaining sample was subject to a prior distribution for the ejecta mass according to CO novae ($lg{M}_{22}^{\mathrm{ej},\mathrm{CO}}\sim \mathcal{N}(-12,1^2)$, Starrfield et al. 2020), resulting in basically zero flux for every object. The fit resulted in an increased upper bound of $M_{\mathrm{ej},\mathrm{CP}}^{22}<1.0\times {10}^{-7}{M}_{\odot }$.

Alternatively, we used the partial pooling approach (Sect. 5.2.3) and inferred the upper bound on the population ejecta mass by determining the spread τ of the distribution of ejecta masses. Ideally, τ would be determined by a combination of the true spread of the ejecta from many observations and their inherent uncertainties. Here, τ is almost purely defined by the measurement uncertainties of the asymmetric flux posteriors. We find that independent of the choice on the prior of τ (cf. Sect. 5.2.3), its posterior populates the range $\tau =0.{52}_{-0.23}^{+0.18}$ (C.I. 99.7%) with a population mean μ₂₂ < 2.0 × 10⁻⁷ M_⊙. This means that the broadest allowed ejecta mass distribution of the population of ²²Na-ejecting CNe that would still be consistent with our data is $lg{M}_{22}^{\mathrm{ej}}\sim \mathcal{N}(-6.7,0.7^2)$. In other words, if the mean of the ejecta mass distribution were any higher than 2 × 10⁻⁷ M_⊙, one or more objects (or the diffuse emission) would have been seen by SPI with at least 3σ above the instrumental background. This value is higher than the tightest constraint from V5113 Sgr in the no pooling context ($M_{22}^{\mathrm{ej}}(\text{V5113 Sgr})<1.8\times {10}^{-8}{M}_{\odot }$ as τ → ∞ and assuming V5113 Sgr is an ONe nova; cf. Table B.1), for example, but smaller than the one purely from diffuse emission. This can be considered the most conservative upper bound on the ²²Na ejecta mass from ONe novae because 1) it does not restrict high ejecta masses as would be case in the complete pooling setting, 2) it takes the diffuse emission and known explosions into account, and 3) it considers the entire Galaxy. Using a log-normal prior for the rate decreases the upper bound on the population mean by less than 5%.

While the rate has only a marginal effect on the final population estimate, the non-detection of individual sources also suppresses the diffuse emission flux and hence lowers the ONe rate: The initial prior mean was ${\mathrm{\mu }}_{R_{\mathrm{ONe}}}=16.7{\mathrm{yr}}^{-1}$ with a standard deviation of ${\sigma }_{R_{\mathrm{ONe}}}=8.3{\mathrm{yr}}^{-1}$. Now, the posterior mean of μ_RONe is found at 6 (complete pooling) and 13 yr⁻¹ (partial pooling), with a standard deviation of 3 and 7 yr⁻¹, respectively, independent of the prior shape. If all CNe happen to eject the same small amount of ²²Na as given in the complete pooling setting, the ONe nova rate in the Milky Way would be 18 yr⁻¹ at most (99.85th percentile). Clearly, this assumption is extreme and the more conservative upper bound on the ONe nova rate from partial pooling is 40 yr⁻¹, also independent of the prior shape.

In the case of ⁷Be and the 478 keV line, the shape of Fig. 9 is almost identical, only the numbers are spread further. We find an upper bound on the ⁷Be ejecta mass from complete pooling of $M_{\mathrm{ej},\mathrm{CP}}^7<0.5\times {10}^{-7}{M}_{\odot }$ and from partial pooling of μ₇ < 2.5 × 10⁻⁷ M_⊙. Neither of these values is constraining compared to theoretical expectations, which predict at most 4 × 10⁻¹⁰ M_⊙ (e.g., Starrfield et al. 2009). However, individual CNe have been found with a higher ⁷Be ejecta mass, of about several 10⁻⁹ M_⊙ (e.g., Molaro et al. 2016; Tajitsu et al. 2016, 2015), which indeed suggests a larger spread in ejecta masses than what is currently discussed.

7. Discussion

7.1. Nova nucleosynthesis

All our upper bounds are consistent with theoretical expectations. Considering the lowest theoretical expectations for 1.35 M_⊙ WDs from Starrfield et al. (2009), for example, of ≈ 3.9 × 10⁻⁷ M_⊙, we find 18 objects whose upper bounds on the ejecta mass would be below this threshold. These 18 objects include V5115 Sgr and V5113 Sgr, which we mentioned above (Sect. 6.1). We find that 9 objects are probably CO novae, 3 are probably ONe novae, and 6 are not determined. Considering the three ONe novae V5115 Sgr, V382 Nor, and V1187 Sco and our ejecta mass estimates ($M_{22}^{\mathrm{ej}}(\text{V5115 Sgr})<1.3\times {10}^{-7}{M}_{\odot }$, $M_{22}^{\mathrm{ej}}(\text{V382 Nor})<2.2\times {10}^{-7}{M}_{\odot }$, $M_{22}^{\mathrm{ej}}(\text{V1187 Sco})<1.3\times {10}^{-7}{M}_{\odot }$; cf. Table B.1), we find that these events occurred on WDs with masses lower than 1.35 M_⊙. This is fully consistent with the independent WD mass estimates from Shara et al. (2018). In general, our upper bounds can be interpreted to mean that most ONe novae in the Milky Way occur on WDs with masses lower than 1.35 M_⊙ because otherwise, SPI would have seen a stronger signal. This is expected because the mass distribution of measured WDs in the Galaxy peaks around 1.13 M_⊙ with a width of 0.12 M_⊙ (Shara et al. 2018). Other theoretical models predict lower ²²Na yields of about 10⁻⁸ M_⊙ (e.g., José & Hernanz 1998). When these models are considered, our results cannot provide any constraints.

The upper bound on the ²²Na ejecta mass from ONe novae translates into a steady-state ²²Na production of 3.3 × 10⁻⁶ M_⊙ yr⁻¹. On a timescale of 10 Gyr, ONe novae hence produced 3.3 × 10⁴ M_⊙ of ²²Ne at most, assuming the ONe nova rate is constant. With the solar abundance ratios by Lodders (2003), the Galaxy would contain a ²²Ne mass of ≈7.5 × 10⁴ M_⊙. Clearly, a full population synthesis model would be required to compare the upper bound on the ²²Ne mass from CNe with the total Galactic content. We nevertheless provide these order-of-magnitude evaluations to show a general broad consistency.

While the results from ⁷Be are far from constraining, we can provide a similar upper bound on the total ⁷Li mass inside the Milky Way. With the upper bound on the ⁷Be ejecta mass from the population of CNe of ${\mathrm{\mu }}_7^{\mathrm{ej}}<2.5\times {10}^{-7}{M}_{\odot }$, we find that on a timescale of 10 Gyr, 1.25 × 10⁵ M_⊙ of ⁷Li has been produced at most. This is about two to three orders of magnitude above the theoretical values from Starrfield et al. (2020), who suggested that about 100 M_⊙ of the 1000 M_⊙ of ⁷Li in the Milky Way is due to CO novae. To further constrain the ⁷Li production in the Milky Way from CNe with soft γ-rays, an instrument with a sensitivity improvement by two orders of magnitude would be required (see also next Section).

7.2. Positrons from classical novae

Given the upper bound on the ²²Na ejecta mass from ONe novae in the Milky Way of ${\mathrm{\mu }}_{22}^{\mathrm{ej}}<2.0\times {10}^{-7}{M}_{\odot }$, we can estimate how much CNe contribute at most to the Galactic positron budget. With the probability of ²²Na to decay through β⁺ decay of $p_{22}^{+}=0.904$, the production rate of positrons from the population of ONe novae is

$$\begin{array}{c}\hfill {\dot{N}}_{e^{+}}=R_{\mathrm{ONe}}{p}_{22}^{+}\frac{{\mathrm{\mu }}_{22}^{\mathrm{ej}}}{m_{22}}<5.5\times {10}^{42}{\mathrm{e}}^{+}{\mathrm{s}}^{-1}.\end{array}$$(27)

Siegert et al. (2016) estimated a total positron annihilation rate in the Milky Way of R_e⁺ ≈ 5 × 10⁴³ e⁺ s⁻¹. Assuming a steady-state annihilation in the Galaxy, this rate has to be sustained by the same production rate. The population of ONe would therefore contribute 10% at most to the total Galactic positron production rate. This is of similar magnitude as other known and validated source types, such as ²⁶Al from massive stars and their supernovae (≈10%; τ₂₆ = 1.05 Myr) and ⁴⁴Ti from core-collapse supernovae (≈10%; τ₄₄ = 86 yr). For the remaining 80% required for steady-state annihilation, several other individual source types, populations, and scenarios have been suggested. However, only indirect evidence for positron annihilation and production in these sources has been found so far (cf. Prantzos et al. 2011, for a review).

The population of CNe in the Galaxy will have a small but important contribution to the total positron production rate, complementing the different nucleosynthesis origins of Galactic positrons with frequent occurrences of CNe and comparably short lifetimes of parent nuclei. These considerations do not exclude the very short-lived light β⁺-unstable nuclei ¹³N or ¹⁸F from contributing significantly to the total content. Because the expected very strong 511 keV flash at the time of the explosion has never been observed, it might even be possible that these positrons also do not annihilate in the expanding nova cloud, but rather escape into the interstellar medium. This would make CNe a dominant source of positrons in the Milky Way.

Theoretical modelling of nucleosynthesis yields predicts values close to our upper bounds. This suggests a bright future for the next generation of soft γ-ray telescopes such as COSI (Tomsick et al. 2019) in nuclear decay and the 511 keV positron annihilation line.

8. Conclusion

We used 15 years of archival INTEGRAL/SPI γ-ray observations to infer nucleosynthesis ejecta masses of decaying ²²Na and ⁷Be from the population of classical novae. In our sample of 97 individual known objects, no signal is detected above a threshold of 3σ.

The best upper bound for an individual CO nova is found for V5668 Sgr owing to a dedicated observation campaign (see Siegert et al. 2018), with a 478 keV line flux of < 1.4 × 10⁻⁴ ph cm⁻² s⁻¹, resulting in a ⁷Be mass of < 3 × 10⁻⁸ M_⊙. For the first time, we considered the diffuse 478 keV line component of CO novae in the Milky Way and found an upper flux bound of < 6.0 × 10⁻⁴ ph cm⁻² s⁻¹, which converts into an average ejecta mass of < 4.1 × 10⁻⁷ M_⊙. We improved these values by using a Bayesian hierarchical model of 478 keV emitting novae, that is, we took a common distribution of ejecta masses into account and inferred a population mean of < 0.5– 2.5 × 10⁻⁷ M_⊙. While these values are hardly constraining compared to theoretical expectations (e.g., Starrfield et al. 2020), we nevertheless showed that proper modelling of sub-threshold sources can provide a significant improvement on the pure stacking of many point sources.

For individual ONe novae, our most constraining upper bounds on the ²²Na ejecta mass are found for V5115 Sgr and V1187 Sco with < 1.3 × 10⁻⁷ M_⊙. The lowest 1275 keV line upper bound is found for V2659 Cyg with < 2.9 × 10⁻⁵ ph cm⁻² s⁻¹, but at an unknown distance. The expected diffuse emission of the 1275 keV line in the Milky Way along the Galactic plane with a strong peak in the bulge provides an upper flux bound of < 3.9 × 10⁻⁴ ph cm⁻² s⁻¹, which converts into a ²²Na ejecta mass bound of < 2.7 × 10⁻⁷ M_⊙. Our Bayesian hierarchical model for the 1275 keV line results in a population mean of < 0.7–2.0 × 10⁻⁷ M_⊙. These values are an improvement over previous studies: (1) the upper bounds on the ejecta masses are 2–5 (depending on the model assumptions) times smaller than from earlier COMPTEL measurements (Jean et al. 2001), and (2) they take into account both diffuse emission and known sources, their individual uncertainties, and uncertainties on circumstantial parameters, such as the nova rate and distances.

The result can be interpreted in a way that if the population mean of ejecta masses were any higher, one or more ONe novae would have been detected with at least 3σ by INTEGRAL/SPI in this 15 yr data set. Still, these upper bounds are on the high side of theoretical expectations (e.g., Starrfield et al. 2009): We can exclude that most ONe novae occur on WDs with masses of about 1.35 M_⊙, and thus finally probe the region of model calculations. When we tighten the uncertainties on nova rate and individual distances, for example, and at the same time might make use of the full data set rather than relying on extracted fluxes, the hierarchical model would probably probe even deeper, that is, set tighter bounds on ejecta masses.

With ²²Na as a β⁺ decayer, we find an upper bound on the positron production rate in the Milky Way from ONe novae of < 5.5 × 10⁴² e⁺ s⁻¹. Compared to the Galactic positron annihilation rate (Siegert et al. 2016), the population of ONe could make up to 10% of the required positron production rate for a steady-state configuration.

Acknowledgments

Thomas Siegert is supported by the German Research Foundation (DFG-Forschungsstipendium SI 2502/1-1 & SI 2502/3-1). He thanks J. Michael Burgess for suggesting to do something important. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1548562. We acknowledge the use of XSEDE resources via Startup allocation TG-PHY200045. Software: OSA/spimodfit (Halloin 2009), numpy (Oliphant 2006), matplotlib (Hunter 2007), astropy (Astropy Collaboration 2013), scipy (Virtanen et al. 2020), Stan/pystan (Carpenter et al. 2017), arviz (Kumar et al. 2019).

References

Appendix A: Additional figures for ⁷Be

The results for ⁷Be from the measurements of the 478 keV line are less constraining than for ²²Na and the 1275 keV line. For completeness, we show the resulting upper bounds for ⁷Be ejecta from individual objects in Fig. A.1 together with theoretical expectations (Starrfield et al. 2020) and one measured value for V5668 Sgr from Molaro et al. (2016, M+16). In Fig. A.2 we show the joint posterior of the diffuse component of the 478 keV line in the Milky Way, in which the strong continuum contribution is already marginalised out. Finally, Fig. A.3 shows the summary of our hierarchical analysis in the case of ⁷Be, with both extremes of complete pooling (τ → 0) and no pooling (τ → ∞). The posterior distribution of τ is completely determined by the uncertainties of the individual flux measurements.

Fig. A.1. Same as Fig. 7 for the case of 7Be and the 478 keV line.

Fig. A.2. Same as Fig. 8 for the 478 keV line.

Fig. A.3. Same as Fig. 9 for 7Be. The x-axis has been extended to show the convergence to the no pooling estimate for large τ.

Appendix B: Complete results table

Our complete results are presented in Table B.1.

All Tables

All Figures

	Fig. 1. Distance distribution of CNe according to Eqs. (9) and (10) through rejection sampling. Shown are 104 samples (grey histogram, corresponding to a timescale of 200 yr, left axis), and an approximation with the Γ-distribution in red. The cumulative distributions are shown in black for the sample and in cyan for the Γ-distribution (right axis).
In the text

	Fig. 2. Diffuse emission template from line-of-sight-integrated nova density distribution, Eq. (15), together with nova sample (cyan points). The contours indicate the De Vaucouleurs profile for the bulge and the exponential disc.
In the text

Fig. 3. Population synthesis of nova events in a Milky Way-like galaxy. Left: year 116.3 of the modelled Poisson process, where a nearby CO nova at (l, b) = (+123° , + 8° ) outshines the remaining galaxy in 478 keV emission (top). At the same time, 1275 keV emission is not enhanced (bottom). The angular resolution is chosen according to the SPI 2.7° (FWHM). Right: time profiles (light curves) of the 478 keV and 1275 keV since the beginning of the synthesis for a duration of ∼150 yr. This timescale is enough to reach convergence and to show characteristic features. The vertical purple line indicates the time shown in the plots on the left.

In the text

Fig. 4. Estimated S/N given the SPI data set (Sect. 3) and the range of plausible 22Na ejecta masses for our sample of 97 novae. Top: the S/N of individual objects does not increase with $\sqrt{T_{\exp }}$ because the exponential decay is faster than the accumulation of 1275 keV photons (cyan lines). The cumulative signal of point sources saturates on a timescale of a few years (dark blue). If all objects are assumed to be ONe novae, the S/N would be about three times higher (pale red). Bottom: comparison between cumulative point-source emission significance and diffuse-emission significance during the INTEGRAL mission timescale. The bands contain the 68th percentile from 1000 realisations of the population synthesis model (boundaries marked by dashed lines, median by solid lines; Sect. 2.4), to distribute ejecta masses among sources and diffuse emission.

In the text

	Fig. 5. Exposure map of the selected data set. The effective areas around 478 and 1275 keV are ∼80 and ∼55 cm^2, respectively.
In the text

Fig. 6. Complete graph of different model variants used in this study. When the source models are used individually, that is, regardless of a common underlying physical process, the most conservative estimates are found in a no pooling setting (Sect. 5.2.1). The opposite extreme assumes that all objects eject the same amount of matter, so that instead of ∼100 different masses, only one is fitted. This resembles complete pooling and defines tight constraints but is subject to bias if the true variance is large (Sect. 5.2.2). A compromise between no and complete pooling is achieved by invoking the hyper-parameters μ and τ that determine the mean and the spread of ejecta masses, if any (partial pooling, dash-dotted box, Sect. 5.2.3).

In the text

Fig. 7. Summary of the no pooling analysis. Shown are the upper bounds (99.85th perc.) on the ejecta mass of 22Na for each object with known (orange) and unknown (black) distance. The dependence on distance lets the bounds on mass appear along the lines of constant flux (dashed lines). Full posteriors are given for objects with the best flux bound (V5115 Sgr), the best mass bound (V5113 Sgr), and whenever the significance is higher than 2σ (V1535 Sco). Theoretical expectations (Starrfield et al. 2009) are indicated by the green bands for 1.25 and 1.35 M⊙ ONe novae. Other theoretical expectations are bound to values below ≈ 10−8 M⊙ (e.g., José & Hernanz 1998). See text for a discussion of individual objects.

In the text

Fig. 8. Posterior distributions for diffuse emission of the 1275 keV line. The flux (top right) is separated into its two determining components, the 22Na ejecta mass (top left) and the ONe nova rate (bottom right, treated as a nuisance parameter). Theoretical predictions are indicated in green from Starrfield et al. (2009). The joint posterior (68.3th, 95.4th, and 99.7th percentile) of ejecta mass and nova rate (bottom left) shows the expected anticorrelated behaviour according to Eq. (17).

In the text

Fig. 9. Summary of the hierarchical model. Shown are the upper bounds on the 22Na ejecta masses (left axis) for individual objects (black) as a function of the population width τ, together with the complete pooling estimate (τ → 0, dashed orange), and the population mean μ (solid orange). Individual objects of interest (best (blue) worst (aqua) bound, diffuse emission (red)) are indicated. As τ → ∞, each upper bound converges to its no pooling value (Table B.1). Given the SPI flux posteriors, a population width of $\tau =0.{52}_{0.23}^{+0.18}$ (right axis, grey shaded) is found. This value is dominated by the flux uncertainties.

In the text

	Fig. A.1. Same as Fig. 7 for the case of 7Be and the 478 keV line.
In the text

	Fig. A.2. Same as Fig. 8 for the 478 keV line.
In the text

	Fig. A.3. Same as Fig. 9 for 7Be. The x-axis has been extended to show the convergence to the no pooling estimate for large τ.
In the text