Novae: An Important Source of Lithium in the Galaxy

Abstract

The source of the Galactic Lithium (Li) has long been a puzzle. With the discovery of Li in novae, extensive research has been conducted. However, there still exists a significant disparity between the observed abundance of lithium in novae and the existing theoretical predictions. Using the Modules for Experiments in Stellar Astrophysics (MESA), we simulate the evolution of nova with element diffusion and appropriately increased the amount of ³He in the mixtures. Element diffusion enhances the transport efficiency between the nuclear reaction zone and the convective region on the surface of the white dwarf during nova eruptions, which results in more ⁷Be to be transmitted to the white dwarf surface and ultimately ejected. Compared to the previous predictions, the abundance of ⁷Be in novae simulated in our model significantly increases. And the result is able to explain almost all observed novae. Using the method of population synthesis, we calculate Li yield in the Galaxy. We find that the Galactic occurrence rate of nova is about 130 yr^-1, and about 110 M_⊙ Li produced by nova eruption is ejected into the interstellar medium (ISM). About 73% of Li in the Galactic ISM originates from novae, and approximately 15%-20% of the entire Galaxy. It means that novae are the important source of Li in the Galactic.

1 Introduction

Lithium (Li) is an extremely fragile element that is destroyed in the hydrogen (H)-capture reaction at temperatures as low as 2 $\times$ $10^{6}$K. Its dominant isotope, ⁷Li, is the decay product of ⁷Be with a half-life of 53.3 days. In the process of H burning, ⁷Be rapidly decays to ⁷Li through electron capture (pp II) after being formed via the proton-proton (pp) chain. ⁷Li itself is also destroyed through proton capture, resulting in very little lithium surviving in the H burning process. This makes it almost impossible to accurately calculate the abundance of Li in most stars during their formation. The production of primordial ⁷Li is a sensitive function of the baryon-to-photon ratio and can be estimated within the framework of standard primordial nucleosynthesis as long as the baryon density is obtained from the initial deuterium abundance or fluctuations of the cosmic microwave background (CMB) (Spergel et al., 2003; Bennett et al., 2013; Planck Collaboration et al., 2020). The expected primordial value is A(Li)$\approx$2.72 dex (Życzkowski et al., 1998; Coc et al., 2014), which is 3 - 4 times higher than the measured values in halo dwarf stars (Spite & Spite, 1982; Planck Collaboration et al., 2020). This difference is often referred to as the Cosmological lithium problem (Fields et al., 2014).

Since the measured Li abundance in meteorites that preserves the protosolar in the interstellar medium (ISM) is A(Li)$\approx$3.3 dex (Asplund et al., 2009; Lodders et al., 2009), there is need for a galactic source to explain the increase from the initial value of 2.72 dex. This identification of the source is commonly known as the Galactic lithium problem. The ISM is the material that fills the space between the stars within a galaxy. It consists primarily of gas (atomic, molecular, and ionized) and dust. Planetary nebulae, supernova explosions, and novae explosions all inject a large amount of chemical elements into ISM. The ISM is crucial for the formation and evolution of stars, and it plays a key role in the chemical enrichment and energy balance of galaxies (Draine, 2011). Various nucleosynthesis processes and sources have been proposed so far. One confirmed source of Li is spallation and fusion processes of galactic cosmic rays (GCRs) in the ISM. The integrated spallation process is estimated to contribute about 10%$\sim$20% to the measured ⁷Li in the entire Galactic lifetime (Prantzos et al., 1993; Romano et al., 2001; Prantzos, 2012). Other proposed sources include spallation processes in the flares of low-mass active stars, red giants (RGs), asymptotic giant branch (AGB) stars and neutrino-induced nucleosynthesis during a type II supernova (Starrfield et al., 1978; Spite & Spite, 1982; Smith & Lambert, 1989; Smith, Verne V. and Lambert, David L., 1990; Hernanz et al., 1996; Romano et al., 1999; Travaglio et al., 2001; Alibés et al., 2002; Prantzos, 2012; Tajitsu et al., 2016; Banerjee, Projjwal et al., 2016; Pignatari et al., 2016; Rukeya et al., 2017). Although Li abundance is indeed enhanced in some studies, but their contribution to the entire Galaxy is too small. Therefore, other sources are still needed for the Galaxy to reach its current value. Subsequently, Romano et al. (2001) and Prantzos, N. et al. (2017) found that low-mass giants are the best candidates for reproducing the late rise of the Li metallicity plateau. However, high A(Li) cannot be sustained for a long period due to convective activity in these stars, and the low percentage of Li-rich giants also indicates this (Casey et al., 2016; Yan et al., 2018; Gao et al., 2022). Therefore, their contribution to the enrichment of Li in the Galactic ISM remains quite uncertain.

