© The Authors 2024

1 Introduction

Lithium is certainly the most intriguing of the chemical elements. Being the third lightest element in the periodic table, it holds a unique position in astrophysics due to its pivotal role in probing the chemical evolution of the Milky Way and the properties of the outer stellar layers. Despite the relative simplicity of its atomic structure, the abundance and distribution of lithium in stellar environments remain enigmatic and continue to challenge the understanding of stellar structure, nucleosynthesis, and galactic chemical evolution. The most abundant isotope of lithium, ⁷Li, is the only nuclide known to be produced in three different astrophysical sites: in the hot early Universe by Big Bang nucleosynthesis (BBN), in Galactic cosmic rays (GCR) through spallation and fusion reactions, and in stars by thermonuclear reactions. The contributions to the total Li abundance of the first two sources can be easily quantified and is found to be less than half of the proto-solar Li value (e.g. Prantzos 2012). At the same time, the type of stars that make the main ‘stellar’ contribution to the Li content of the Galaxy remains a matter of debate.

The production mechanism of the light elements Li, Be, and B was an enigma to the founding fathers of stellar nucleosynthesis, who coined the name ‘x-process’ for it (Burbidge et al. 1957). The pioneering studies of Reeves et al. (1970) and Meneguzzi et al. (1971) showed that the x-process is due to GCR interacting with the interstellar medium (ISM) and can produce the totality of proto-solar (meteoritic) ⁶Li, ⁹Be, ¹⁰B, and most of ¹¹B but only ~20% of ⁷Li.

Cameron & Fowler (1971) suggested that ‘... potentially very large amounts of ⁷Li may be produced by the ⁷Be transport mechanism of Cameron (1955)’ in red giant stars. Subsequent theoretical (Scalo et al. 1975; Sackmann & Boothroyd 1992; Forestini & Charbonnel 1997) and observational (Smith & Lambert 1989, 1990) studies showed that this could occur in asymptotic giant branch (AGB) stars. Specifically, in the ‘hot bottom’ – temperatures T~(40–60) × 10⁶ K – of their convective envelopes, where the fusion of ³He with ⁴He produces ⁷Be, which is convectively transported to the surface while beta-decaying within a timescale of 53 days to ⁷Li.

In the meantime, Spite & Spite (1982) and Spite et al. (1984) found that the Li abundances in dwarf halo stars of Population II display a remarkably constant value for various temperatures (unlike in Population I stars), and they should thus reflect the original Li content of those stars. Moreover, the lithium abundance in the hotter halo dwarfs appeared to be very uniform (independent of metal abundance, spatial velocity, eccentricity, and galactic stellar orbit), arguing in favour of the Big Bang origin of lithium, as calculated by Wagoner et al. (1967). Up to a metallicity [Fe/H]~−1.2 (see Fig. 3 in Rebolo et al. 1988), the ‘Spite plateau’ abundance of Li is approximately one-tenth of the meteoritic one, and combined with the known contribution of the GCR, it suggests that the dominant source of Li is stellar.

Beyond AGBs, other candidate stellar sources of Li include

• novae (Arnould & Norgaard 1975; Starrfield et al. 1978; José & Hernanz 1998), where explosive H-burning occurs episodically at the surface of a white dwarf accreting material from a companion star;

• core-collapse supernovae (CCSNe), where energetic μ- and τ- neutrinos from the core excite the He nuclei of the outer layers and induce reactions between them and ³He or ³H, leading to the production of ⁷Li (Woosley et al. 1990);

• low-mass red giant stars, via extra deep mixing and the associated ‘cool bottom processing’, with the amount of Li produced depending critically on the details of the extra mixing mechanism (Sackmann & Boothroyd 1999).

The Galactic chemical evolution of Li has mostly been studied with one-zone models, adopting one or all of the aforementioned stellar sources and considering the Spite plateau value of Li as the primordial one (e.g. Abia & Canal 1988; Abia et al. 1993, 1995; D’Antona & Matteucci 1991; Matteucci et al. 1995; Travaglio et al. 2001; Romano et al. 2001; Prantzos 2012). Though no robust conclusions have been reached as to the main stellar source, it was shown that the yields of all candidate sources are insufficient to explain the proto-solar Li by a large factor (see Prantzos 2012, and references therein). Moreover, the observational constraints are weak, since Li is prone to depletion in stellar envelopes, and its photospheric abundance does not trace the initial one, except in A-type stars (see Charbonnel et al. 2021; Randich & Magrini 2021, and references therein). The ‘upper envelope’ of the Li observations can be used, in principle, as a tracer of the true Li evolution (Rebolo et al. 1988), or at least as a lower limit describing that evolution, but it is difficult to define it in practice (Lambert & Reddy 2004).

In the past years, novae attracted considerable attention as a source of Li (Cescutti & Molaro 2019; Grisoni et al. 2019; Matteucci et al. 2021; Romano et al. 2021) due to the detection of substantial abundances of ⁷Be and, perhaps, ⁷Li in their ejecta (Tajitsu et al. 2015; Izzo et al. 2015; Molaro et al. 2016; Izzo et al. 2018; Molaro et al. 2023, and references therein). On the other hand, the recent work of Kemp (2022) provided for the first time theoretical delay time distributions (DTDs) for novae through population synthesis models of binary star evolution. Kemp et al. (2022a) used those DTDs and observationally inferred Li yields to study the Li evolution with a one-zone model and to show that novae can be the only stellar source of Li.

In light of these recent developments, we reassess the problem of the Galactic evolution of Li, and we compare our results to observations of field stars from the GALAH survey and of open clusters from the Gaia-ESO survey. In Sect. 2 we present the observational data. In Sect. 3 we first discuss in some depth the nova DTDs (which are metallicity dependent) and the corresponding Li yields. Then we use a one-zone model to show that AGBs are negligible contributors to Li with current Li yields, while novae can be the sole stellar source. We also evaluate the contribution of BBN, GCR, and novae to the Galactic Li content. In Sect. 4 we use a multi-zone model with radial migration to discuss the impact of the adopted nova DTDs on the evolution of Li in the various Galactic zones, paying particular attention to the metal-rich stars formed in the inner disc. We summarise our results in Sect. 5.

2 Observational data

2.1 Upper lithium envelope

In exploring the chemical evolution of the Milky Way, it is crucial to take into account the evolution of Li in the stars that are used to constrain the models. The abundance of this fragile element is subject to significant depletion within stars, and stars rarely exhibit in their photospheres the Li content they inherited from their protostellar cloud (see Borisov et al. 2024, and references therein). To track the original abundance of Li in stars in the Galaxy, we need to turn our attention to the stars that might preserve it as close as possible to the initial content, and the best candidates are warm stars on the hot side of the so-called Li-dip with the effective temperature above ~6800 K (Randich et al. 2020; Charbonnel et al. 2021). These hot and relatively massive stars (masses above ~1.3 M_⊙ at solar metallicity) have shallow convective envelopes, which prevents the transport of Li to the inner hot layers where it can be destroyed by the proton capture reaction ⁷Li(p,α)α. Unfortunately, these stars are less numerous than G-type stars due to the initial mass function (IMF), and they are also rarely included in surveys, as they are fast rotators for which it is more difficult to determine abundances.

Unfortunately, the selection of hot dwarfs is complicated by the fact that at lower metallicities such stars are scarce or simply absent, because of the strong dependence of their lifetimes on mass and metallicity, with more metal-poor stars being older, on average, than metal-rich stars. At low metallicity ([Fe/H] below ~−1.5 dex), we are thus left with the main sequence or turnoff stars that have undergone significant Li depletion along their evolution. For this reason, we expect an increase in the difference in A(Li) between the predictions of Galactic evolution models and the observed stellar Li values, especially at low metallicities.

2.2 Selection of the sample

We selected dwarf field stars from GALAH DR3 (Buder et al. 2021) in a wide metallicity range from −2.84 ≤[Fe/H] ≤ 0.5 (lower limit on metallicity corresponds to the most Fe-poor star in GALAH DR3). We applied the quality flag on stellar parameters flag_sp=0 (for T_eff and log g). When set to zero, it also eliminates possible binaries. In addition, we applied a flag on iron and lithium abundance: flag_fe_h=0 and flag_li_fe=0. As recommended by Buder et al. (2021), we selected stars with S/N ≥ 30. We kept only those objects that have a relative uncertainty of age σ_age/age≤ 10% and low interstellar extinction A_G ≤ 0.2. The latter allowed us to avoid the influence of interstellar extinction on the T_eff and L determination, which leads to incorrect age determination.