For many years, classical novae have been proposed as feasible sites for Li production (Arnould & Norgaard, 1975; Starrfield et al., 1978). The novae eruption is the result of unstable hydrogen-burning on the CO or ONeMg white dwarfs (WDs) surfaces, which accrete hydrogen-rich material from their main sequence (MS) or red giant (RG) phase companions in low-mass, close binary systems (Starrfield et al., 2009; Jose, 2016; Starrfield et al., 2016; Rukeya et al., 2017). When the companion star fills its Roche Lobe, the hydrogen-rich material flows through the inner Lagrange point to the WD and accretes onto its surface. This material accumulates and is compressed until thermonuclear runaways (TNR) is triggered, resulting in mass ejection that ultimately pollutes the interstellar environment. Both theoretical and observational evidence confirms that novae contribute many nucleosynthetic isotopes to the Galactic ISM (José & Hernanz, 1998; Starrfield et al., 1998). For example, the elements ⁷Li and ⁷Be.

⁷Be is generated through the ³He($\alpha$,$\gamma$)⁷Be reaction at a temperature of 150 million K (Hernanz et al., 1996). To avoid destruction, ⁷Be needs to be transported to cooler regions quickly through convection, as described in the Cameron–Fowler (CF) mechanism (Cameron & Fowler, 1971). When these cooler regions are subsequently ejected, the absorption lines of ⁷Be during nova outbursts can be observed (José & Hernanz, 1998), eventually decaying into ⁷Li. This view was later confirmed by observations by Izzo et al. (2015) who detected a potentially detectable ⁷Li I $\lambda$6708 Å absorption line in the spectrum of the Nova Centauri 2013 (V1369 Cen), providing observational evidence for the presence of Li in nova explosions that had been predicted since the mid-1970s but only recently discovered. Subsequently, the observation studies continuously detected the progenitor nucleus ⁷Be or ⁷Li in the prominent post-outburst spectra of classical novae (Tajitsu et al., 2015, 2016; Molaro et al., 2016; Izzo et al., 2018; Selvelli et al., 2018; Molaro et al., 2020, 2022, 2023).

For a long time, people have continuously used theoretical models to calculate the accurate Li yield from nova explosions. Hernanz et al. (1996) and José & Hernanz (1998) calculated the theoretical values for different WD masses and mixing ratios of ejected material in their nova models. Then Rukeya et al. (2017) expanded on the nova grid model and computed more detailed mass and Li yields of the ejecta. They concluded that the Li produced in novae accounts for 10% of the total Galactic ISM ($\sim$ 150 M_⊙) production. However, Starrfield et al. (2024) predicted that the abundance of Li ejected from CO novae explosions is A(Li) $\approx$ 6, but most observed novae showed A(Li) $\textgreater$ 7 (i.e. A(Li)=Log(N(⁷Li)/N(H)+12)). The actual abundance of ejected Li is one order of magnitude higher than the theoretical predictions (Izzo et al., 2015; Molaro et al., 2016; Izzo et al., 2018; Arai et al., 2021; Molaro et al., 2022).

This suggests that traditional nova models may have deficiencies in their physical mechanisms, and the Galactic lithium problem unresolved. Subsequently, researchers have attempted different calculation methods by adjusting the mass distribution of binary systems, eruption time intervals, criteria for eruptive episodes, metallicity, and accretion material mixing ratios to compute Li yields under different scenarios (José & Hernanz, 1998; Starrfield et al., 2009; Cescutti & Molaro, 2019; José et al., 2020; Starrfield et al., 2020). Kemp et al. (2022a) demonstrated that using the yield values given in Molaro et al. (2023), then novae can definitely be the main factories of lithium in the Galaxy while using the theoretical values, they could not. Therefore, finding a reasonable nova model and Li production mechanism that aligns with the observations and theoretical values is of great significance in addressing the Galactic lithium problem and even the Cosmological lithium problem. It is well known that the chemical composition of the ejecta in nova explosions depends on the surface chemical element composition of the WDs at the moment of the eruption. Several factors can influence the chemical composition of the WD’s surface, such as metallicity, accretion efficiency, and element diffusion (Kippenhahn et al., 1980; Dupuis et al., 1992; Zhu et al., 2021). The effects of nuclear reaction rate, accretion efficiency and metallicity on the surface material of WDs have been explored by Starrfield et al. (2016, 2020) and Kemp et al. (2022a, b). However, element diffusion has been rarely considered. Element diffusion is a dynamic process that alters the distribution of chemical elements within a star. It primarily results from the combined effects of pressure, temperature, and material concentration. Kovetz & Prialnik (1985) investigated that effective diffusion between the accreted layer and the inner regions for accreting WD, and they found that the diffusion can lead to enrichment of CNO (and other metals) in the ejected envelop. However, CNO enhancements are possible in CO novae only in the absence of the He buffer (Iben et al., 1991). The presence of helium layer would prevent the diffusion of CNO and metals (such as ³He), although the average abundance of helium is almost the same as that of solar (Gehrz et al., 1998; Yaron et al., 2005). Therefore, the ³He inside the WD cannot be effectively brought to the envelope where TNR occurs by element diffusion, and there are other sources of ³He in the envelope (Iben et al., 1991). In our model, ³He comes from the accreted material from the donor star but not the from the interior of WD. The role of element diffusion in our model is that ⁷Be produced by ³He($\alpha$,$\gamma$)⁷Be among TNR can be taken more efficiently to WD surface, and then can ejected more easily. Therefore, a large amount of ⁷Be can be detected in the ejecta of the nova, ultimately decaying into ⁷Li.

In this paper, we assume the novae as main sources of Li. In particular, Section 2 describes the details of the nova model. Section 3 presents the ⁷Li yields predicted by our models and the comparative analysis of the observation and theoretical results. The summary is given in Section 4.

2 Models

Our nova model has been using version-12778 of the MESA stellar evolution code, which constructs CO and ONeMg nova models based on its white dwarf and nova modules (Paxton et al., 2011, 2013, 2015, 2018, 2019). The advantage of MESA is its ability to calculate nova outbursts by dividing the outer shell mass of WD into approximately 1000 or more cells, resulting in smaller errors. At each cell, such as temperature, density and isotopic composition are calculated. In MESA, when the luminosity of a star ($L$) exceeds the super-Eddington luminosity ($L_{\rm Edd}$), it will trigger mass loss. The mass loss rate is

$$\dot{M}=-2\eta_{\rm Edd}\frac{(L-L_{\rm Edd})}{v_{\rm esc}^{2}}$$ (1)

where $v_{\rm esc}=\sqrt{2GM/R}$, $L_{\rm Edd}=(4\pi GcM)/\kappa$. $M$ and $R$ are the mass and radius, while $\kappa$ is the Rosseland mean opacity at the WD’s surface. The scaling factor is taken $\eta_{\rm Edd}=1$ (Denissenkov et al., 2013). The cells will be ejected when $L>L_{\rm Edd}$. In this paper, it is set that the nova ejection begins when the total WD luminosity ($L$) is greater than $10^{4}$ times the solar luminosity ($\rm L_{\odot}$), and ends when less than $10^{3}$ $\rm L_{\odot}$. Therefore, the MESA code can be used to simulate multi-cycle nova and construct a large-scale nova model grid from the accretion phase to expansion, explosion and ejection phases (Paxton et al., 2011, 2013, 2015).

For CO WD models, the nuclear network selects pp_and_cno_extras_net, while the ONeMg WD models uses h_burn_net. These nuclear networks include the CNO burning cycle and the proton-proton reaction chains (pp chain), the latter of which include ³He($\alpha,\gamma$)⁷Be, ⁷Be($e^{-},\nu$)⁷Li, ⁸B($\gamma,p$)⁷Be, ⁷Li($p,\alpha$)⁴He, and ⁸B($e^{+},\nu$)⁸Be^∗(2$\alpha$), which are sufficient nuclear synthesis to treat ⁷Be and ⁷Li nucleosynthesis (Starrfield et al., 1978). It is worth noting that the degree of mixing also affects the mass and isotopic abundance of the ejected material (Hernanz et al., 1996). José & Hernanz (1998), Denissenkov et al. (2014) and Rukeya et al. (2017) have used the NOVA and MESA codes to calculate the nova eruption model, assuming that the degree of mixing could be 25%-50%. That is, 25% of WD material and 75% of solar material, or 50% of WD material and 50% of solar material (Lodders et al., 2009). We adopted the latter, which is more widely used. However, the ⁷Be simulated by the previous models did not explain the observed values. Subsequent studies suggested that there is a corresponding relationship between the ⁷Be and ³He (Boffin et al., 1993; Hernanz et al., 1996; Molaro et al., 2020; Denissenkov et al., 2021).