Finally, to focus on dwarf stars and avoid possible contamination from (sub)giants, we made the preliminary selection of stars with log g≥3.5 and then followed a similar selection algorithm as in Borisov et al. (2024, see their Sect. 4.1.3). Briefly recapitulate the basic mechanism of the algorithm: we use a grid of stellar evolutionary tracks with the determined positions of the ZAMS and TAMS (terminal-age main sequence) in the Hertzsprung-Russell Diagram (HRD) for each track. For each chemical composition of the grid, we define a polygon with ZAMS and TAMS points as its vertices. We assume that a star is a dwarf if its position in the HRD (T_eff and Gaia DR3-based luminosity L) falls into this polygon. To take into account the discreteness of the grid, we perform bicubic interpolation of vertices on a star’s [Fe/H] and [α/Fe]. Here we used the BaSTI tracks (Pietrinferni et al. 2004, 2006, 2024) as they have wide coverage in T_eff space. As a result, the final sample of field stars contains 2136 objects with 3D non-local thermodynamic equilibrium Li abundances. The hottest stars of the sample are on the hot side of the Li-dip and correspond to the ‘warm group’ referred to in Gao et al. (2020). We used the ages of stars provided in the GALAH value-added catalogue and determined with the code BSTEP by Sharma et al. (2018) using PARSEC+COLIBRI isochrones (Marigo et al. 2017). The wide age coverage of the GALAH sample, which ranges from 25 Myr to 12.8 Gyr, allowed us to challenge the predictions of our Galactic evolution models over almost the entire history of the Milky Way.

The vast majority of the sample field stars are from the solar vicinity: 95% are within 1 kpc from the Sun, while the median distance is ~590 pc. This did not allow us to compare our multi-zone model predictions (see Sect. 4) with observations at different galactocentric distances R_GC. For this reason, we also selected open clusters from Romano et al. (2021) that cover a wide range of R_GC from 5.2 to 15 kpc. In the paper, the authors provide mean values of A(Li) – not corrected for non-local thermodynamic equilibrium effects – for the stars that have suffered minimal Li depletion. In most cases (especially in the case of clusters older than 0.1 Gyr), this corresponds to stars from the hot side of the Li-dip that are, on average, slightly hotter than the Li-richest stars from the GALAH sample. Otherwise, the measurements are based on the Li content of pre-main sequence stars.

3 One-zone chemical evolution model for the solar neighbourhood

3.1 Description of the model

Our simple, one-zone model, is introduced in detail in Prantzos et al. (2018) and here we present its main features. The local disc is built by infall of gas at an exponentially decreasing rate and a characteristic time-scale of 8 Gyr, where the star-formation rate (SFR) Ψ is given by a Schmidt-Kennicutt law: $$\mathrm{\Psi }(t)=C{\mathrm{\Sigma }}_G(t),$$(1)

where ∑_G is the local gas surface density and the coefficient C is chosen to obtain a gas fraction of ~20% at the end of the simulation, compatible with the one presently evaluated in the solar neighbourhood. We assume that the duration of the system is 12 Gyr.

We use the metallicity-dependent stellar lifetimes and yields of Cristallo et al. (2015) for Low and Intermediate Mass stars (1–7 M_⊙) and of Limongi & Chieffi (2018) for the massive ones (M >13 M_⊙), which include mass loss and rotation. No Li is produced in those models (see Sect. 3.3). Stars above 25 M_⊙ are assumed to collapse into black holes and enrich the interstellar medium only with their wind ejecta. We interpolate the yields in the region 7–13 M_⊙. For the massive stars, we use an Initial Distribution of Rotational Velocities (IDROV), calibrated on observations of abundances in halo and disc stars (see Prantzos et al. 2018, for details). We adopt the stellar IMF of Kroupa (2001) in the mass range 0.1–120 M_⊙ and we calculate chemical evolution with the Single Particle Population (SSP) method.

We adopt metallicity-dependent yields of Type Ia supernovae (SNe Ia) from Leung & Nomoto (2018) and we consider the DTD of SN Ia as discussed in Maoz & Graur (2017), that is, a power law in time with index x = −1.1. We imposed a lower cut-off at 40 Myr, which corresponds to the shortest possible time between the formation of a binary system and the formation of a white dwarf in it. At a given moment t, we computed the SN Ia rate as a convolution of the SFR Ψ(t) and the SN Ia DTD: $$R_{\mathrm{I}\mathrm{a}}(t)={\int }_0^t\mathrm{DTD} \left(t^{\mathrm{\prime }}\right)\mathrm{\Psi }(t-t^{\mathrm{\prime }})d{t}^{\mathrm{\prime }}.$$(2)

It should be noted that the one-zone model presented here was developed only for the disc of the Milky Way. Despite this fact, it can serve as a rough approximation for capturing the general evolutionary trends in the Galactic halo since the halo phase is characterised by metallicities [Fe/H]≤−1, a duration of approximately 1 Gyr, and a ~constant [α/Fe] ratio during that period. These features are also met in the one-zone model of Prantzos et al. (2018) adopted here.

3.2 Production of Li: GCR and nova yields and DTDs

We model the evolution of Li starting with an initial abundance of A(Li)=2.75 dex that corresponds to the SBBN value of Pitrou et al. (2018). This Li is injected into the system during the evolution through the infall of primordial composition.

We compute Li production by GCR as discussed in detail in Prantzos (2012), in a fully self-consistent manner: (a) we consider the energy released by supernovae (10% of the kinetic energy of 1.5 × 10⁵¹ erg from both CCSNe and SNe Ia); (b) we assume that this energy is distributed across an injection spectrum of GCR (power-law in momentum space) accelerated near each SNe; (c) we consider the GCR spectra at equilibrium after injection spectrum is propagated with an escape length Λ=10 g/cm², (d) we take into account the composition of the GCR as resulting from the composition of the stellar winds in the case of CCSNe (properly weighted with the adopted IMF) and from the composition of the ISM in the case of SNe Ia; (e) the GCR particles of H–He–C–N–O with the aforementioned equilibrium spectra and composition interact with the ISM H–He–C–N–O with energy-dependent cross-sections and produce the isotopes ⁶Li, ⁷Li, ⁹Be, ¹⁰B and ¹¹B through spallation of CNO nuclei with protons and alphas as well as the fusion of two α-particles (in the case of Li isotopes). The resulting ‘yields’ $y_6^{GCR}$ and $y_7^{GCR}$ (in M_⊙) are time-dependent because they involve the composition of the ISM and of the GCR. The mono-isotopic ⁹Be is used as a test of the GCR composition, since observations of [Be/Fe] versus [Fe/H] in the past 30 years suggest that this composition has little changed over time (Prantzos 2012, and references therein)¹.

Nova eruptions are produced when H-rich material from a companion star is piled up on the surface of a white dwarf (composed of either C-O or O-Ne) and reaches ignition in degenerate conditions, at sufficiently high pressures. A fraction of the H-rich material is convectively mixed with the C-O-rich material of the white dwarf and H-burning takes place at temperatures of 200–300 MK through the hot CNO cycle. Nuclear reactions between the accreted ³He and ⁴He produce unstable ⁷Be that decays into ⁷Li with a half-life of ~53.22 days. The fragility of ⁷Be makes it almost impossible to survive at such high temperatures. But the rapid evolution and high mass loss rate during the nova eruption, allows ⁷Be to be transported away from the hot region and be ejected in the ISM.

Although novae were proposed a long time ago as a key stellar source for Li production, the quantitative evaluation of its importance was hampered by two factors: the very uncertain Li yields of novae and the unknown DTD of those sources. Various empirical DTDs for novae – with a nova rate proportional to the formation rate of whited dwarfs plus some time-delay – were used in several papers (D’Antona & Matteucci 1991; Romano et al. 1999, 2001; Grisoni et al. 2019), while Rukeya et al. (2017) adopted a population synthesis code to model the nova rate. The very recent work of Kemp (2022) – also presented in Kemp et al. (2022a,b) – provides a very promising theoretical framework for the nova DTD, and we adopt it here.

In Fig. 1, we display the metallicity-dependent DTDs of Kemp et al. (2022b) for both novae events and ejected masses because they have not the same DTDs. In that respect, the case of novae differs from the one of SN Ia, in which only the DTD of SN Ia rate is required, while the ejected mass is taken always to be approximately the Chandrasekhar mass. We show DTDs of nova events and ejecta for 1 M_⊙ of stars formed and for a few selected metallicities Z in the two top panels, colour-coded with their metallicities. The curves are fitted with power-laws and the corresponding slopes appear also in the panels. Those slopes, and in particular the ones for the ejecta, are systematically smaller – in absolute value – than the slope of x=-1.1 at late times for SN Ia (with grey curves in both figures). This implies that nova production of Li will rise at late times more rapidly than Fe production from SN Ia, and thus the [Li/Fe] ratio is expected to rise with metallicity in those times, as is discussed in Sect. 4.1.

In the third panel of Fig. 1, we display the average mass of the nova ejecta, obtained by dividing the curves of the second panel by those of the first one. This average mass is time-dependent (increasing with time), because of the systematic difference in the slopes of the previous panels. The time-integrated average mass of the ejecta (over 12 Gyr) is ⟨M_EJECTA⟩) ~ 5 × 10⁻⁵ M_⊙ per nova event and varies little between the various metallicities. Compared to the results of recent hydrodynamical nova explosions from Starrfield et al. (2024), these masses are several times higher, especially for the higher metallicities. In contrast, the 3D hydrodynamical models of José et al. (2020), which appear within grey shaded areas in the figure, obtain values in the range of those provided in Kemp et al. (2022a). Those values are within the range of observations (indicated on the right of the third panel), which display a rather wide diversity.

The work of Kemp et al. (2022a) adopts a novel approach in the study of the nova impact on the Galactic chemical evolution (GCE), by introducing the correspondence of ejecta DTDs to different regions of the parameter space of a white dwarf mass M_WD versus ${\dot{M}}_{accretion}$. This allows them to map physics-dependent yields from hydro-simulations to the relevant ejecta rates of their own calculations. It turns out, however, that the Li yields of current nova models are notoriously insufficient to produce the meteoritic value of Li. For that reason, Kemp et al. (2022a) adopt Li mass fractions from observations (see below) and they evaluate Li yields by multiplying those mass fractions with their own ejecta masses; we adopt that approach here. The obtained yields are provided in the bottom panel of Fig. 1 and are considerably higher than those of the hydro models of Starrfield et al. (2024) or other recent studies (Gao et al. 2024; José et al. 2020; Rukeya et al. 2017). On the other hand, despite the high ejecta masses found in the new 3D hydro-simulations of José et al. (2020), the authors find a serious decrease (by more than a factor of 50) in the Li amount with respect to their previous calculations.

At this point, it should be noted that current hydrodynamical models of nova explosions (1D, 2D or 3D) suffer from various shortcomings due to the poorly understood mixing processes (shear, diffusion, Kelvin-Helmholtz instabilities) that operate in the interface between the core (CO or ONe WD) and the accreted envelope (H-rich). The mixing is responsible for the metallicity enhancement observed in the nova ejecta, and its modalities determine critically the amount of surviving ⁷Be (see José et al. 2020).

The discrepancy of current hydro-models and observations regarding the Li mass fraction in novae appears very clearly in Fig. 2. Only the most massive ONe white dwarf models of Starrfield et al. (2024) produce a Li abundance compatible with the lowest range of the observed values. A surprising exception comes from the very recent work of Gao et al. (2024) who find ⁷Be mass fractions from CO nova compatible with the observed ones. The authors claim that this is due to the high mass fraction of ³He assumed in their calculations (an idea suggested in Molaro et al. 2020) combined to their treatment of element diffusion, while acknowledging that Denissenkov et al. (2021) find the opposite effect, namely that an enhanced amount of ³He reduces the ⁷Be mass fraction. However, the ejecta masses in the models of Gao et al. (2024) are quite small and so are their corresponding Li yields (see bottom panel in Fig. 1).

In those conditions, it appears to us that the theoretical Li yields are by far insufficient to explain the proto-solar Li value at present. For that reason, we adopt here the mean ⁷Be abundance observed in novae, estimated by Molaro et al. (2022, 2023) (dotted horizontal line in Fig. 2) as A(⁷Be)=7.34 ± 0.47, which corresponds to the Be (and, thereof, Li) mass ratio over H of $X_{\mathrm{L}\mathrm{i}}/X_{\mathrm{H}}=7\times {10}^{\mathrm{A}(\mathrm{B}\mathrm{e})-12}=\left({1.53}_{-1.01}^{+2.99}\right)\times {10}^{-4}$. To convert this ratio into the required Li mass fraction, the corresponding mass fraction of H is needed and this depends on the nova model. Nova modellers assume mixtures of white dwarf material (H-free) and accreted material (with solar H), of either 25–75% or 50–50 % (e.g José & Hernanz 1998; Starrfield et al. 2024), leading to X_H ≈ 0.5 for the former case and X_H ≈ 0.35 for the latter. The corresponding Li mass fractions in the ejecta are then $X_{\mathrm{L}\mathrm{i}}=\left({7.7}_{-5.1}^{+14.9}\right)\times {10}^{-5}$ and $X_{\mathrm{L}\mathrm{i}}=\left({5.4}_{-3.5}^{+10.5}\right)\times {10}^{-5}$. In order to calculate the Li production rate here, we used the former value multiplied by the nova ejecta DTDs (second panel in Fig. 1) after folding them with the SFR. On the other hand, the time-integrated average values of the ejecta reported previously (5 × 10⁻⁵ M_⊙) multiplied by the above values lead to an average value of the Li yield of a nova of (3–4) × 10⁻⁹ M_⊙, 15–20 times higher than the maximum Li yields obtained in all recent studies (Starrfield et al. 2024; Gao et al. 2024; José et al. 2020; Rukeya et al. 2017), as indicated in the bottom panel of Fig. 1.

Fig. 1 Derived delay time distributionss for novae event rates $\dot{N}$ (top) and ejecta mass rates $\dot{M}$ (second from top) per 1 M⊙ of stars formed and for a few selected metallicities, adopted from Kemp et al. (2022b). The slopes X of power-law fits appear colour-coded for each metallicity. The grey curves display the corresponding DTDs for SN Ia, with a single slope x=-1.1 above t=0.3 Gyr; the SN Ia curves are multiplied by 10 000 in the top panel and divided by 100 in the second panel. Third panel: average ejecta mass of Kemp et al. (2022a) compared to the highest ejecta mass of the models of Starrfield et al. (2024) indicated by upper limits, to hydrodynamical models of José et al. (2020) within grey shaded areas, and to available observations from S07 (Schwarz et al. 2007), M22 (Molaro et al. 2022), P22 (Pandey et al. 2022), and DM15 (Das & Mondal 2015). Bottom panel: average Li yield when assuming the ejecta mass of the third panel and an average Li mass fraction ⟨XLi⟩OBS from observations (see Fig. 2 and text). The grey arrows indicate the highest Li yield obtained in any of the cited studies.

Fig. 2 Relative abundances of Li/H in nova explosions. The shaded area shows the range of observed values and the red horizontal line is the average value, as reported in Molaro et al. (2022). Curves represent the highest values found in the calculations of Starrfield et al. (2024, in red), Rukeya et al. (2017, in green) and Gao et al. (2024, in blue), in explosions of CO and ONe white dwarfs of various masses. The results of José et al. (2020) are considerably below the scale of the figure.

3.3 Results based on the one-zone model

In this paper, we focus on novae as the key stellar source of Li production. However, within the context of the one-zone model, we initially explore the possibility of AGB stars as a potential source of Li. This possibility has been studied in the past (e.g. Travaglio et al. 2001; Romano et al. 2001; Prantzos 2012) and rejected, since the then available AGB yields were too small to explain the pre-solar (meteoritic) Li value. Cristallo et al. (2015) do not provide Li yields and we assess here briefly that point in the light of the metallicity- and mass-dependent AGB yields of Li from Karakas & Lugaro (2016).

In our approach, we assume that AGB producers of Li cover the range of initial masses from 1.5 M_⊙ to M_max. We compute then models with different values of M_max and flat yields X_Li in the corresponding mass range. The value of A(Li) at solar metallicity and birth-time (4.56 Gyr ago) is shown colour-coded in Fig. 3 as a function of both parameters, M_max and X_Li. The dashed line in the figure indicates the combinations of M_max and X_Li that provide the meteoritic value of A(Li)=3.31 at [Fe/H]=0. The required Li yields are slightly higher than 10⁻⁷ M_⊙ in the range from 1 M_⊙ to M_max. Such high yields are found only in the supersolar-metallicity AGB stars of Ventura et al. (2020) but only for masses higher than 7 M_⊙, while they are much smaller at lower masses. In fact, below 5 M_⊙ stars display no hot bottom burning and rather destroy their initial Li in those models (depending on the metallicity).