The importance of ³He for novae was first studied by Schatzman (1951) in the context of a theory of novae powered by thermonuclear detonations. The initial ³He abundance for the accreted matter, could potentially come from the donor star (Shara, 1980; Townsley & Bildsten, 2004). As the donor star are ascending the red giant branch the convection dredges up ³He enriched material to the surface which is later expelled into the ISM by wind or is accreted by WDs in nova systems (Shen & Bildsten, 2009; Denissenkov et al., 2021). Dantona & Mazzitelli (1982) and Iben & Tutukov (1984) found that the mass fraction of accreted $X_{0}$(³He) can reach values as high as $4\times 10^{-3}$ during the evolution of systems with low-mass donors. Furthermore, some observational studies have shown the existence of relatively high levels of ³He/H in certain planetary nebulae (Rood et al., 1992; Balser & Bania, 2018). If ³He abundance of the donor’s envelope is much higher than that in the Sun, this would almost align the model predictions of ⁷Be production with observations (Molaro et al., 2020). However, observations of the ISM indicate that the abundance of ³He in the Galaxy cannot be too high (Dearborn et al., 1996; Romano & Matteucci, 2003). Denissenkov et al. (2021) demonstrated that an excess of ³He in accreted material would actually reduce the production of ⁷Be. Their work showed that when $X_{0}$(³He) exceeds $3\times 10^{4}$, it could lead to the early onset of the TNR, result in a reduction in peak temperature and accreted mass, and thereby suppress the production of ⁷Be. While this level of ³He-induced ⁷Be abundance increase would be slightly elevated, most of it would still be consumed at high temperatures and not transported to the WD surface to generate sufficient ⁷Li, thus not closing the gap with observations. Therefore, an effective transport mechanism is required to safely transport ⁷Be to the WD surface and allow it to survive. Element diffusion provides an effective pathway for this.

Elemental diffusion in stellar interior is mainly driven by a combination of pressure gradients (or gravity), temperature gradients, compositional gradients, and radiation pressures. By solving the Burgers equation (Burgers, 1969), Thoul et al. (1994) proposed a general method of arranging the entire system of equations into a single matrix equation, so that the relative concentrations of various species have no approximate values and the number of elements considered is not limited. Therefore, this method is applicable to various astrophysical problems. In MESA, the method of Thoul et al. (1994) is used to calculate the diffusion of chemical elements within stars (Paxton et al., 2015, 2018). This diffusion effect can bring internal elements to the surface. We speculate that when WD accretes enough hydrogen-rich material from companion star to reach the critical point, TNR occurs. The thermal instability caused by the thermonuclear runaway leads to the appearance of a convective zone, which extends from approximately the middle of the combustion shell to the surrounding non combustion helium layer (Cameron & Fowler, 1971). The convective zone can bring ⁷Be from the stellar interior to the surface, which is the CF mechanism. Element diffusion improves the mixing efficiency of the entire convection zone, bringing more ⁷Be to the surface. When a nova explosion occurs, a large amount of ⁷Be on the surface is ejected and decays into ⁷Li in the ISM through a short half-life. This work utilize MESA to construct a large-scale nova grid by setting different basic parameters and calculated the ⁷Li production of the nova.

MESA is one-dimensional code, and only simulate spherically symmetric nova outburst. However, based on the images of NOVA V5668 SAG combining the optical and radio observations from the Hubble Space Telescope and the Very Large Array telescope, Mukai & Sokoloski (2019) suggested that the geometry distribution of nova ejecta is not spherical symmetry but is the shape of an equatorial torus. It indicates that ejecta from novae has very complicated structure. It is unclear which physical process led to such an unsymmetric structure. Up until now, there have been few instances of 2D or 3D simulations depicting nova outbursts. Therefore, it should be noted that the MESA nova model also has its limitations.

3 result