In the same figure, we show the average yields of AGBs with initial masses from 1.5 M_⊙ to M_max and different metallicities according to Karakas & Lugaro (2016). For each value of M_max, the average yield is computed as an IMF-weighted mean: $$⟨X_{\mathrm{L}\mathrm{i}}⟩=\frac{{\int }_{1.5{\mathrm{M}}_{\odot }}^{{\mathrm{M}}_{max}}{X}_{\mathrm{L}\mathrm{i}}(\mathrm{M}){\mathrm{M}}^{-\alpha }\mathrm{d}\mathrm{M}}{{\int }_{1.5{\mathrm{M}}_{\odot }}^{{\mathrm{M}}_{max}}{\mathrm{M}}^{-\alpha }\mathrm{d}\mathrm{M}},$$(3)

where α = 2.3 according to Kroupa (2001). The average yield of AGBs obtained that way is far below the one required to produce the meteoritic Li value at [Fe/H]=0, by almost 2 orders of magnitude on average. Thus, AGB stars cannot be the main stellar source of Li, at least with presently available yields. In fact, the discrepancy in that respect is even larger than in Prantzos (2012), who had adopted older AGB yields.

Using the same one-zone GCE model, we studied the case of nova as a sole stellar source of Li in the Galaxy by using the metallicity dependent DTDs from Kemp et al. (2022b) for nova rates of events and ejected masses (as shown in Fig. 1). From our modelling, we find (Fig. 4) that the nova mass fraction that best fits the solar meteoritic value of Li is X_Li=1.0 × 10⁻⁴, which is consistent with the observational average estimate of Fig. 2 within ~0.25σ of A(⁷Be). However, the large dispersion of A(⁷Be) gives room for variation in the mixing percentage, as indicated with the shaded area in Fig. 4. The four curves were obtained with four different values of Be (Li) mass fraction: the best-fit value that we found, the mean value based on the estimation by Molaro et al. (2022) and 25% mixing, and the maximum possible values of the mass fraction from hydrodynamical simulations by José & Hernanz (1998) and Starrfield et al. (2024). The mean observational value and maximum possible hydrodynamical values are clearly insufficient to make the meteoritic Li value in the framework of our model (see also Kemp et al. 2022a).

In Fig. 5 (top panel), one sees that GCR reproduce correctly the meteoritic abundance of ⁶Li a solar metallicity, a key test regarding the treatment of the GCR component. Following Prantzos (2012), we present the relative contributions of the various Li sources in Fig. 5 as a function of time in the framework of this one-zone model. Those contributions correspond to the time-integrated masses of Li ejected during the evolution by each of its three different sources (BBN, GCR, and novae) and were calculated as follows: $$C_{BBN}(t)={\int }_0^t{\dot{m}}_{INF}\left(t^{\mathrm{\prime }}\right)X_7^{BBN}d{t}^{\mathrm{\prime }},$$(4)

where ${\dot{m}}_{INF}(\mathrm{t})$ (M_⊙/Gyr) is the rate of infall with primordial composition; $$C_{GCR}(t)={\int }_0^t{R}_{SN}\left(t^{\mathrm{\prime }}\right)y_7^{GCR}\left(t^{\mathrm{\prime }}\right)+y_6^{GCR}\left(t^{\mathrm{\prime }}\right)]d{t}^{\mathrm{\prime }},$$(5)

where R_SN(t) is the total supernova rate (CCSN + SN Ia) per Gyr and $y_6^{GCR}$ and $y_7^{GCR}$ (in M_⊙) are the (time and metallicity dependent) GCR ‘yields’, as discussed in Sect. 3.2; and $$C_{NOVA}(t)={\int }_0^t{\dot{m}}_{NOVA}(z\left(t^{\mathrm{\prime }}\right),t^{\mathrm{\prime }}){⟨X_{\mathrm{L}\mathrm{i}}^{NOVA}⟩}_{OBS},d{t}^{\mathrm{\prime }}$$(6)

where ${\dot{m}}_{NOVA}(\mathrm{z}(\mathrm{t}),\mathrm{t})$ is the nova ejecta rate in M_⊙/Gyr, obtained by the convolution of the SFR and the metallicity-dependent nova ejecta DTDs of Kemp et al. (2022b), and $⟨X_7^{NOVA}⟩$ is the average mass fraction of Li observed in novae, as discussed in Sect. 3.2.

During the early period of the Galactic evolution, primordial Li dominates the total Li abundance; this Li is in the form of Li-7 and is contained both in the original gas and in the gas of primordial composition continuously accreted by the Galaxy. The contribution of GCR (approximately one-third of that in the form of Li-6 and the other two-thirds as Li-7, by mass) reaches ~0.1 % of the total around [Fe/H]~−2.5, while the contribution of nova reaches that level around [Fe/H]~−2. The sum of those two exceeds the primordial contribution only around [Fe/H]~−0.5 when the metallicity was slightly above one-third solar. The relative contributions to the proto-solar (meteoritic) Li are approximately: 28% from primordial nucleosynthesis, 17% from GCR, and 55% from nova. Thus, although the stellar component definitely dominates the abundance of Li at Sun’s formation, the other two components almost match it. But 4.5 Gyr later (today), the stellar contribution of novae is increased by 10% and the BBN contribution is decreased by a similar amount. Today, the nova component clearly dominates, having made approximately two-thirds of the Li present in the Galaxy.

We note that in Prantzos (2012), the primordial contribution to Li in the proto-Sun was found to be smaller than here (12% instead of 28%), mainly because in that work the adopted primordial value of Li was A(Li)=2.6 instead of A(Li)=2.75 considered here and the meteoritic value A(Li)=3.4 instead of 3.3 here.

Fig. 3 Average Li yields of AGB stars in the range 1.5 to Mmax as a function of Mmax (according to Eq. (3)) from Karakas & Lugaro (2016). The data are colour-coded according to the initial stellar metallicity. They are compared to the average yield required to obtain the meteoritic Li value at [Fe/H]=0 (dashed curve).

Fig. 4 Li abundance as a function of [Fe/H] with nova as a sole stellar source and with different mass fractions of Li in the nova ejecta. Our ‘best-fit’ value (red dashed line) is slightly higher than the mean value of the observational analysis of Molaro et al. (2022) (red solid line), and the pinkish shaded area corresponds to the dispersion of XLi around that mean value. The theoretical nova Li yields of José & Hernanz (1998) and Starrfield et al. (2024) are also used (purple and cyan curves, respectively). The five-pointed stars show the abundance of the open clusters, while the underlying number density plot represents the GALAH sample. The solar symbol indicates the meteoritic value.

Fig. 5 Evolution of Li abundance (top) and relative contributions (by mass) of different Li sources (bottom) as a function of [Fe/H]. In the top panel, the evolution of total abundance is shown in black, and the solar symbols indicate the solar meteoritic values for 6Li and 7Li. The other curves (in both panels) show the components from the BBN, novae, and GCR (two curves for GCR show components both for 6Li and 7Li). In the bottom panel, the solar symbols indicate the contributions at solar metallicity.

4 Multi-zone chemical evolution models for the disc of the Milky Way

We extend our analysis to the study of the whole disc by using a multi-zone model, which offers more insight into the chemical evolution of the Galactic disc. Such models use more constraints, both global (e.g. the present-day total mass of gas and stars, rates of star formation, infall, supernova, and nova, etc.) and local ones (e.g. radial profiles of gas, stars, and metals). These constraints result from the variation of various physical processes across the disc and enable a more detailed study of its history.

We use the code described in the recent work of Prantzos et al. (2023), which is an update of Kubryk et al. (2015). Here we present briefly its main features. The Galactic disc is formed by the infall of primordial gas into the gravitational well of an evolving dark matter halo, with a final mass of ~10¹² M_⊙. The gas infall timescales are shorter in the inner regions, leading to the inside-out formation of the disc. Stars are formed from molecular gas, the amount of which is calculated from the amount of total gas through semi-empirical prescriptions. Stars undergo radial displacements due to epicyclic motions (‘blurring’) as well as variations in their guiding radius (‘churning’). Radial migration is modelled as a diffusion phenomenon, with coefficients derived from N-body simulations. As initially discussed in Kubryk et al. (2013), radial migration influences not only ‘passive tracers’ of chemical evolution (i.e., long-lived stars preserving the chemical composition of their birth gas in their photospheres), but also ‘active agents’ of chemical evolution, namely long-lived nucleosynthesis sources, such as SN Ia producing iron, low-mass stars producing s-process elements and neutron-star mergers as producers of r-elements. Lithium has now been added to that list since novae are considered the main producers here. The stellar IMF and the nucleosynthesis part of the model are the same as in the one-zone model described in the previous section and we focus here only on novae as stellar sources of Li.