According to Yaron et al. (2005), the WD mass ($M_{\rm WD}$), accretion rate ($\dot{M}$), WD core temperature ($T_{\rm c}$), and composition of the accreted material constitute the four basic input parameters that determine the nova eruption (Starrfield et al., 2016; Jose, 2016). Yaron et al. (2005) constrained the parameters of the nova eruption and set a three-dimensional restricted space for the nova eruption conditions. Our model set the $T_{\rm c}$ at 3$\times$10⁷K and selected CO WD model mass: 0.5 $\rm M_{\odot}$, 0.6 $\rm M_{\odot}$, 0.7 $\rm M_{\odot}$, 0.8 $\rm M_{\odot}$, 0.9 $\rm M_{\odot}$, 1.0 $\rm M_{\odot}$, 1.1 $\rm M_{\odot}$, 1.2 $\rm M_{\odot}$; and ONeMg WD model mass: 1.0 $\rm M_{\odot}$, 1.1 $\rm M_{\odot}$, 1.2 $\rm M_{\odot}$, 1.3 $\rm M_{\odot}$; with accretion rates ($\dot{M}$) of 1 $\times$ $10^{-11}$ $\rm M_{\odot}$ yr^-1, 1 $\times$ $10^{-10}$ $\rm M_{\odot}$ yr^-1, 1 $\times$ $10^{-9}$ $\rm M_{\odot}$ yr^-1, 1 $\times$ $10^{-8}$ $\rm M_{\odot}$ yr^-1, 1 $\times$ $10^{-7}$ $\rm M_{\odot}$ yr^-1. In addition to the above four parameters, it is expected that material transferred from the companion (solar-like) will mix with the outer layers of the WD. This model adopts the widely accepted composition of 50% WD material and 50% solar material. Solar component data is taken from Lodders et al. (2009). Molaro et al. (2020) argued that the variety of high ⁷Be/H abundances in nova could be originated in a higher than solar content of ³He in the donor star. Observations show relatively high levels of ³He/H in some planetary nebulae, but observations of ISM indicate that the abundance of ³He in the Galaxy cannot be too high (Rood et al., 1992; Dearborn et al., 1996; Romano & Matteucci, 2003; Balser & Bania, 2018). According to the $X_{0}$(³He) range: 2.96$\times$10^-5$\thicksim$2.96$\times$10^-3 used by Denissenkov et al. (2021), we set $X_{0}$(³He) to be 4$\times$10⁴.

3.1 Novae models as stellar sources of Lithium

Figure 1: Nova explosion process with $M_{\rm WD}$ = 1$\rm M_{\odot}$ and $\dot{M}$ = 1 $\times$ $10^{-10}$ $\rm M_{\odot}$ yr^-1. Profiles of ⁷Be during the nova explosion show in panels. The y axis represents the boundary radius of WD during explosion and the x axis represents the lifetime of the nova explosion. Green shadow shows the convective zone of different models.

It is well known that during the occurrence of TNR, there are intense nuclear reactions, primarily dominated by hydrogen burning due to the WD accreting hydrogen-rich material from the companion. During hydrogen burning, a large amount of ⁷Be is produced through ³He($\alpha$,$\gamma$)⁷Be, with some being convectively transported to the WD’s surface while the rest is depleted at high temperatures. In our models with element diffusion, the transport efficiency within the WD is enhanced. Figure 1 shows the evolution of elements during nova explosion. The left two panels clearly display the effect of element diffusion. There are some blank discontinuous area on the surface convection zone (green shadow) of the model without element diffusion, indicating weak or even no mixing activity in these areas, as shown in the left-upper panel. In Figure 1, the mixing area within the convection zone of the WD model with element diffusion becomes more continuous and dense, and the mixing activity becomes stronger, as shown in the left-lower panel. In the nova model with element diffusion, a large amount of ⁷Be generated by nuclear reactions is effectively transported by this strong mixing effect to the low-temperature region of the surface, thereby avoiding depletion at high temperatures. On the contrary, the model without element diffusion cannot send out more ⁷Be, which is consumed by high temperature.

Meanwhile, the improvement of ³He can generate more ⁷Be to contribute to the transport to the surface as shown by the right-top panel in Figure 1. The ³He increase greatly enhances the yield of ⁷Be in envelope, which is shown by pink color area. However, the increase of ³He causes the convective region to move inward, resulting in the failure to bring ⁷Be to a sufficiently low temperature surface. This result is consistent with Denissenkov et al. (2021), who stated that excessively high ³He actually reduces the ⁷Be of nova ejection. However, combining the element diffusion and the increase of ³He, the large amount of ⁷Be produced via enhancing ³He can be efficiently transfered to WD surface, which are clearly shown by the right-lower panel.

Figure 2: Comparison with observed Li abundances and the numerical simulations along with WD masses. Solid lines of different colors represent CO nova models with different accretion rates, while dashed lines represent ONeMg nova models. Red bordered circles denote the observed abundances of CO WD novae, and the squares denote the ONeMg WD novae. Since evaluating the quality of each WD is difficult, we simply plotted the error bars for each observed nova. The black, yellow, and red dotted lines correspond to theoretical predictions; highest values for each WD mass in José & Hernanz (1998), the “123-321 models” of José et al. (2020), and 25%–50% mixing model of Starrfield et al. (2024), respectively. The blue pentagram represents the prediction of nova M31N 1926-07c and the blue triangle represents the nova T CrB.

We find that the abundance of ⁷Be in the model combining the element diffusion and the increase of ³He can be increased by about 4 - 10 times. It means that ⁷Be abundance calculated in this work is consistent with the observational values (Izzo et al., 2015; Tajitsu et al., 2015; Molaro et al., 2020, 2023).

In order to compare our results with the observations, we selected 13 published observational samples of CO and ONeMg novae with relatively accurate measurements of parameters such as mass and Be abundance (9 CO novae and 4 ONeMg novae) (e.g. Tajitsu et al. (2015, 2016); Molaro et al. (2016); Izzo et al. (2018); Selvelli et al. (2018); Molaro et al. (2020); Arai et al. (2021); Molaro et al. (2022, 2023)). We calculated the ⁷Be abundance ejected in a single nova eruption using the criteria that the eruption begins when the total WD luminosity ($L$) is greater than $10^{4}$ $\rm L_{\odot}$ and ends when it is below $10^{3}$ $\rm L_{\odot}$.

Figure 2 presents a comparison of the Be abundance produced by our simulated nova models with other models and matching to observations. The solid lines in different colors represent the CO nova models in the mass range of $0.7\rm M_{\odot}-1.0\rm M_{\odot}$, with accretion rates from 1 $\times$ $10^{-7}$ $\rm M_{\odot}$ yr^-1 to 1 $\times$ $10^{-11}$ $\rm M_{\odot}$ yr^-1. The bold dashed lines depict the results of the ONeMg nova models in the mass range of $1.1\rm M_{\odot}-1.3\rm M_{\odot}$, with accretion rates from 1 $\times$ $10^{-7}$ $\rm M_{\odot}$ yr^-1 to 1 $\times$ $10^{-9}$ $\rm M_{\odot}$ yr^-1. The results simulated by other literatures also are given, such as the black and red dotted lines, representing models constructed by José & Hernanz (1998); José et al. (2020) with 50% mixing and different WD masses, and the yellow dotted line representing models constructed by Starrfield et al. (2024) with 25% and 50% mixing and different WD masses. It is worth mentioning that our model provide ⁷Be yields for a range of accretion rates for each $M_{\rm WD}$, rather than sampling only at a fixed accretion rate of 2 $\times$ $10^{-10}$ $\rm M_{\odot}$ yr^-1 as in the other models mentioned above. It is evident that the ⁷Be abundance produced by our model is significantly higher compared to other models and covers almost all CO and ONeMg nova observational samples. In addition, we have made predictions for the two upcoming nova T CrB and M31N 1926-07c eruptions in 2024 and marked the expected abundance of ejected ⁷Be in Figure 2. This result indicates that element diffusion enhances the transport efficiency between the interior and surface convective region of the WD, transporting sufficient ⁷Be to the surface region. Our results are consistent with observations and successfully explain the high ⁷Be abundance observed in the spectra after nova eruptions.

Figure 3: ⁷Li yields vs. WD mass for all CO WDs in the grid.

Due to the intense shock produced by nova eruption, the ash of TNR is partially ejected into the ISM (Tajitsu et al., 2015; Izzo et al., 2015; Molaro et al., 2022, 2023; Starrfield et al., 2024). The ejected ⁷Be during nova eruption decays into⁷Li. This contributes to the total Li elements in the Galactic ISM. Therefore, we consider the abundance of ejected ⁷Be to be ⁷Li. In our simulations, Li mass produced by a nova outbursts is calculated by $M_{\rm Li}=\int\dot{M}X({\rm{}^{7}Li}){\rm d}t$, where $\dot{M}$ is the mass-loss rate of WD during nova outburst (See Eq.1), $X(^{7}{\rm Li})$ is the mass fraction of ⁷Li in the nova ejecta. Figure 2 shows $X(^{7}{\rm Li})$. Compared to the previous theoretical results in José & Hernanz (1998), José et al. (2020) and Starrfield et al. (2024), the Li abundances in our simulations are higher and are closer to the observed values. Figure 3 gives the Li yields for a nova outbursts. Our results are consistent with those of Starrfield et al. (2024). Based on Figures 2 and 3, our model may underestimate the mass-loss rate $\dot{M}$, which is described by Eq.1. Observations indicate that the range of ejection masses falls between 10^-6 to 10^-4  M_⊙, and José & Hernanz (1998) proposed an average ejection mass of $2\times 10^{-5}$  M_⊙ for novae. In our models, we calculate the range of ejection masses falls between $9.17\times 10^{-7}$ to $5.93\times 10^{-5}$ M_⊙. This result is slightly smaller than the observed value. As mentioned by Denissenkov et al. (2021), the increase of ³He in the accretion layer of our WD model can lead to an earlier TNR, a decrease in peak temperature and accreted mass, and ultimately result in a smaller ejection mass. By utilizing the Li and Be abundances and the ejected mass of WDs output by MESA, we can obtain the Li mass for each burst ejection.