Fig. 6 Evolution of total SFR (orange), as well as rates of novae (black), SNe Ia (red), and CCSNe (blue) in the Milky Way predicted by our model. In the legend, the final model values (at 12 Gyr) are indicated. The SFR values are on the right axis (linear scale). Points with error bars show the present-day estimated values (see Sect. 4.1 for references).

4.1 Results of the multi-zone model

Some global results of the multi-zone model are displayed in Fig. 6. The total SFR (linear scale on the right) is rather high (almost 10 M_⊙/yr) for the first Gyr, being fuelled mostly by the gas of the inner disc, while it declines smoothly at late times to a value ~2.4 M_⊙/yr, compatible with the observed ones (e.g. Elia et al. 2022; Chomiuk & Povich 2011; Licquia & Newman 2015) reported the values of 2.0 ± 0.7, 1.9 ± 0.4, and 1.65 ± 0.19 respectively (all three values are shown in Fig. 6). The same trend is followed by the CCSN rate (logarithmic scale on the left), while the SN Ia rate stays almost constant during the whole evolution, because of the DTD of SN Ia, as weighted by the SFR history of the various zones. In both cases, the final values of those rates are compatible with the present-day values of those quantities, reported in the figure (R_CCSN = 1.79 ± 0.55 yr⁻¹, Rozwadowska et al. 2021 and R_Ia = 0.72 ± 0.23 yr⁻¹, Maoz et al. 2014).

Most importantly for our study, the total nova rate – which is dominated by the inner regions and declines only by a factor of two in the last 8 Gyr – reaches a final value of 36 novae yr⁻¹, in fair agreement with the observationally inferred ones (e.g. Shafter 2017, De et al. 2021, and Kawash et al. 2022 provide estimations of ${\text{R}}_{\text{nova }}={50}_{-23}^{+31},{44}_{-9}^{+20}$, and 26 ± 5 novae yr⁻¹, respectively). It is the first time – to our knowledge – that the nova event DTD of Kemp et al. (2022b) is used in a multi-zone model tailored to reproduce the properties of the Milky Way disc (Prantzos et al. 2023) and reproduces correctly the present-day nova rate of our Galaxy (within observational uncertainties). This is encouraging, both regarding the work presented for the one-zone model (in which the present-day theoretical nova rate could not be checked by local observations) and what follows in this section.

Although the Galactic rate of nova events appears to be ~ constant in time in Fig. 6, this is not true for the Galactic rate of the mass of nova ejecta: the slope of their DTD with time (Fig. 1) is considerably smaller in absolute value than the one of the nova event DTD. As a result, the nova ejecta rate increases with time, especially in the inner Galactic zones, leading to a continuous increase of the Li abundance, in contrast to the outer zones (see Fig. 7).

In Sect. 3.3 we adjusted the Li mass fraction of novae to reproduce the value of Li/H in the model 4.56 Gyr ago (observed meteoritic value). In the multi-zone model, the same should be done for the zone of Galactocentric distance R_GC=8 kpc, in the absence of radial migration. However, several observational indices suggest that the Sun was not formed in that zone but rather in the inner disc and subsequently migrated to its current position.

Indeed, observations show that the gaseous abundances of several elements within ±1 kpc from the Sun are solar to within ±0.04 dex (Cartledge et al. 2006). Also, observations of young B-stars both in the field and in the nearby star-forming region of Orion (Nieva & Przybilla 2012), find solar abundances with small dispersion. Finally, recent observations of Ritchey et al. (2023) find clear evidence that the dispersion in the metallicities of neutral interstellar clouds in the solar neighbourhood is small (±0.10 dex) and only slightly larger than the typical measurement uncertainties. Taking into account the observed abundance gradient, the above results suggest that the Sun migrated to its present-day position from its birthplace, located a couple of kpc inwards in the disc (e.g. Nieva & Przybilla 2012). Theoretical arguments for such a displacement, invoking radial migration, have been proposed (e.g. Wielen et al. 1996; Minchev et al. 2014; Kubryk et al. 2015). As discussed in Prantzos et al. (2023), we find in the framework of this model that [Fe/H]=0 is obtained 4.56 Gyr ago in the gas of the zone at Galactocentric radius R_GC ~6 kpc. We then use the results for that zone for the calibration of the nova Li mass fraction.

The results appear in Fig. 7 where the evolution of the gas abundances for Fe (top) and Li (bottom) is presented for zones at various galactocentric distances. For the zone at 6 kpc, we find a best-fit value X_Li = 7.2 × 10⁻⁵, that is, about two-thirds of the one determined in the one-zone model but still well within the range of observed values reported for novae in Molaro et al. (2022). Taking into account the uncertainties involved in our model (including, among others, the actual extent of the radial migration of novae in the Galaxy) as well as in the observations, we find that the result is quite satisfactory and strengthens the confidence in the role of novae as major stellar contributors to Li enrichment in the Milky Way.

An interesting feature is that in the inner disc, the Fe abundance practically stops increasing at late times, while the abundance of Li increases continuously to very high values. The former is due to the large amounts of gas returned from low-mass stars, which are abundantly produced in the first few gigayears from the intense star formation in those regions; this gas, released at late times, is metal-poor and dilutes the local metallicity. In contrast, Li is less diluted because the gas returned from low-mass stars is Li-free, but the infalling gas contains primordial Li. Moreover, the DTD of nova ejecta declines more slowly with time than the SN Ia DTD (see second panel in Fig. 1); as a result, the strong early activity of those inner regions produces a larger ratio of nova/SN Ia ejecta at late times, enhancing more and more the Li abundance compared to the one of Fe.

Another interesting feature, related to the above, is that for the same metallicity, Li is more abundant in the outer zones than in the inner ones. This, rather counter-intuitive, result was first obtained in Prantzos et al. (2017) and is due to the fact that the metallicity [Fe/H] rises very rapidly early on (within a few gigayears) in the inner zones, because of the intense star formation and the large contribution of CCSNe. During that period, the Li abundance rises slowly, because a small fraction of novae have the time to release their Li. In contrast, the outer zones reach similar levels of [Fe/H] much later, at a time when a larger fraction of nova have operated. Interestingly, the opposite happens for products of massive stars (alpha-elements), which have very short DTDs, as we discuss below.

In Fig. 7, we have also plotted the corresponding data for main sequence stars from GALAH (see Sect. 2.2 for details on the sample selection) and open clusters (OC, blue asterisks, Romano et al. 2021). With the exception of the highest look-back time (most uncertain ages), the [Fe/H] values of GALAH correspond well to the evolution of Fe in our model zones of 6–10 kpc for the last 6 Gyr and to zones inside 6 kpc for older ages. However, the vast majority of the Li values of GALAH at all ages (and even more so for stars older than 2.5 Gyr) are far below the BBN value, a clear signature of Li depletion in stellar envelopes. This conclusion is independent of the formation zone of the GALAH main sequence stars, which are all observed locally.

In the framework of our model, the four OCs with the highest [Fe/H] and ages 1–2 Gyr are formed in the region of R_GC = 5 to 6 kpc. However, the observed highest Li abundance of those OCs is substantially smaller than the model abundances in those zones and ages: it corresponds to stars formed outside R_GC = 8 kpc. This discrepancy is again attributed to the depletion of Li on the surfaces of those stars, formed in the inner disc, in the last couple of Gyr.

The inserts in Fig. 7 zoom in on the youngest ages 0.01–0.05 Gyr, providing a different picture: the abundances of Fe and Li_max of those clusters correspond to the model predictions for the zones 8–10 kpc. Those very young OCs were formed locally – having little time to migrate – and for the same reason they had little Li depletion in their envelopes.

Figure 8 (top) displays the commonly studied A(Li) versus [Fe/H] diagram for the various model zones. Because of the high initial A(Li), the various curves evolve very closely to each other and start exceeding the primordial value around [Fe/H]~−1, but they display no significant divergence (higher than an observational uncertainty of 0.1 dex) up to the solar value. As discussed before, at a given metallicity, it is the outer zones that have the higher Li abundance in their gas. Above solar metallicity, the rise of A(Li) concerns only the inner zones, because at 8 kpc, gas metallicity stops at [Fe/H]~0.1 dex, (comparable to one of the young local stars and ISM). In those inner zones, the abundance of Li continues increasing, more rapidly than the one of Fe, which is affected by the aforementioned dilution effect from the Fe-poor ejecta of low-mass stars.