The known novae are mainly concentrated in the range of WD mass $M_{\rm WD}$=0.8 M_⊙ - 1.0M _⊙ and accretion rates around $1\times 10^{-10}$ $\rm M_{\odot}$ yr^-1. In this range, our estimating Li yield reaches to about $2\times 10^{-10}$ M_⊙. This yield is slightly lower than the observations of the latest recurrent novae RS Oph (Izzo et al., 2015; Molaro et al., 2023), but it has already reached the same order of magnitude. It can be demonstrated that the nova model constructed in this paper can explain a large part of the observations in the Galaxy.

3.2 Contribution of ⁷Li produced by novae to the Galactic

Theoretically, Rukeya et al. (2017) and Starrfield et al. (2024) estimated the specific contribution of ⁷Li produced by novae to the Galactic, and they found that this contribution is approximately 10%. However novae have been proposed as the main factories of lithium in the Galaxy is something that was proposed in the 1970s by Arnould & Norgaard (1975) and Starrfield et al. (1978), and then demonstrated in Izzo et al. (2015), Tajitsu et al. (2015) and Molaro et al. (2016, 2023), based on the detection of lithium and beryllium in a few novae. As shown by Figure 3, compared with the observations, the theoretical predictions underestimate Li yields produced by nova eruption.

In order to estimated the contribution of ⁷Li produced in our nova models to the Galactic ISM, we use the method of binary population synthesis (BPS), which has been applied by our group (Lü et al., 2006, 2009, 2013; Rukeya et al., 2017; Yu et al., 2019; Lü et al., 2020; Yu et al., 2021; Gao et al., 2022). BPS is a robust approach to evolve a large number of stars (including binaries) so that we can explain and predict the properties of a population of a type of stars (Han et al., 2020). In the population synthesis method for binary systems, input parameters include the star formation rate (SFR), the initial mass function (IMF) of the primaries, the initial mass ratio distribution, the initial orbital distribution, the eccentricity distribution, and the metallicity Z of binary systems. Here, with help of the rapid binary evolution (BSE) code, originating from Hurley et al. (2002), we can rapidly evolve a large-sample binaries into nova binaries, in which $M_{\rm WD}$ and the $\dot{M}$ are given.

Figure 4: Li mass contributions from different sources as a function of Galactic age. The orange line represents our work, while the green and gray dashed lines represent the model results for Rukeya et al. (2017) and Kemp et al. (2022a). The shaded area is the error range of the latter. The results of AGB are represented by purple dotted lines.

Using MESA, we construct a grid for each nova eruption with $M_{\rm WD}$, $\dot{M}$ and ⁷Li yields. For every nova simulated using the BSE code, by a bilinear interpolation of the above two physical quantities ($M_{\rm WD}$ and $\dot{M}$) in MESA, one can then calculate the ⁷Li yields of every nova.

Following Lü et al. (2006), in the method of population synthesis, we use the initial mass function of Miller & Scalo (1979) for the mass of the primary components and a flat distribution of mass ratios (Kraicheva et al., 1989; Goldberg & Mazeh, 1994). A logarithmic flat distribution of initial separations between 10 and 10⁶ R_⊙ is used.

In the binary evolution, the common envelope (CE) evolution is critical important for the formation of nova systems. In BSE code, a combined parameter $\alpha_{\rm CE}\times\beta_{\rm CE}$ is used to calculate common envelope evolution, where $\alpha_{\rm CE}$ is the fraction of binary binding energy which is spent to expel the CE, and $\beta_{\rm CE}$ is a parameter for the envelope structure of the donor. As discussed in Rukeya et al. (2017), $\alpha_{\rm CE}\times\beta_{\rm CE}$ has a weak effect on ⁷Li yields produced by novae. In this work, we take $\alpha_{\rm CE}\times\beta_{\rm CE}=1.0$. Simultaneously, circle orbits and solar metallicity for all binary systems are taken. To calculate the production rate of ⁷Li, we assume a constant star formation rate of 5 M_⊙ yr^-1 over the past 13 Gyr. In the case of a constant star formation rate, one simulates that $10^{6}$ binary systems, in which the primaries are more massive than 0.8 $\rm M_{\odot}$, which causes a statistical error for Monte Carlo simulations lower than 1% in classical novae.