Compared to the observations, there is a fair agreement of the various curves (within observational uncertainties) with the open cluster data when [Fe/H]>−0.5 dex. In the highest metallicities, up to [Fe/H]=0.3 dex, the data seem to be well explained by the evolution of the innermost zones (<4 kpc), with no need to assume any Li depletion in the stellar envelopes. However, this interpretation is incorrect. In Fig. 7, the age and [Fe/H] of those super-solar metallicity clusters correspond to their formation at R_GC~ 5 kpc, while the corresponding observed A(Li) is lower than predicted for that zone at that time, and it is thus attributed to depletion in the stellar envelope. This idea was already invoked in Guiglion et al. (2019), who used different data (from AMBRE, lacking stellar ages and hot enough stars) and different models (with initial Li from the Spite plateau, no radial migration and different prescriptions for nova DTD). This interpretation of Li abundances at supersolar metallicities is consistent with the ages of those clusters, around 1–2 Gyr, which is consistent with sizeable Li depletion in their stars. In contrast, no depletion is required in order to interpret the data around solar [Fe/H], since those clusters are very young, less than 0.2 Gyr old (see also Dantas et al. 2022). This is consistent with recent stellar evolution models (e.g. Dumont et al. 2021a,b).

In addition to the open cluster data, Fig. 8 also includes data for field stars from the GALAH survey. These field stars show a larger scatter in A(Li) (top panel) and [Li/Fe] (bottom panel, discussed below) compared to the open clusters due to factors such as large span in age (from 25 Myr to 12.8 Gyr) and T_eff (meaning that the sample also contains G- and K-type stars that deplete Li much faster than their F-type counterparts). While the GALAH data are not directly comparable to the models, they can thus be seen as providing useful observational limits.

The bottom panel of Fig. 8 provides another, original, view, of the Li versus metallicity relation. [Li/Fe] decreases with increasing [Fe/H] as does [α/Fe] (not shown). However, the reason in the case of Li is not the different DTDs of the corresponding sources (short DTD for α-elements from CCSNe versus long DTD for Fe from SNe Ia) but simply the initial abundances (primordial for Li versus zero for Fe, so that [Li/Fe] can only decrease with time). Around solar metallicity, [Li/Fe] flattens because the main sources of the two elements (novae for Li and SN Ia for Fe) have similar DTDs. Finally, at the highest metallicities of the innermost zones [Li/Fe] rises up again because of the differential dilution discussed before (higher for Fe than for Li).

Furthermore, as already mentioned, at a given metallicity [Li/Fe] is higher in the outer zones with their younger (low [α/Fe]) population than in the inner ones with their older stars (high [α/Fe] population); this behaviour is opposite to the one of [α/Fe] versus metallicity: at a given metallicity, the [α/Fe] ratio is lower in an inner than an outer zone, as it is well known from models of the chemical evolution of the Milky Way (Prantzos et al. 2023). This difference, due to the vastly different DTDs of the sources of those elements (very short for α-elements, quite long for Li with the adopted nova ejecta DTDs), is clearly illustrated in Fig. 9, where the [Li/α] ratio is plotted versus [Fe/H]. The model curves are distinguishable now even at low metallicities, in contrast to the situation for [Li/Fe] (bottom panel in Fig. 8).

The observed double-branch behaviour of [α/Fe] versus metallicity (with the high [α/Fe] values at a given [Fe/H] corresponding to the thick disc and the low ones to the thin disc) is interpreted in our model through the combined effects of radial migration, secular evolution and different DTDs for α elements (produced by massive stars) and Fe (mostly produced by SN Ia), as discussed in detail in Kubryk et al. (2015); Prantzos et al. (2023): stars from the inner, rapidly evolving regions (with high [α/Fe] values) are found today in the local disc, which has undergone a slower evolution (leading to lower [α/Fe] values for the same metallicities)². A similar question arises obviously for Li since the DTDs of novae (Li producers) are more extended than those of SN Ia (Fe producers). However, the observational situation is not yet clear. Li abundances in the local thick disc are definitely lower than those of the thin disc, but the upper Li envelope is very poorly defined in the latter case. Starting from the Spite plateau, Ramírez et al. (2012) report a flat A(Li) versus [Fe/H] behaviour up to solar metallicity, while Guiglion et al. (2016) and Fu et al. (2018) find an increase of the upper Li envelope with metallicity and Bensby & Lind (2018); Bensby et al. (2020) find a decrease. The various difficulties plaguing the separation of the discs (by age, kinematics, or chemistry, respectively) and the subsequent evaluation of the upper Li envelopes are summarised in Smiljanic (2020).

From the theoretical side, Prantzos et al. (2017) made the first attempt to quantify the effect with a multi-zone model with radial migration. They found that, if the high Li abundance of BBN is adopted along with a long-lived source of Li (as here), the upper envelope of the thick disc stars should lie above the observations up to a metallicity [Fe/H]~−0.5, and this is also the case here (compare their Fig. 3 with the top panel of our Fig. 8). This implies that up to that metallicity, it is Li depletion that determines the observed upper Li envelope, requiring appropriate stellar physics for its interpretation (Borisov et al. 2024).

In conclusion, the interpretation of the upper envelope of the Li abundance versus metallicity requires, from the observational side, stars that are hot enough (preferentially above the Li-dip, see discussion in Charbonnel et al. 2021) and an estimate of their ages and, from the theoretical side, models of Li depletion in stellar envelopes for various metallicities and multi-zone models of Galactic evolution with radial migration.

Fig. 7 [Fe/H] (top) and A(Li) (bottom) as a function of age for the GALAH DR3 stars (green circles) and open clusters (five-pointed stars colour-coded according to their Galactocentric radius). The solid curves (colour-coded the same way as the open cluster symbols) show the model evolution of A(Li) in different radial zones (4, 5, 6, 8, 10, 12, and 14 kpc, from top to bottom). The insets in both panels focus on the very last 0.5 Gyr of the Galactic evolution. The Sun is formed at Galactocentric radius RGC 6 kpc, 4.56 Gyr ago.

Fig. 8 A(Li) and [Li/Fe] as a function of [Fe/H] (top and bottom panels respectively) for GALAH DR3 stars (circles) and open clusters from Romano et al. (2021) (five-pointed stars). The age of the stars and open clusters is colour-coded. In both panels, the solid curves show model A(Li) in different radial zones (4, 5, 6, 8, 10, 12, 14 kpc from bottom to top; the colours from dark red to orange are as in Fig. 7).

Fig. 9 [Li/α] as a function of [Fe/H] from GALAH, with thin disc stars (blue), thick disc (orange), and halo (purple). We selected the halo stars as those with [Fe/H]≤−1 dex regardless of their [α/Fe], while we distinguished between thin and thick discs on the basis of their [α/Fe] by the Gaussian mixture model. Model curves are colour-coded as in Fig. 8.

Fig. 10 Radial profiles of A(Li), [Fe/H], [α/H], and [O/H] in the galactic gas of our model. Different solid lines show the profiles at different ages (lookback times) and have colours according to their age. The five-pointed stars show the position and abundances of open clusters according to Romano et al. (2021), with age colour-coded. For Berkeley 73 and NGC 2158, we show additional [Fe/H] estimations from Netopil et al. (2016) with smaller five-pointed stars. The grey dashed lines show the solar birth radius and abundances, while the solar symbol shows the present-day location.

4.2 Li across the Galactic disc

Figure 10 shows the predicted radial profiles of A(Li), [Fe/H], [α/H], and [O/H] in the Milky Way gas at different stages of the evolution in our model (every 2 Gyr in age). The final (present-day) gradient of [Fe/H] is found to be d[Fe/H]/dR_GC ~−0.07 dex/kpc in the range 5 < R_GC (kpc) < 15 (Prantzos et al. 2023), which is slightly steeper than values estimated from Cepheid data: −0.06 dex/kpc (Genovali et al. 2014) or − 0.05 dex/kpc (Luck 2018). We notice here that the abundance profile of Fe (or any other element) is not necessarily characterised by a unique slope over the whole radial range. The shape of the abundance profile, flat in the innermost 4 kpc, which is steeply declining in the intermediate region 4 < R_GC (kpc) < 14 and is less steep in the outer disc, resembles the one obtained in Magrini et al. (2009), who adopted very different prescriptions for the star-formation, infall, and SN Ia rates. A break of the slope in the inner disc is also obtained in other studies (e.g. Minchev et al. 2014).

For Li, we find a present-day gradient of d[Li/H]/dR_GC ~ –0.06 dex/kpc, slightly steeper than in Prantzos (2012) who found –0.05 dex/kpc and slightly flatter than in Romano et al. (2021) with –0.07 dex/kpc. Those models, however, used different prescriptions for the nova rate and had no radial migration. We note that in contrast to the Fe gradient, which becomes flatter with time (steeper with age) because of the inside-out disc formation, the Li gradient starts quite flat (because of the initial Li from BBN) and steepens considerably at late times because of the adopted nova DTDs.

The, now indisputable, impact of radial migration on the evolution of the Galactic disc, makes it difficult to infer many of its past properties. This concerns, in particular, the determination of the birth radius of disc stars from observations and therefore the determination of the evolution of the abundance gradient in the ISM (i.e. at the birthplace of stars that are currently observed locally). Minchev et al. (2018) attempted to estimate the birth radii by projecting locally observed mono-age stars in radius along hypothesised age gradients and using the metallicity distributions of those mono-age populations (from the HARPS-GTO survey) to constrain the result. Recently, Zhang et al. (2023) used the result of Minchev et al. (2018) to infer the birthplace of stars and open clusters (OCs) with measured Li abundances from observations with the multi-subject optical fibre facility FLAMES. They found that the highest Li abundances of local and young field stars and OCs show no sign of depletion in the stellar envelopes, while stars born in the inner disc are older on average and are Li-depleted.

Subsequent studies on the evolution of the [Fe/H] gradient at birthplace from local observations showed that the gradient does not evolve monotonically (Lu et al. 2022b) and displays two recent fluctuations, attributed to recent star-formation episodes (Ratcliffe et al. 2023); moreover, the limitations and drawbacks of the method have been underlined already in Minchev et al. (2018) and in more detail in Lu et al. (2022a). We thus refrain from projecting our sample of field stars to their birth radii.

In Fig. 10, we overplot on our model results data for young OCs from Romano et al. (2021). We compute [α/H] as a error-weighted mean of [Mg/H] and [O/H] based on values provided by Van der Swaelmen et al. (2023). For those clusters that do not have measurements of [O/H], we adopt [α/H]=[Mg/H]. We do not show [Fe/H] uncertainties in this figure since they are lower than the size of the symbols.

As has already been discussed, Li is subject to depletion in stellar interiors. With this in mind, we see a relatively good agreement between the data and the model predictions for most OCs. This is the case for the very young OCs in the solar neighbourhood, with [Fe/H], [α/H], and [O/H] values corresponding to the model predictions in the last Gyr of the evolution: their A(Li) is slightly lower than the model prediction, as expected if they had known some Li depletion in their atmosphere. For stars in the inner disc, the agreement between the model and observations is fairly satisfactory for the heavy elements and corresponds to the model metallicity gradient; however, the required amount of Li depletion (taking into account the model predictions) is larger than for the local OCs, as expected from their higher age.

The oldest OC, NGC 2243 with the age of ~4.4 Gyr, is at the same time the most iron-poor and Li-poor one. Its values of [Fe/H], [O/H] and A(Li) correspond well to the radius and age of the Galactic model, with no need for radial migration or Li depletion in its atmosphere. On the other hand, from the stellar evolutionary models that we computed and described in Borisov et al. (2024), we expect Li depletion from the initial abundance of at least ΔA(Li)=0.25–0.30 dex at this age for a star with the composition of the cluster and T_eff range of the stars that were used to compute the OC parameters (6064–6314 K, see Romano et al. 2021 for details). Taking into account the observational uncertainty of ±~0.1 dex in A(Li), we think that the properties of this cluster are well described by the model.

In the outer disc, the situation appears less satisfactory, since all observed abundances are higher than the model predictions. This discrepancy may be due to an underestimate of the rate of evolution in those outer regions because of the model prescriptions. Alternatively, it may be due to the possibility that those OCs have migrated radially from their birthplace in the inner disc where metallicity was higher than in their current position today (see e.g., Dantas et al. 2023, and references therein). Below, we scrutinize each one of the three cases of OCs in the outer disc.

We see a significant overabundance, both in [Fe/H] and α-elements, for the two outer OCs, NGC 2158 (Age_N2158=1.55 Gyr) and Berkeley 73 (Age_B73=1.4 Gyr). The GES values of metallicity are [Fe/H]_N2158=-0.16 and [Fe/H]B73=-0.26 dex, which is ~0.2–0.25 dex higher than the model’s predictions. The difficulty is alleviated if we adopt the lowest values from the literature ([Fe/H]_N2158=-0.32±0.08, Netopil et al. 2016 and $$[\mathrm{F}\mathrm{e}/\mathrm{H}{]}_{\mathrm{B}73}=-{0.41}_{-0.09}^{+0.26}$$, Perren et al. 2022; those values are shown with smaller stars in Fig. 10). For Li, one should consider a higher initial value than the one appearing in the figure, to account for the depletion expected by the cluster age. However, the uncertainties in A(Li), combined with the possibility of some amount of radial displacement (a couple of kpc is reasonably expected during the ~1.5 Gyr ages of those OCs) again alleviate the discrepancy between the model and observations.

Finally, the very young open cluster Berkeley 30 (Age_B30=0.3 Gyr), shows an overabundance of [Fe/H] and [α/H] with respect to the model. This could indicate an overestimation of the metallicity for this cluster or highlight the limitations of the model in accounting for specific local conditions or processes affecting individual clusters. According to the models of Dumont et al. (2021b,a), the expected Li depletion at this age is very small, ΔA(Li)≈0.07 dex. Thus, despite the challenges with [Fe/H] and [α/H], the cluster’s A(Li) is consistent with the model predictions well within the observational uncertainties.

Overall, we think that the position of the studied OCs in the five-dimensional phase space of (R_GC, age, [Fe/H], A(Li) and [O/H]) are well described by the model, especially if on considers all the observational uncertainties related to Li determination, nicely resumed in the discussion of Romano et al. (2021). However, it should also be noted that for the youngest OCs, the observed [Fe/H] values appear flatter than predicted by the model. This discrepancy could suggest that the model overestimates the steepness of the gradient for young populations. However, given the small number of OCs in this regime and the uncertainties involved in our model, we refrain from drawing firm conclusions.

5 Summary

In this work, we have reassessed the evolution of Li in the Galaxy, focussing on the role of novae. Although that source was proposed a long time ago as a key stellar source for Li production, the quantitative evaluation of its importance was hampered by two factors: the very uncertain Li yields of novae and the unknown DTD of those sources. Recent hydrodynamical models do not help with the former problem, as they have found very low Li mass fractions in the nova ejecta and thus insufficient Li yields, more than 20 times lower than required by simple theoretical expectations. This includes the surprising recent results of Gao et al. (2024), who find high Li mass fractions but very low ejecta masses and thus low Li yields. In contrast, observations of ⁷Be abundances in nova explosions in the past ten years have found sufficiently high values for the Li mass fraction in the nova ejecta. As for the nova DTDs, the recent work of Kemp et al. (2022b) provides a very promising theoretical framework (see Sect. 3.2 and Fig. 1).

We first used a well-tested one-zone model of GCE to explore the impact of two potentially important stellar sources of Li, namely AGB stars and novae. For the former source, we adopted the mass and metallicity dependent yields of stars with mass 1.5–8 M_⊙ from Karakas & Lugaro (2016), who find Li production by hot bottom burning in stars more massive than ~4 M_⊙. Our analysis shows that AGB stars fall short (by a large factor) of reproducing the Li meteoritic abundance at the Sun’s formation, 4.5 Gyr ago. Unless the Li production is seriously underestimated in current AGB calculations, those stars cannot be considered as an important Li source, as already discussed in several analogous studies in the past (Travaglio et al. 2001; Romano et al. 2001; Prantzos 2012).

In contrast, our investigation of Li production from novae has yielded interesting findings. Our simple one-zone model shows that by using the Kemp et al. (2022b) DTDs and within the bounds of Li mass fraction dispersion from nova observations, nova can reproduce the Li meteoritic abundance. This confirms previous findings using similar Li yields but different nova rates and different one-zone models of GCE (Cescutti & Molaro 2019; Grisoni et al. 2019) and, of course, the results of Kemp et al. (2022a) with the same nova DTDs. In the framework of our one-zone model, we find that the typical mass fraction in nova ejecta should be X_Li ~ 10⁻⁴, which is rather on the upper range of current nova observations.

We found that the primordial contribution from BBN dominates the total Li abundance until a metallicity of [Fe/H]~−0.5. The relative contributions to the proto-solar (meteoritic) Li are found to be approximately 28% from BBN, 17% from GCR, and 55% from nova. Although the stellar component dominates the abundance of Li at the Sun’s formation, the other two components almost match it. Presently, however, the stellar contribution of novae clearly dominates, having made approximately two-thirds of the Li present in the Galaxy.