Figure 4 shows the evolutionary trajectory of ⁷Li produced by novae in the Galaxy at a constant star formation rate. According to the results of population synthesis, we obtained an average ⁷Li yield of about 6.4 $\times$ $10^{-11}$ $\rm M_{\odot}$ for each nova eruption, with approximately 8,000 eruptions per nova binary, which is consistent with the previous results (Shara et al., 1986; Molaro et al., 2020). We then select nova systems from binary systems, the primary is a white dwarf (including CO white dwarf and ONeMg white dwarf) and is accreting material from its secondary (a main sequence or red giant star) via Roche lobe or stellar wind. The result is that among the $10^{6}$ binary systems, 17368 of them evolve into nova systems, accounting for approximately 1.7%. This is consistent with Rukeya et al. (2017) ($\sim$1.5%) and Hurley et al. (2002) ($\sim$1.9%). Based on this, we can conclude that the nova rate is 130 per year, which is within the predicted range (Sharov, 1972; della Valle & Livio, 1994; Shafter, 1997, 2002, 2017). Among the $10^{6}$ binary systems, 1.7% of them evolve into nova systems, resulting in an annual ⁷Li production of approximately 8.4 $\times$ $10^{-9}$ $\rm M_{\odot}$. Considering the age of the Galaxy to be $1.3\times 10^{10}$ years, the contribution of novae to ⁷Li is estimated to be around 110 $\rm M_{\odot}$, compared to the total Li content of approximately 150 $\rm M_{\odot}$ in the interstellar medium of the Galaxy (Hernanz et al., 1996; Molaro et al., 2016). Therefore, novae contribute approximately 73% of the ⁷Li in the Galactic ISM.

In Figure 4, the solid red line represents the ⁷Li production at the age scale of the Galaxy, which is significantly higher than the results obtained by Rukeya et al. (2017). Kemp et al. (2022a) used the binary population synthesis code binary_c, to simulate the novae in the Galactic and, by using observed Li production, effectively explained the Li abundance in the early Sun. The ⁷Li production calculated in our model falls within the predicted error range for nova ⁷Li production, as demonstrated by Kemp et al. (2022a), proving that novae are indeed the most significant contributors to Li in the Galaxy. In addition, we also calculate the contribution of asymptotic giant branch (AGB) stars using the same method. The hot bottom burning mechanism of AGB stars can bring Li elements to the surface, which will evolve into the planetary nebulae stage and ultimately be ejected into the ISM (Smith, Verne V. and Lambert, David L., 1990; Romano et al., 2001). After calculation, each AGB star produces approximately $10^{-9}$ $\rm M_{\odot}$ of ⁷Li, resulting in a final contribution to ⁷Li of around 1 $\rm M_{\odot}$, which is only 1%. This result has been consistent with previous research (Romano et al., 2001; Doherty et al., 2014; Rukeya et al., 2017).

In our simulation, the lithium mass produced by novae is $\sim$ 110 M_⊙. It indicates that the novae can provide about 70% of the contribution to the Galactic ISM. However, compared to the total Li mass in the entire Galaxy ($\sim$ 1000 M_⊙) (Molaro et al., 2023), Li produced by the novae is limited (about 15%-20%), and most of Li may originate from the Big Bang.

4 summary

We used the MESA code to construct nova models with element diffusion to calculate the abundance of ⁷Be on the surface of the WD during the occurrence of a TNR, as well as the production of ⁷Li after the eruption. Element diffusion enhances the efficiency of transport between the nuclear reaction zone and the convective region on the surface of the WD, allowing more material to be transported to the surface. Based on theoretical and observational evidence, we appropriately increased the amount of ³He produced during the mixing of material from the donor star and the accretor star, leading to more ⁷Be being produced. A transport channel within the WD is formed due to the increased transmission efficiency caused by element diffusion. This channel effectively transports the ⁷Be from the hydrogen-burning region of the WD to the convective envelope, where it decays into ⁷Li. Therefore, a large amount of ⁷Be produced by the nuclear reaction process in the WD can be transferred to the surface and ultimately ejected. The abundance of ⁷Be on the surface of the WD during the occurrence of a TNR in our model is consistent with the values obtained from the observed nova samples.

When the ⁷Be ejected by the nova eruption decays into ⁷Li in a short period of time, we use population synthesis methods to calculate the production of ⁷Li from the nova eruption. We find that about 1.7% binary systems in the Galaxy can evolve into nova binaries. Every nova binary evenly can occur about 8,000 eruptions, and the occurrence of nova eruption is about 130 yr^-1. Each eruption can produce about 6.4 $\times$ $10^{-11}$ M_⊙ of ⁷Li. There are about 110 M_⊙ Li produced by novae, which is about 73% of total Li mass in the interstellar medium of the Galaxy, and approximately 15%-20% of the entire Galaxy. It means that novae are the primary source of Li in the Galactic ISM, and also one of the important sources of lithium in the entire Galaxy.