We then extended our analysis by using a multi-zone model, presented in Prantzos et al. (2023), that takes into account the radial migration of stars and various long-lived nucleosynthesis sources, including novae. A crucial test regarding the nova DTDs of Kemp et al. (2022b) was provided by the current nova rate, which was found to be R_Novae ~ 36/yr in our model and is in fair agreement with current observations.

In the framework of our multi-zone model, we find that the Sun was born about 2 kpc inwards from its current position, namely at galactocentric radius R_GC~6 kpc. The nova rate in that zone includes novae with locally formed progenitors but also novae with progenitors that migrated from other Galactic regions. As a result, the mass fraction of Li in nova ejecta required to produce the meteoritic Li at R_GC=6 kpc and t_lookback=4.56 Gyr ago was found to be ~7×10⁻⁵, close to the average range of observationally inferred values in Molaro et al. (2023).

We found large differences in Fe and Li abundances between the various zones at any given time (Fig. 7). However, for a given metallicity, Li abundances vary little between the various zones (Fig. 8) and are at the limits of observational uncertainties. This is due to the interplay between various factors: high early star-formation activity in the inner zones versus low activity in the outer ones and different DTDs for the Fe ejected by SN Ia and the Li ejected by novae, with the latter DTD being more extended in time than the former. This is at the origin of a counterintuitive result first obtained in Prantzos et al. (2017), namely that at the same metallicity, Li is more abundant in the outer zones than in the inner ones if its source is sufficiently long-lived.

We found that the most metal-rich OCs, with metallicities [Fe/H]=0.1–0.3 dex, are formed in the region 5–6 kpc (Fig. 7 top) and that their stars have suffered modest Li depletion of 0.2–0.3 dex during the 2–4 Gyr of their evolution. In general, taking into account the observational uncertainties, we find that the position of the studied OCs in the five-dimensional phase space of (R_GC, age, [Fe/H], A(Li), and [O/H]) are well described by the model.

Despite the success of the theoretical framework for nova provided recently by Kemp (2022), it remains to be seen how Li observations can help constrain the proposed nova ejecta DTDs, which differ substantially from the nova event DTDs. The latter are constrained by the observed Galactic nova rate, but the former require good knowledge of the upper envelope of Li observations at various radii and good understanding of the Li depletion in stars of various masses, ages, and metallicities. Current (APOGEE, GALAH, RAVE, LAMOST) and forthcoming (4MOST, WEAVE) surveys will contribute to the observational side (see Nepal et al. 2023, and references therein). In our recent work on Li depletion in stars of halo metallicities (Borisov et al. 2024), we found excellent agreement between model predictions and observed Li trends in post-turnoff stars in the globular cluster NGC 6752. This opens perspectives for refining the estimates of initial Li abundance in metal-rich globular clusters (e.g. 47 Tuc and NGC 6496), thus providing additional constraints on the intermediate metallicity range between Population II and Population I stars. This integrated approach will lead to a deeper understanding of Li evolution from primordial nucleosynthesis to present-day observations.

Acknowledgements

We are grateful to Dr. Alex Kemp for providing us with the set of his metallicity-dependent DTDs of nova events and ejecta; also to Jordi José, Paolo Molaro and Carlos Abia, as well as to the referee Gabriele Cescutti, for a careful reading of the manuscript and very useful comments and suggestions that improved considerably the final text. We acknowledge support from the Swiss National Science Foundation (SNF; Project 200021-212160, PI CC).

References

All Figures

Fig. 1 Derived delay time distributionss for novae event rates $\dot{N}$ (top) and ejecta mass rates $\dot{M}$ (second from top) per 1 M⊙ of stars formed and for a few selected metallicities, adopted from Kemp et al. (2022b). The slopes X of power-law fits appear colour-coded for each metallicity. The grey curves display the corresponding DTDs for SN Ia, with a single slope x=-1.1 above t=0.3 Gyr; the SN Ia curves are multiplied by 10 000 in the top panel and divided by 100 in the second panel. Third panel: average ejecta mass of Kemp et al. (2022a) compared to the highest ejecta mass of the models of Starrfield et al. (2024) indicated by upper limits, to hydrodynamical models of José et al. (2020) within grey shaded areas, and to available observations from S07 (Schwarz et al. 2007), M22 (Molaro et al. 2022), P22 (Pandey et al. 2022), and DM15 (Das & Mondal 2015). Bottom panel: average Li yield when assuming the ejecta mass of the third panel and an average Li mass fraction ⟨XLi⟩OBS from observations (see Fig. 2 and text). The grey arrows indicate the highest Li yield obtained in any of the cited studies.

In the text

Fig. 2 Relative abundances of Li/H in nova explosions. The shaded area shows the range of observed values and the red horizontal line is the average value, as reported in Molaro et al. (2022). Curves represent the highest values found in the calculations of Starrfield et al. (2024, in red), Rukeya et al. (2017, in green) and Gao et al. (2024, in blue), in explosions of CO and ONe white dwarfs of various masses. The results of José et al. (2020) are considerably below the scale of the figure.

In the text

	Fig. 3 Average Li yields of AGB stars in the range 1.5 to Mmax as a function of Mmax (according to Eq. (3)) from Karakas & Lugaro (2016). The data are colour-coded according to the initial stellar metallicity. They are compared to the average yield required to obtain the meteoritic Li value at [Fe/H]=0 (dashed curve).
In the text

Fig. 4 Li abundance as a function of [Fe/H] with nova as a sole stellar source and with different mass fractions of Li in the nova ejecta. Our ‘best-fit’ value (red dashed line) is slightly higher than the mean value of the observational analysis of Molaro et al. (2022) (red solid line), and the pinkish shaded area corresponds to the dispersion of XLi around that mean value. The theoretical nova Li yields of José & Hernanz (1998) and Starrfield et al. (2024) are also used (purple and cyan curves, respectively). The five-pointed stars show the abundance of the open clusters, while the underlying number density plot represents the GALAH sample. The solar symbol indicates the meteoritic value.

In the text

Fig. 5 Evolution of Li abundance (top) and relative contributions (by mass) of different Li sources (bottom) as a function of [Fe/H]. In the top panel, the evolution of total abundance is shown in black, and the solar symbols indicate the solar meteoritic values for 6Li and 7Li. The other curves (in both panels) show the components from the BBN, novae, and GCR (two curves for GCR show components both for 6Li and 7Li). In the bottom panel, the solar symbols indicate the contributions at solar metallicity.

In the text

	Fig. 6 Evolution of total SFR (orange), as well as rates of novae (black), SNe Ia (red), and CCSNe (blue) in the Milky Way predicted by our model. In the legend, the final model values (at 12 Gyr) are indicated. The SFR values are on the right axis (linear scale). Points with error bars show the present-day estimated values (see Sect. 4.1 for references).
In the text

Fig. 7 [Fe/H] (top) and A(Li) (bottom) as a function of age for the GALAH DR3 stars (green circles) and open clusters (five-pointed stars colour-coded according to their Galactocentric radius). The solid curves (colour-coded the same way as the open cluster symbols) show the model evolution of A(Li) in different radial zones (4, 5, 6, 8, 10, 12, and 14 kpc, from top to bottom). The insets in both panels focus on the very last 0.5 Gyr of the Galactic evolution. The Sun is formed at Galactocentric radius RGC 6 kpc, 4.56 Gyr ago.

In the text

Fig. 8 A(Li) and [Li/Fe] as a function of [Fe/H] (top and bottom panels respectively) for GALAH DR3 stars (circles) and open clusters from Romano et al. (2021) (five-pointed stars). The age of the stars and open clusters is colour-coded. In both panels, the solid curves show model A(Li) in different radial zones (4, 5, 6, 8, 10, 12, 14 kpc from bottom to top; the colours from dark red to orange are as in Fig. 7).

In the text

	Fig. 9 [Li/α] as a function of [Fe/H] from GALAH, with thin disc stars (blue), thick disc (orange), and halo (purple). We selected the halo stars as those with [Fe/H]≤−1 dex regardless of their [α/Fe], while we distinguished between thin and thick discs on the basis of their [α/Fe] by the Gaussian mixture model. Model curves are colour-coded as in Fig. 8.
In the text

Fig. 10 Radial profiles of A(Li), [Fe/H], [α/H], and [O/H] in the galactic gas of our model. Different solid lines show the profiles at different ages (lookback times) and have colours according to their age. The five-pointed stars show the position and abundances of open clusters according to Romano et al. (2021), with age colour-coded. For Berkeley 73 and NGC 2158, we show additional [Fe/H] estimations from Netopil et al. (2016) with smaller five-pointed stars. The grey dashed lines show the solar birth radius and abundances, while the solar symbol shows the present-day location.

In the text