© The Authors 2023

1 Introduction

The chemical evolution of galaxies describes the changes in the composition of the interstellar medium (ISM) produced by subsequent generations of stars. In this context, the analytic models are a powerful tool used to predict, among other phenomena, the evolution of the metallicity and the production of chemical elements on short timescales in different galactic systems.

Generally, analytic solutions for the chemical evolution of galaxies have been presented for the so-called ‘simple model’ (Tinsley 1980), which assumes: (i) an initial mass function (IMF) not dependent on time, (ii) the gas is well mixed at any time of the galaxy evolution (instantaneous mixing approximation, IMA), and (iii) the lifetime of stars with mass m ≥ 1 M_⊙ is negligible compared to the timescale of the stars with mass m < 1 M_⊙, whose longer lifetimes are behind the approximation of no contribution to the chemical enrichment of the ISM (see also the discussion in Matteucci 2012, 2021). This assumption constitutes the instantaneous recycling approximation (IRA). In this framework, several analytic solutions for the evolution of the gas phase metallicity have been presented adding more complexity to the system: namely, the infall of gas, galactic winds, radial gas flows, and interactions between galaxies and galactic fountains (Matteucci & Chiosi 1983; Lacey & Fall 1985; Clayton 1988; Edmunds 1990; Recchi et al. 2008; Spitoni et al. 2010; Spitoni 2015; Lilly et al. 2013; Peng et al. 2015; Kudritzki et al. 2015, 2021). Under the IRA approximation, Spitoni et al. (2017) presented the first analytic solution for the time evolution of the metallicity, gas mass fraction, and total mass assuming an exponential infall rate of gas, although no predictions were made for the iron abundance. Beverage et al. (2021); Spitoni et al. (2020, 2021a) used this analytic solution to model the properties of the star forming and quenched galaxies. One of the major limitations of the models described above is the impossibility of obtaining realistic predictions for those elements produced on longer timescales, such as iron.

Supernovae (SNe) of Type Ia are considered to be major producers of ⁵⁶Fe, although a smaller fraction of this element is produced by core-collapse SNe. This was first demonstrated by Greggio & Renzini (1983) and then by Matteucci & Greggio (1986), who showed that with a normal IMF suitable for the Milky Way, the Type Ia SNe contribute up to 70% of the ⁵⁶Fe enrichment in the solar vicinity. Therefore, understanding Type Ia SN progenitors is key for modelling iron production.

Historically, two main channels have been proposed for the formation of the Type Ia SNe: (i) the single-degenerate (Whelan & Iben 1973, SD) and (ii) the double-degenerate (Iben & Tutukov 1994, DD) scenarios. In the first channel, the primary, intermediate mass star of a binary system evolves to produce a carbon oxygen white dwarf (WD) with a close companion. When the secondary star evolves, it fills its Roche Lobe promoting accretion on to the WD, which grows in mass and can explode after reaching the Chandrasekhar limit (roughly 1.4 M_⊙). In the DD channel, when the secondary star fills its Roche Lobe, the companion WD does not accrete the incoming material, forming instead a common envelope engulfing the two stars. This common envelope is eventually lost from the system, leaving behind a close DD system. This system looses orbital energy by emitting gravitational waves (Lorén-Aguilar et al. 2005), leading to the final merging of the two WDs. If the total mass of the system exceeds the Chandrasekhar limit, explosion can occur. Both channels imply a time delay between the formation of the progenitors and the final explosion that can be much larger than lifetime of massive stars. Thus, the instantaneous recycling approximation is no longer accurate for those elements produced mainly in Type Ia SNe, such as iron, requiring a more advance approach that accounts for the distribution of the mentioned time delay.

In the last years, more channels involving sub-Chandrasekhar mass have been proposed for Type Ia SNe to explain some peculiar cases (Nomoto 1982; Iben & Tutukov 1991; Pakmor et al. 2012). These alternatives, however, produce negligible effects on the ⁵⁶Fe production (see the extensive discussion in Palla 2021 and references therein).

By relaxing the IRA approximation (Chiosi & Matteucci 1982), the detailed chemical evolution model originally proposed by Chiappini et al. (1997) can trace the dichotomy in the [α/Fe] versus [Fe/H] diagram (the so-called low- and high-α sequences) observed in the Galactic disc (Lee et al. 2011; Haywood et al. 2013, 2015; Recio-Blanco et al. 2014; Anders et al. 2014; Nidever et al. 2014; Hayden et al. 2014, 2015; Bovy et al. 2016; Gaia Collaboration 2023) and subsequently explained by Kobayashi et al. (1998, Kobayashi et al. 2006); Fenner et al. (2002); Noguchi (2018); Spitoni et al. (2019b, 2023); Lian et al. (2020). These works assume that the Galaxy has been formed by one or more separated accretion episodes, modelled by decaying-time exponential-like infalls of gas.

Vincenzo et al. (2017) have provided a numeric solution for chemical evolution of the Galactic disc by assuming the IRA for chemical elements produced by massive stars and the delay time distribution formalism for the iron. They show the results for the solar neighbourhood for the one-infall model for the single degenerate (SD) scenario proposed by Matteucci & Recchi (2001), as well as that considering the bimodal DTD of Mannucci et al. (2006), in which one half of the Type Ia SNe are produced in prompt explosions (~ few Myr), while the remaining show a wider time distribution (~ hundreds of Myr). Vincenzo et al. (2017) found a good agreement between their predictions and the Bensby et al. (2014) data for [O/Fe] and [Si/Fe] versus [Fe/H] assuming the DTD of Matteucci & Recchi (2001).

Weinberg et al. (2017) avoided the numerical integration of the models with DTD by presenting the analytic solution for the evolution of the iron produced in Type Ia SN events. These authors were able to model the evolution of [Fe/H] and [α/Fe] for three different star formation histories (constant, exponentially declining, linear-exponential). Their analysis, however, is restricted to a specific prescription for the delay time distribution of Type Ia SNe, which decreases exponentially with time. Similarly, Pantoni et al. (2019); Lapi et al. (2020) found the analytic solutions for the evolution of iron considering an exponential DTD.

In this work, we present a new analytic solution for the chemical evolution of galactic systems accounting for the enrichment from Type Ia SNe. Compared to the numerical approach, the analytical solutions have the advantage of providing the exact abundance values at any evolutionary time, with no approximation errors and in a more direct fashion than the recurrent iteration over previous time steps. We consider a prescription for the DTD that extends these used in previous works (Weinberg et al. 2017; Pantoni et al. 2019; Lapi et al. 2020). We test our solution with the DTDs proposed by Matteucci & Recchi (2001), Greggio (2005), Mannucci et al. (2006), Totani et al. (2008), Pritchet et al. (2008) and Strolger et al. (2004, 2005). The paper is organised as follows: in Sect. 2, we present the chemical evolution model, the prescription for the IRA approximation and the adopted formalism for the Type Ia SNe enrichment (detailed in Appendix A). In Sect. 3 and Appendix B, we present the analytic solutions for different DTDs and apply them to the one and two infall scenarios. In Sect. 4, we prove that the new solution can be a handy tool for studies in Galactic archaeology. Using our analytic solution, here we study the chemical dichotomy of the disc and compare it with that observed by APOGEE DR17 (Abdurro’uf et al. 2022). We include similar prescriptions of the detailed two-infall models proposed by Spitoni et al. (2019b, 2021b) which were designed to reproduce APOKASC (Silva Aguirre et al. 2018) and APOGEE DR16 (Ahumada et al. 2020) data, respectively. In Sect. 5, we model the iron and silicon abundances of a Milky Way-like galaxy from its SFR. Our conclusions and plans for future works are summarised in Sect. 6.

2 Model prescriptions

In this section, we present the main model prescriptions. For a more complete discussion of Galactic chemical evolution assumptions and ingredients, we refer to the review by Matteucci (2021) and the book by Matteucci (2012).

2.1 Useful quantities with IRA

Under the assumption of the IRA and the instantaneous mixing approximation, the returned mass fraction, R, indicates the total amount of mass restored to the ISM by a single stellar generation after ~ 10 Gyr (Tinsley 1980). Given an IMF ϕ(m), R can be computed as: $$R=\frac{{\displaystyle {\int }_{1M_{\odot }}^{100M_{\odot }}\left(m-M_{\text{R}}\right)\varphi \left(m\right)\text{d}m}}{{\displaystyle {\int }_{0.1M_{\odot }}^{100M_{\odot }}m\varphi \left(m\right)\text{d}m}},$$(1)

where M_R is the mass of the stellar remnant. Similarly, the yield per stellar generation for the X element, 〈y_X〉, is defined as: $$\langle y_{\text{X}}\rangle =\frac{1}{1-R}{\displaystyle {\int }_{1M_{\odot }}^{100M_{\odot }}m{p}_{\text{X}}\left(m\right)\varphi \left(m\right)\text{d}m},$$(2)

where p_X(m) is the ratio between the ejected mass of the element i, and newly produced by a star of mass m. As Eq. (2) shows, only those stars with masses larger than 1 M_⊙ contribute to the chemical enrichment of the ISM, while the (1 – R) term in the denominator accounts for the amount of mass locked up in stars of lower mass. Thus, 〈y_X〉 can be understood as the ratio between the ejected mass and remnant mass for the X element in a single stellar generation. Both for R and 〈y_X〉, the choice of the lower mass limit in their definitions does not significantly change their values (Tinsley 1980).

Fig. 1 Delay time distribution functions normalised to their own maximum value adopted in this work: MR01 (Matteucci & Recchi 2001, blue line), WIDE and CLOSE G05 (Greggio 2005, purple and pink lines, respectively), MVP06 (Mannucci et al. 2006, green line), S05 (Strolger et al. 2005, yellow line), and T08 (Totani et al. 2008, cyan line).

2.2 Type Ia SNe and the DTD formalism

Greggio (2005) proposed a new formalism for the Type SN Ia rate based on the concept of the delay time distribution, namely the functional form which indicates how the SN Ia progenitors die as a function of time considering an instantaneous starburst, that is, a single stellar population. Given a SFR and a delay time distribution DTD(t), the SN Ia rate at time t is obtained as the following integral: $${ℛ}_{\text{Ia}}\left(t\right)=C_{\text{Ia}}{\displaystyle \underset{{\tau }_1}{\overset{\min \left(t,{\tau }_2\right)}{\int }}{\text{DTD}}_{\text{Ia}}\left(\tau \right)\psi \left(t-\tau \right)\text{d}\tau ,}$$(3)

where τ₁(τ₂) is the minimum (maximum) time for the explosion of a Type Ia SN, and the normalisation constant C_Ia is set to reproduce the observed present time Type Ia SN rate.

In this work, we provide an analytic solution of the chemical evolution model considering different realisations of the DTD (Figs. 1 and 2): the one computed for the single degenerate scenario by Matteucci & Recchi (2001, hereafter MR01), those proposed by Greggio (2005, G05) for the WIDE and CLOSE double degenerate scenario, the DTDs derived empirically by Mannucci et al. (2006, MVP06), Totani et al. (2008, T08), and Pritchet et al. (2008, P08), as well as the Gaussian DTD proposed by Strolger et al. (2005, S05), from the observed cosmic Type Ia SN rate. The terminology for the G05 DTDs is related to the distribution of the separation of the DD systems, which can be more or less populated at the low values, as resulting from respectively a less or more efficient transfer of orbital energy to the potential energy of the envelope. Correspondingly, the distribution of gravitational delays turns out more skewed towards the short delays in the CLOSE DD scheme, leading to steeper DTDs. We refer to Greggio (2005) for a detailed technical description of the two cases.

We consider the ample variety of formulations mentioned above because they all are based on astrophysical arguments. We note that different DTDs can account for the observed Type Ia SN rates in external galaxies, within the current uncertainties (Botticella et al. 2017; Greggio & Cappellaro 2019). Thus, the analytic formulations presented in this work for these DTDs can be used to construct models for the evolution of the iron abundance (or other element produced mainly in Type Ia SNe) in other galaxies besides the Milky Way.

The DTD of MVP06 is described as a combination of a Gaussian and an exponential distribution, leading to the bimodality that characterises this DTD. On the contrary, S05 suggested a single Gaussian distribution with no prompt Type Ia SNe that peaks at 3–4 Gyr. The T08 and P08 DTDs are described by power law relations of the form t⁻¹ and t^−1/2, respectively, in which the T08 traces the slope of the WIDE G05 DTD. For a detailed explanation of the functional forms of the MR01 and G05 DTDs, we refer to Matteucci et al. (2006), Greggio (2005) as well as Sects. 2.1 and 2.2 in Bonaparte et al. (2013).

Figure 2 shows the cumulative fraction of Type Ia SNe from a single stellar population as a function of time. We consider the value of 150 Myr for the definition of the upper limit of the prompt regime (approximately an intermediate value among the time intervals given by Acharova et al. 2022; Aubourg et al. 2008; Maoz & Badenes 2010). As we may note, there are significant differences among the prompt fractions of the DTDs: in the SD scenario, the fraction of prompt Type Ia SNe is ~10% (Bonaparte et al. 2013), while for the WIDE and CLOSE DD G05 scenarios these prompt fractions are ~6% and ~22%, respectively. For the MVP06 DTD, this prompt fraction is substantially higher (50%). The empirical T08 and P08 DTDs have prompt fractions similar to that of the WIDE DD case (approximately 6% and 9%, respectively).

In order to integrate analytically the chemical evolution model for the complex DTDs, we approximate them by a combination of truncated Gaussian, exponential, and inverse of time (~(t – τ₀)⁻¹) functions. Hence, the general expression for the fit of a DTD is: $$\begin{array}{l}\text{DTD}\left(t\right)={\displaystyle \sum _{i=1}^{N_{\text{G}}}{A}_{\text{G,}i}\text{\hspace{0.17em}exp\hspace{0.17em}}\left(-\frac{{\left(t-\tau {\prime }_i\right)}^2}{2\sigma {\prime }_i^2}\right)\cdot \mathrm{𝟙}{}_{\left[{\tau }_{1,\text{G,}i},{\tau }_{2,\text{G,}i}\right]}}\left(t\right)+\hfill \\ +{\displaystyle \sum _{i=1}^{N_{\text{E}}}{A}_{\text{E,}i}\text{\hspace{0.17em}exp}}\left(\frac{-t}{{\tau }_{\text{D,}i}}\right)\cdot \mathrm{𝟙}{}_{{}_{\left[{\tau }_{{}_{1,\text{E,}i}},{\tau }_{{}_{2,\text{E,}i}}\right]}}\left(t\right)+\hfill \\ +{\displaystyle \sum _{i=1}^{N_{\text{I}}}{A}_{\text{I,}i}}\frac{{\tau }_{\text{I}}}{t-{\tau }_0}\cdot \mathrm{𝟙}{}_{{}_{\left[{\tau }_{{}_{1,\text{I,}i}},{\tau }_{{}_{2,\text{I,}i}}\right]}}\left(t\right),\hfill \end{array}$$(4)

where we introduce τ_Ι = 1.00 Gyr just to keep the same units in all the amplitudes A_x,i. N_G, N_E, and N_I are the number of Gaussian, exponential, and (t – τ₀)⁻¹ functions that characterise the DTDs, respectively¹. The value of the indicator function 1 is one if the argument is within the interval of the sub-index and zero otherwise.

We note that within the formalism of Eq. (4), the S05 and T08 DTDs emerge by setting A_E,i, A_I,i = 0 and A_G,i, A_E,i = 0, respectively. Similarly, a proper combination of A_G,i, A_E,i with A_I,i = 0 results in the MVP06 DTD. For those DTDs with a more complex functional form, such as these of MR01, G05 and P05, the values of the amplitudes A_G,i, A_E,i, A_I,i are determined by a restricted least-squares fitting that minimises the difference between the original DTD and Eq. (4) (see Table 1). Since the implementation of this method is rather technical, we refer to Appendix A for a more detailed explanation of this procedure.

Fig. 2 Cumulative fraction of Type Ia SNe as a function of time for the DTDs used in this work: MR01 (Matteucci & Recchi 2001, blue line), WIDE and CLOSE G05 (Greggio 2005, purple and pink lines, respectively), MVP06 (Mannucci et al. 2006, green line), S05 (Strolger et al. 2005, yellow line), and T08(Totani et al. 2008, cyan line). Dashed vertical line denotes the reference time limit for the prompt of 150 Myr.

2.3 Chemical evolution equations

We consider a one-zone chemical evolution model, assuming the following form for the Kennicutt-Schmidt law (Schmidt 1959; Kennicutt 1989) law for the SFR: $$\psi \left(t\right)=v_{\text{L}}{\sigma }_{\text{gas}}\left(t\right),$$(5)

where v_L is the star formation efficiency (SFE) and has the dimension of (Gyr⁻¹). As in Spitoni et al. (2017) and Vincenzo et al. (2017), we consider galactic winds proportional to the SFR: $$W\left(t\right)=\omega \psi \left(t\right).$$(6)

In the scenario proposed by several works in literature (e.g. Chiosi 1980; Boissier & Prantzos 2000; Schönrich & Binney 2009; Andrews et al. 2017), the galaxy has been formed out by the accretion of distinctive exponential infall events. Here, we provide analytic solutions for the chemical evolution of a system built up by N infalls, in which the total gas accretion rate can be expressed as $$I\left(t\right)={\displaystyle \sum _{j=1}^N{A}_j\text{\hspace{0.17em}exp\hspace{0.17em}}\left(-\frac{\text{Δ}{t}_j}{{\tau }_j}\right)\theta \left(\text{Δ}{t}_j\right),}$$(7)

where the jth infall starts at time t_j and is characterised by the timescale, τ_j; while the amplitude, A_j, tunes the amount of gas accreted due to the jth infall. In order to simplify the notation, we denote t – t_j as Δt_j and the Heaviside step function as θ. Vincenzo et al. (2017) provided the following analytic expression for the star formation history of a galactic system formed by the accretion of N separate infalls characterised by exponential rate decays in presence of the IRA and the Schmidt (1959) law for the SFR: $$\begin{array}{l}\psi \left(t\right)=v_{\text{L}}{\displaystyle \sum _{\begin{array}{c}j=1\\ {\tau }_j\ne {\alpha }^{-1}\end{array}}^N\frac{A_j{\tau }_j}{\alpha {\tau }_j-1}\left[\text{exp\hspace{0.17em}}\left(-\text{Δ}{t}_j/{\tau }_j\right)-\text{exp\hspace{0.17em}}\left(-\alpha \text{Δ}{t}_j\right)\right]\theta \left(\text{Δ}{t}_j\right)}\hfill \\ +v_{\text{L}}{\sigma }_{\text{gas}}\left(0\right)\text{\hspace{0.17em}exp\hspace{0.17em}}\left(-\alpha t\right)\theta \left(t\right)\hfill \\ +v_{\text{L}}{\displaystyle \sum _{\begin{array}{c}j=1\\ {\tau }_j={\alpha }^{-1}\end{array}}^N{A}_j\text{\hspace{0.17em}Δ}{t}_j\cdot \text{exp\hspace{0.17em}}\left(-\alpha \text{Δ}{t}_j\right)\cdot \theta \left(\text{Δ}{t}_j\right),}\hfill \end{array}$$(8)

where α = v_L(1 + ω − R), as indicated in Table 2. Compared to Eq. (13) in Vincenzo et al. (2017), we include the additional summation term in ~ A_jΔt_j that corresponds to the particular case in which τ_j → α⁻¹. Although this term is necessary to provide the full general solution, we can ignore it hereafter since we do not make use of any τ_j = α⁻¹ in this work. The contribution of this term, however, can be found in Appendix B. Finally, the equation for the evolution of the surface gas density for the X-element σ_X(t) with the Type Ia SN contribution is expressed as: $$\frac{\text{d}{\sigma }_{\text{X}}\left(t\right)}{\text{d}t}=\stackrel{\text{IRA\hspace{0.17em}terms}}{\overbrace{-\alpha {\sigma }_{\text{X}}\left(t\right)+\langle y_{\text{X}}\rangle \left(1-R\right)\psi \left(t\right)}}+\stackrel{\text{Type\hspace{0.17em}Ia\hspace{0.17em}SNe\hspace{0.17em}term}}{\overbrace{\langle m_{X\text{,Ia}}\rangle {ℛ}_{\text{Ia}}\left(t\right)}},$$(9)

where 〈m_X,Ia〉 is the mass of the element X synthesized by each single Type Ia SN explosion.

This is a first order inhomogeneous differential equation whose solution can be computed analytically if the SFR ψ(t) takes the form given by Eq. (8). This solution is summarised in Appendix B.

3 Results

In Sect. 3.1, we present the new analytic solutions for the temporal evolution of different chemical elements considering the iron produced by different DTD prescriptions. We also test the effects of different DTD prescriptions for the one-infall scenario on the [α/Fe] versus [Fe/H] abundance ratios (Sect. 3.2) and on the metallicity distribution functions (Sect. 3.3). We refer to Appendix B for the detailed explanation of the analytic form of the solutions and to the ChEAP² (Chemical Evolution Analytic Package) repository for its implementation in the PYTHON language.

3.1 The new analytic solution

In Appendix B, we present the analytic expression for the Type Ia SN rates R_Ia(t) and the surface mass density of the element X, which can be written as the sum of the contribution of the IRA and Type Ia SNe enrichment as: $${\sigma }_{\text{X}}\left(t\right)={\sigma }_{X,\text{IRA}}\left(t\right)+{\sigma }_{X,\text{Ia}}\left(t\right).$$(10)

Similarly, σ_X,Ia (t) can be separated into terms that depend on the Gaussian (σ_X,Ia,G), exponential (σ_X,Ia,E), and inverse time (σ_X,Ia,I) DTDs: $${\sigma }_{X,Ia}\left(t\right)={\displaystyle \sum _{}^{N_{\text{G}}}{\sigma }_{X,Ia,G}\left(t\right)+}{\displaystyle \sum _{}^{N_{\text{E}}}{\sigma }_{X,Ia,E}\left(t\right)+{\displaystyle \sum _{}^{N_{\text{I}}}{\sigma }_{X,Ia,I}\left(t\right).}}$$(11)

Table 2 summarises all the parameters considered in the proposed chemical evolution model, distinguishing between the “galaxy model”, “IRA”, “type Ia SNe & DTD” quantities. Furthermore, we provide some useful definitions to simplify the analytic expressions.

Once σ_X(t) is known, we compute the abundance ratio [X/Fe] as: $$\left[X/\text{Fe}\right]={\text{log}}_{10}\left(\frac{{\sigma }_{\text{X}}}{{\sigma }_{\text{Fe}}}\right)-{\text{SV}}_{\text{Fe}}^{\text{X}},$$(12)

where ${\text{SV}}_{\text{Fe}}^{\text{X}}$ is a scaling factor derived from the solar reference values of Asplund et al. (2009). For the particular case of the iron, its abundance is computed as: $$\left[\text{Fe}/\text{H}\right]={\text{log}}_{10}\left(\frac{{\sigma }_{\text{Fe}}}{0.75\cdot {\sigma }_{\text{gas}}}\right)-{\text{SV}}_{\text{H}}^{\text{Fe}},$$(13)

where σ_gas is given by Eqs. (5) through (8). Based on Big Bang nucleosynthesis, we assume the hydrogen comprises the 75% of the gas mass (factor 0.75 in Eq. (13)). Also, we assume the H abundance in mass does not change significantly during the Galactic evolution.

3.2 Testing the new analytic solutions: effects of the DTD on the one-infall model

In this section, we show the effects of different DTD prescriptions on the simplest case of the one-infall scenario at t₁ = 0 Gyr. This model has been widely used in the past to describe the thin disc of our Galaxy (Spitoni et al. 2015, 2019a; Grisoni et al. 2017, 2018). As in Vincenzo et al. (2017), we study the evolution of oxygen, silicon, and iron by assuming 〈y_O〉 = 1.022 × 10⁻², 〈y_Si〉 = 8.5 × 10⁻⁴, 〈y_Fe〉 = 5.6 × 10⁻⁴ and the returned fraction, R = 0.285. These values are derived from the Kroupa et al. (1993) IMF and the collection of nucleosynthesis yields suggested by Romano et al. (2010, see also Vincenzo et al. 2016). For the stellar yields of Type Ia SNe, we made use of those from Iwamoto et al. (1999).

We consider the same set of parameters as in Vincenzo et al. (2017). This implies a total surface mass density in the solar neighbourhood of σ_tot = 54 M_⊙ pc⁻² and also used widely in other works (e.g. Spitoni et al. 2015; Spitoni & Matteucci 2011), a star formation efficiency of v_L = 2 Gyr⁻¹, an infall timescale for the gas mass accretion of τ₁ = 7 Gyr, and a mass-loading factor of ω = 0.4. In contrast to Vincenzo et al. (2017), we set the value of C_Ia (see Eq. (3)) by imposing the predicted present-day SN Type Ia rate of 0.54 ± 0.12 events per century (Li et al. 2011) in the disc region (3 kpc < R < 20 kpc). Thus, we assume the solar neighbourhood is representative of that annular region. No initial amount of gas is assumed for the Milky Way (σ_gas = 0 M_⊙ pc⁻²Gyr⁻¹).

Figure 3 shows the temporal evolution of Type Ia SN rates for the seven DTDs considered in this study. The rate R_Ia(t) computed with the S05 DTD peaks at later evolution time with respect to the majority of the other distributions because of the lack of prompt Type Ia SNe. However, the variation of R_Ia among the DTDs does not depend only on the shape of the DTD but also on its “convolution” with the star formation history (see Eq. (3)).

In Fig. 4, we show the predicted [O/Fe] (left panel) and [Si/Fe] (right panel) versus [Fe/H] for the seven DTD prescriptions. As already noted in Matteucci et al. (2009), the DTD of S05 shows the longest plateau in the [α/Fe] versus [Fe/H] diagram because no prompt Type Ia SNe is present. The predicted abundance ratios agree with the Matteucci et al. (2009) and Vincenzo et al. (2017) results. We have checked that the fraction σ_X/σ_Fe predicted by our solution, for a generic chemical element X, is equal to the ratio of the yields 〈y_X〉/〈y_Fe〉 in the limit t → 0 Gyr.

Fig. 3 Predicted Type Ia SN rates for the different DTDs considered in this work (see Fig. 1) assuming the one-infall model introduced in Sect. 3.2. The area in the vertical axis refers to that of the Galactic disc region (3 kpc < R < 20 kpc). The rates are expressed in SN events per century and normalized (through the CIa) to reproduce the present time Type Ia SN rate in the Milky Way (Li et al. 2011, filled red circle and error bar).

3.3 The metallicity distribution function

The proposed analytic solution can be used to predict the metallicity distribution function³ (MDF). Given the star formation rate ψ(t), the amount of stellar mass dM formed in an area dA in the time interval [t, t + dt] is $$\text{d}M=\psi \left(t\right)\cdot \text{d}A\cdot \text{d}t.$$(14)

By definition of initial mass function ϕ(m), the total number of stars with masses m_*,min ≤ m < m_* resulting from a single formation event is $$N_{\star }=k{\displaystyle {\int }_{m_{\star ,\min }}^{m_{\star }}\varphi \left(m\right)\text{d}m},$$(15)

where the value of the constant k is determined by the total mass of the population, M_*,tot, as: $$k=\frac{M_{\star ,\text{tot}}}{{\displaystyle {\int }_{m_{\star ,\min }}^{m_{\star ,\max }}m\cdot \varphi \left(m\right)\text{d}m}}.$$(16)

We note that the integration limits in Eq. (16) includes all the possible stellar masses, without excluding sources more massive than m_* as in Eq. (15). Thus, for a population of mass, M_*,tot = dM, the total number of stars within m_*,min < m < m_* per unit of area, dN_*/dA, is $$\frac{\text{d}{N}_{\star }}{\text{d}A}=\psi \left(t\right)\frac{{\displaystyle {\int }_{m_{\star \text{,min}}}^{m_{\star }}\varphi \left(m\right)\text{d}m}}{{\displaystyle {\int }_{m_{\star ,\min }}^{m_{\star ,\max }}m\cdot \varphi \left(m\right)\text{d}m}}\cdot \text{d}t=f\cdot \psi \left(t\right)\cdot \text{d}t.$$(17)

Since the IRA approximation assumes all stars more massive than the Sun die immediately, we must consider m_* = 1 M_⊙ in the numerator of Eq. (17). For m_*,min and m_*,max we use the values 0.1 M_⊙ and 100 M_⊙, respectively, to be consistent with the R and 〈y_X〉 values adopted in this study.

We can construct the MDF from Eq. (17) by integrating dN_* within the limits of each bin in metallicity. This procedure, however, generally requires solving transcendental equations to get the integration limits as a function of the metallicity. This is specially complicate when several infalls are included, since one may need to account for multiple branches of t([Fe/H]). A more practical approach is performed thanks to the following numerical integration $$\begin{array}{l}{\text{MDF}}_i={\displaystyle {\int }_{\text{A}}{\displaystyle {\int }_0^{t_{\text{G}}}\text{d}{N}_{\star }}}\hfill \\ =f\cdot {\displaystyle {\int }_{\text{A}}{\displaystyle {\int }_0^{t_{\text{G}}}\psi \left(t\right)\cdot {\mathrm{𝟙}}_{[{\left[\text{M}/\text{H}\right]}_i,{\left[\text{M}/\text{H}\right]}_{i+1})}\left(\left[\text{M}/\text{H}\right]\left(t\right)\right)\cdot \text{d}t\cdot \text{d}A},}\hfill \end{array}$$(18)

where the left-hand term is the number of stars in the ith metal-licity bin [M/H]_i ≤ [M/H](t) < [M/H]_i+1, A is the integration area and t_G = 13.8 Gyr. Since Eq. (9) has no spatial dependence, we can substitute the integral in dA by the total area, A.

In Fig. 5, we show the MDFs predicted by the one-infall model for the different DTDs. We can see the MR01, CLOSE G05, MVP06, and T08 DTDs results in similar MDFs, while the S05 DTDs shows a peak at [Fe/H] ≈ −0.6 dex. The MDFs computed with the WIDE G05 and P08 DTDs peak at lower metallicities (−0.06 dex and −0.21, respectively), the latter showing a wider distribution in [Fe/H]. In order to explain these discrepancies, we explore the age-metallicity relation and the evolution of the fraction of iron produced by Type Ia SNe (Figs. 6 and 7, respectively). As can be seen in Fig. 6, the CLOSE G05, MVP06, MR01 and T08 DTDs increase the metallicity up to solar values during the first ~2 Gyr to continue afterwards with in a more steady evolution up to [Fe/H] ≈ 0. 25 dex. As Fig. 7 shows, within the initial 1.4 Gyr, the Type Ia SN explosion becomes the dominant iron producing mechanism for the mentioned DTDs, especially for the CLOSE G05 DTD. On the contrary, this transition occurs 2 Gyr later for the S05 and P08 DTDs. The WIDE G05 corresponds to an intermediate case, showing a more quenched iron production after the first ~2 Gyr. The age-metallicity relations for the WIDE G05 and P08 DTDs can explain the peaks at lower metallicities in their MDFs: since in these scenarios the synthesis of iron is slower, most of stars are formed at lower metallicities compared to the other DTDs, requiring more time to reach the plateau value.

The S05 DTD presents the more complex age-metallicity relation, with three different regimes: during the first 2.5 Gyr the IRA mechanism drives the production of iron up to [Fe/H] ≈ −0.6 dex. At this metallicity, the [Fe/H] versus age curve flattens contributing to the peak observed in the MDF. At later times, the Type Ia SNe accelerate the synthesis of iron during the next 3.5 Gyr to the saturation at [Fe/H] ≈ 0.25 dex, showing a more extended plateau compared to the other DTDs.

Fig. 4 Evolution of [O/Fe] (left panel) and [Si/Fe] (right panel) versus [Fe/H] predicted by the one-infall model for seven different DTDs (solid curves). See Sect. 3.2 for the details of the model.

Fig. 5 Normalised MDFs predicted by the one-infall model (see Sect. 3.2) for different DTDs. Colour code is the same as in Fig. 4.

4 Galactic archaeology with the analytic solution: Disc bimodality in the chemical space

From the chemical point of view, the Galactic disc shows two substructures in the [α/Fe] versus [Fe/H] plane: the so-called high-α sequence, classically associated with an old population of stars (thick disc), and the low-α sequence, characterised by the younger stars of the thin disc (Fuhrmann 2004; Reddy et al. 2006; Bensby et al. 2014; Lee et al. 2011; Haywood et al. 2013; Adibekyan et al. 2013). This dichotomy has been confirmed by the analysis of APOGEE data (Nidever et al. 2014; Hayden et al. 2015; Ahumada et al. 2020; Queiroz et al. 2020; Abdurro’uf et al. 2022), the Gaia-ESO survey (e.g. Recio-Blanco et al. 2014; Rojas-Arriagada et al. 2016, Rojas-Arriagada et al. 2017), AMBRE (Mikolaitis et al. 2017; Santos-Peral et al. 2021), GALAH (Buder et al. 2019, 2021), LAMOST (Yu et al. 2021), and Gaia DR3 (Gaia Collaboration 2023; Recio-Blanco et al. 2023). Spitoni et al. (2019b, 2021b) revised the classical two-infall chemical evolution model (Chiappini et al. 1997) to reproduce the trends in the [α/Fe] versus [Fe/H] diagram observed in the APOKASC (Silva Aguirre et al. 2018) and APOGEE DR16 (Ahumada et al. 2020) samples, respectively. In this Galactic formation scenario, the thick and thin disc components have been formed by to distinct episodes of gas accretion. These authors suggest the presence of a ~4 Gyr time delay between these episodes in order to reproduce the high- and low-α sequences, imposing precise asteroseismic ages as a constraint.

Here, we apply the new analytic solution introduced in Sect. 3 in the framework of the two-infall model in order to reproduce the new APOGEE DR17 data (Abdurro’uf et al. 2022) for the abundance ratio [Si/Fe] versus [Fe/H] in the annular region 7.2 kpc< R < 9.2 kpc (i.e. R_⊙ ± 1 kpc). As in Spitoni et al. (2021b), we impose a signal-to-noise ratio S/N > 80, a surface gravity log g < 3.5 and vertical height |Z| < 1 kpc in our selection.

Using the same formalism and notation introduced in Sect. 2.3, this infall rate can be written as: $$ℐ\left(t\right)=\stackrel{\text{1st\hspace{0.17em}infall,\hspace{0.17em}high-}\alpha }{\overbrace{A_1\theta \left(t\right){\text{e}}^{-t/{\tau }_1}}}+\stackrel{\text{2nd\hspace{0.17em}infall,\hspace{0.17em}low-}\alpha }{\overbrace{A_2\theta \left(\text{Δ}{t}_2\right){\text{e}}^{-\text{Δ}{t}_2/{\tau }_2}}}.$$(19)

We impose the present total surface mass density (sum of high- and low α sequence contributions) of σ_tot(tG) = 47.1± 3.4 M_⊙ pc⁻², suggested by McKee et al. (2015) for the local disc.

Initially, we evaluate our analytic solution for the two infall model by assuming the MR01 scenario for the DTD, in which the parameters of the infall are adapted to mimic the observed [Si/Fe] versus [Fe/H]. For the MR01 DTD, these parameters are τ₁ = 0.4 Gyr, τ₂ = 7.0 Gyr, t₁ = 0 Gyr, t₂ = 3 Gyr, A₁ ≈ 35.128 M_⊙ pc⁻² Gyr⁻¹ and A₂ ≈ 10.207 M_⊙ pc⁻²Gyr⁻¹. This combination of A₁ and A₂ implies a second infall four times more massive than the first one. The star formation efficiency v_L is set to 0.75 Gyr⁻¹ and the loading factor for the wind is ω = 0.8.

As the left panel of Fig. 8 illustrates, we can recover the characteristic “loop” in the low-α sequence already found in the detailed chemical evolution models with delayed gas infalls of Calura & Menci (2009); Spitoni et al. (2019b); Palla et al. (2020); Romano et al. (2020); Cescutti et al. (2022). This delayed infall creates the low-α sequence by bringing pristine metal-poor gas into the system, which dilutes the metallicity of interstellar medium while keeping [α/Fe] abundance almost unchanged. When star formation resumes, the Type II SNe produce a steep increment in the [α/Fe] ratio. At later times, the pollution from the Type Ia SNe raises the metallicity and decreases [α/Fe]. This sequence creates a loop in the [α/Fe] versus [Fe/H] diagram that overlaps with the region spanned by the APOGEE DR17 (Abdurro’uf et al. 2022) data.

It is important to underline that in the Spitoni et al. (2021b) model no Galactic winds have been considered to fit the APOGEE DR16 data. On the contrary, we impose a significant mass loss due to the Galactic winds (ω = 0.8) to reproduce the APOGEE DR17 data. Possibly it is due to the nucleosynthetic prescriptions for massive stars used in that work, which, in line with François et al. (2004), include a modification of the Woosley & Weaver (1995) yields in order to mimic the data available in the solar vicinity. Motivated by this explanation, we multiply 〈y_Si〉 and 〈m_Si,Ia〉 by a factor of 0.85 to reproduce the [Si/Fe] versus [Fe/H] diagram observed with APOGEE data (Fig. 9). Using these yields, we can model the chemical evolution track with a lower wind-loading factor (ω = 0.2), while for the rest of the parameters we consider τ₁ = 0.13 Gyr, v_L = 0.8 Gyr⁻¹, τ₂ = 6.75 Gyr , t₂ = 3.5 Gyr, A₁ ≈ 52.886 M_⊙ pc⁻² Gyr⁻¹, and A₂ ≈ 6.898 M_⊙ pc⁻² Gyr⁻¹ (mass ratio between the two infalls of 5.3).

In Fig. 8, we show the age-metallicity relation predicted by our model, where the effects of the dilution produced by the delayed infall is clear. The consequent gap in the star formation rate has a significant effect on the Type Ia SN rate (local minimum at age ≈ 10.5 − 10.0 Gyr in the lower right panel in Fig. 8). A similar feature in the age-metallicity relation has been found by Nissen et al. (2020) in the analysis of the HARPS spectra of local solar-like stars. They note that the distribution of stars in the age-metallicity relation has two distinct populations with a clear age dissection. The authors suggest these two sequences may be interpreted as an evidence of two gas accretion episodes onto the Galactic disc, with a quenched star formation between them. This is in agreement with the scenario proposed by Spitoni et al. (2019b) and with the results shown here. By analysing subgiant stars of LAMOST, Xiang & Rix (2022) identify two distinct sequences in the stellar age-metallicity distribution separated at age ~8 Gyr. Similarly, Sahlholdt et al. (2022) propose an age–metallicity relation characterised by several disconnected structures, which could be linked to different star-formation regimes throughout the Milky Way disc evolution.

Figure 8 shows the comparison of the MDFs predicted by our two-infall model and that observed in the APOGEE DR17 data. We can see that, although both distributions have similar median values, their shapes differ. This difference between the predicted and the observed MDFs is more significant in the super metal-rich regime ([Fe/H] ≳ 0.1 dex, Santos-Peral et al. 2021), where our two-infall model sub-estimates the number of sources at that metallicity. This discrepancy can be explained by the effect of the radial migration from the inner Galaxy: stars born in the central high-metallicity regions have experienced a change in their angular momentum due to the interaction with the non-axisymmetric structures of the disc, such as the bar and the spiral arms (Sellwood & Binney 2002; Schönrich & Binney 2009; Minchev et al. 2011). Such migrated population shapes the metal-rich tail of the MDF, increasing its skewness as reported by Hayden et al. (2015) and latterly confirmed by Loebman et al. (2016); Martinez-Medina et al. (2016, 2017). Since Eq. (12) does not include any term associated with the radial migration, our analytic solution predicts a lower number of super-solar metallic-ity stars. However, by blurring the distribution using a Gaussian Kernel Density Estimator of width 0.1 dex, we obtain a smooth distribution whose shape agrees better with that of the APOGEE MDF in the super-solar regime.

We evaluate the dependence of the two-infall chemical evolution model on the DTD by repeating the previous analysis with the WIDE G05 DTD. We discarded the use of the CLOSE G05 DTD for this test because no adequate combination of infall parameters has been found. In short, this DTD provides a large fraction of prompt events, so that the Fe enrichment occurs very fast and the MDF results overpopulated at high metallicities for all the realistic options of infall parameters tested (see Appendix C). Similarly, we use a different set of parameters compared to the MR01 case because no satisfactory common parameters have been found. For the WIDE G05 DTD, the parameters of the infall that better reproduce the APOGEE DR17 data are τ₁ = 0.13 Gyr, τ₂ = 6.75 Gyr, t₁ = 0 Gyr, t₂ = 3.5 Gyr, A₁ ≈ 52.887 M_⊙ pc⁻² Gyr⁻¹ and A₂ ≈ 6.898 M_⊙ pc⁻² Gyr⁻¹ (also equivalent to a mass ratio between infalls of 5.3), while the starformation efficiency and the wind-loading factor have been set to v_l = 0.8 Gyr⁻¹ and ω = 0.2, respectively. As in the MR01 case, we consider a rescaled version of the silicon yields by applying a factor 0.8 to the nominal values presented in Sect. 3.2. As we can see in Fig. 10, the resulting chemical evolution track is able to reproduce the observed [Si/Fe] versus [Fe/H]. In the high-α regime, the solution with the WIDE G05 DTD shows a similar trend to that found with the MR01 DTD, while for the low-α sequence it requires an earlier second infall and a more extended tail loop to trace the chemical evolution. Compared to Fig. 8, the model with the WIDE G05 DTD results in a more metal rich MDF, with a peak at [Fe/H] ≈ 0.09 dex and a larger discrepancy with the median metallicity of the APOGEE DR17 sample (Δ [Fe/H] = 0.08 dex).

Fig. 6 Age-metallicity relation for the one-infall model (see Sect. 3.2) for different DTDs. Colour code is the same as in Fig. 4.

Fig. 7 Fraction of iron produced by Type Ia SNe as a function of age (left vertical axis) assuming the MR01 (blue line), WIDE and CLOSE G05 (purple and pink, respectively), MVP06 (green), S05 (orange), T08 (cyan), and P08 DTDs (brown). The right vertical axis indicates the fraction of iron produced by Type II SNe.

Fig. 8 Predictions of the two-infall model assuming the MR01 DTD. Left panel: comparison of the [Si/Fe] versus [Fe/H] predicted by our two-infall model in the solar neighbourhood (see model details in Sect. 4) with the observational data of APOGEE DR17 (Abdurro’uf et al. 2022) in the Galactic region between 7.2 and 9.2 kpc. The colour coding throughout the model track stands for the ages of the SSPs formed during the Galactic disc evolution. The white circle represents the [Si/Fe] versus [Fe/H] value at the beginning of the second infall (tG – t2 = 10.8 Gyr). Upper-right panel: MDF expected from the two-infall model shown in Fig. 8 (blue histogram) compared to the MDF of the APOGEE DR17 sample (red empty histogram). The vertical lines indicate the median values of each distribution. Middle-right panel: age–metallicity (black line) and [Si/Fe] versus age (blue line) relations predicted by the same model as the upper panel. Lower-right panel: Temporal evolution of the Type Ia SN rate for the whole Galactic disc is indicated with the blue line.

Fig. 9 Same as left panel in Fig. 8, but considering an alternative model with a smaller loading factor for the wind (ω = 0.2), reducing both 〈ySi〉 and 〈mSi,Ia〉 for silicon by a factor of 15%.

5 The Milky Way-like disc in the cosmological context

Numerical simulations constitute an important tool for the study of the formation and evolution of galaxies (Vogelsberger et al. 2020). They allow for comparisons of the structure, kinematics, and chemical composition inferred from the observational data with the models. Some simulations, however, are limited by their lack of chemical information. We can overcome this restriction by modelling the chemistry in this simulations with the analytic solution presented in this work.

In this section, we propose applying our analytical model to one simulated Milky-Way like galaxy, nicknamed GALACTICA, introduced in Park et al. (2021). GALACTICA is extracted from a zoom-in hydrodynamical simulation in a cosmological context (i.e. the region of interest, including the galaxy and its host dark matter halo, has been re-simulated at much higher resolution), using the same spatial resolution (~40 pc) and the same sub-grid models than the NewHorizon simulation (see Dubois et al. 2021). Using the SFH of GALACTICA, we aim to model the chemistry of that Galaxy by applying our analytic solution to Eq. (9). Figure 11 shows the SFR of GALACTICA (same as in Fig. 10 of Park et al. 2021) normalised to a total integrated mass of 1 M_⊙ and re-scaled in time to set its age at t_G = 13.8 Gyr (3% older than in the original work).

To fix some fitting problems, we oversampled the binned SFR using a third order spline interpolator evaluated in a grid of time nodes of step size Δt = 0.005 Gyr (black curve in Fig. 11). We smoothed this curve by performing a convolution with a Gaussian kernel of width 0.15 Gyr. The resulting SFR is multiplied by a renormalisation factor to keep the total mass of 1 M_⊙ fixed (blue curve). As Fig. 11 shows, this process redistributes the stellar mass near the most prominent peaks and leads to a less spiky stellar formation history.

Analogously to Eq. (8), we propose a fit for the smoothed SFR by a superposition of functions ψ_k(t) of the form: $${\psi }_k\left(t\right)=C_k\cdot \theta \left(t-t_k\right)\left[\text{exp\hspace{0.17em}}\left(-\frac{\left(t-t_k\right)}{{\tau }_k}\right)-\text{exp\hspace{0.17em}}\left(-\alpha \left(t-t_k\right)\right)\right],$$(20)

where the free parameters are: amplitudes, C_k, timescales, τ_k, and offsets, t_k. The parameter a is given as input and set to 2.23 Gyr⁻¹ (equivalent to considering v_L = 2 Gyr⁻¹ and ω = 0.4). In contrast to Eq. (19), we did not include the ~ exp(−αt) term because we find better a fit without it. Among the free parameters, only the amplitudes can be determined exactly by the linear least-squares fitting method. On the contrary, the timescales and offsets require more advanced optimisation techniques whose convergence can be very slow. For this reason, we propose the following alternative for τ_k and t_k: (i) we select the position of the local maxima of the interpolated SFR indicated (t_max) in Fig. 11 (vertical dashed lines). (ii) We consider the set of timescales from τ_k = 1 Gyr to τ_k = 5 Gyr (step size of 1 Gyr), adding τ_k = 8, 10 Gyr. (iii) Differentiating both sides of Eq. (20), we find the maximum of ψ_k(t) is located at t_max = t_k + τ_k ln (ατ_k)/(ατ_k – 1). Thus, we solve for t_k for all the possible combinations (τ_k, t_max). If the resulting offset is negative, we substitute it by the half of the minimum positive t_k. (iv) We construct a basis of ψ_k functions with all the combinations of (τ_n, t_k), but excluding the cases in which both τ_n and tmax are larger than 6 Gyr. This results in a set of 60 functions ψ_k(t). (v) Finally, using the standard least-squares fitting method we compute the values of the amplitudes C_k.

The resulting fitting function is illustrated in Fig. 11 (dashed red curve). Although some deviations are observed, it traces the general trend of the smoothed SFR. In order to get realistic chemical evolution tracks, we multiply the SFR by a factor to get the present-day stellar density of σ_*(t_G) = 33.4 ± 3 M_⊙pc⁻² (McKee et al. 2015).

In Fig. 12, we show the [Si/Fe] versus [Fe/H] abundance ratios predicted by our analytic model from the fitting SFR mentioned above. We note that for ages older than 10 Gyr, the SFH traced by the fitting function shows three peaks without significant extended quenching periods between them. Their imprint on the [Si/Fe] versus [Fe/H] is characterised by a “smooth” evolution, with a mild dilution signature associated with the second peak of star formation. Nevertheless, some burst features can be found, as previously discussed in the analysis of the high-α sequence of the Milky Way-like zoom-in cosmological VINTERGATAN simulation by Agertz et al. (2021). At more recent ages, the subsequent infalls of pristine gas produce a depletion of [Fe/H], especially in the age intervals 7−6 Gyr and 3−2 Gyr, while it increases during the following extended periods of low star formation activity. As expected, the “loop” features become more prominent and extended at recent times, creating a low-α structure at super-solar metallicities.

It is worth mentioning that our analytic solution has been computed imposing a constant value for the star formation efficiency, v_L. As discussed in Spitoni et al. (2023), two chemical evolution models constrained by the same SFH can lead to different enrichment of silicon and iron just imposing a less massive gas infall and a higher star formation efficiency. In their Fig. 5, they show that for this case, the dilution is substantially diminished. Thus, with the presented analytic solution, we maximise the dilution effect (through a massive infall of pristine gas) since v_L cannot increase during the Galactic evolution nor mimic the strong star burst phases.

We checked that in the GALACTICA simulation the peaks of the star formation (Fig. 11) are associated with a rapid increase of the gas mass. Therefore, the scenario proposed by our model, in which the SFH is the result of 60 subsequent events of gas infall, is valid. In any case, according to several chemical evolution models (Spitoni et al. 2019b; Palla et al. 2020; Lian et al. 2020) and chemo-dynamical simulations (Agertz et al. 2021; Khoperskov et al. 2021; Vincenzo & Kobayashi 2020), the dilution effect originated by the accretion of pristine (or mildly chemical enriched) gas is expected to dominate the chemical enrichment of the low-α sequence.

Fig. 10 Similar to Fig. 8, but assuming the WIDE G05 DTD. The yields 〈ySi〉 and 〈mSi,Ia〉 for silicon have been re-scaled by a factor of 80%.

Fig. 11 SFR of the GALACTICA simulated galaxy explored in Park et al. (2021) normalised to a total mass of 1 M⊙. The grey histogram represents the reported SFR, while the solid blue curve corresponds to the smoothed version of the interpolated curve (solid black line). The fit of the smoothed curve is represented by the dashed red line. The positions of the maxima, tmax, used in the determination of tk (see the text) are indicated by the vertical dashed lines.

Fig. 12 Chemical information extracted from the star formation history shown in Fig. 11. Upper panel: [Si/Fe] versus [Fe/H] diagram. The evolution track is colour-coded with the stellar ages while the star markers indicate the maxima in the SFR (tmax). The evolutionary track is blurred with a 2D Gaussian kernel of width 0.05 dex (greyscale contour plot). Lower panel: age-metallicity (solid black line) and [Si/Fe] versus age (solid grey line) relations predicted by the same model as the upper panel. Dashed vertical lines denote the ages associated with tmax (star markers in the upper panel and dashed vertical lines in Fig. 11) and follow the same colour code convention as the upper panel.

6 Conclusions and future work

In this work, we present a new analytic solution to the Galactic chemical evolution model which can be used with different prescriptions of the DTD, including the single and double degenerate scenarios. We provide some examples of possible applications of our solution and our main conclusions are summarised as follows:

• We prove that our solution can constitute a useful tool for Galactic archaeology as it is able to interpret the chemical APOGEE DR17 disc stars. The analytic solution can reproduce the expected chemical evolution of the α and iron-peak elements. In particular, we compare the pattern in the [Si/Fe] versus [Fe/H] plane observed by APOGEE DR17 with these predicted by two different models: one assuming a single-degenerate DTD and another that considers the double-degenerate DD scenario. In both cases, we find the low-α sequence can be explained by a delayed gas infall, in agreement with the results of detailed numerical models, but considering different Galactic and infall parameters;

• The super-solar metallicity regime observed in APOGEE DR17 is poorly reproduced by our solution since the considered chemical evolution model does not include radial migration terms. However, the blur of the predicted MDF with a Gaussian kernel of width 0.1 dex improves the comparison;

• According to our tests, it is not possible to discern the best DTD for reproducing the data. With the suitable realistic combination of parameters, both the MR01 and the WIDE G05 DTDs can predict the two-sequence pattern seen in the [Si/Fe] diagram, as well as the approximated shape for the MDF. In order to break this degeneracy, more constraints based on accurate stellar ages, more precise stellar yields and gas infall timing among others are required;

• By modelling the chemistry of a simulated Milky Way-like galaxy from its star formation history, we exploit the applicability of our solution in a cosmological context. The study presented here for the GALACTICA simulation constitutes a preliminary work which will be extended with galaxies of different morphology and formation history.

In future, we plan to include in our solution the contribution of periodic perturbations, such as those caused by the spiral arms (Spitoni et al. 2019a; Poggio et al. 2022; Palicio et al. 2023) and bars (Palicio et al. 2018). We also aim to extend our analytic solution to 2D and 3D models by including gas flows, the transport of metals, as well as radial migration.

Since the analytic solution presented here can be used to model dwarf galaxies, it is possible to perform Bayesian fits of Local Group galaxies with our solution (Johnson et al. 2022). Similarly, the chemistry of galaxies with different morphology can be addressed. For instance, for early-type galaxies it should be possible to compare the predictions for 〈[α/Fe]〉 with the results of the MaNGA survey (Liu 2020). Moreover, it will possible to characterise the star-forming objects which obey to scaling-relations, such as the main sequence star formation (Spitoni et al. 2020, Spitoni et al. 2021a) providing their [α/Fe] evolution over time.

Acknowledgements

We want to thank the helpful comments from P. de Laverny. P. A. Palicio acknowledges the financial support from the Centre national d’études spatiales (CNES). E. Spitoni and A. Recio-Blanco received funding from the European Union’s Horizon 2020 research and innovation program under SPACE-H2020 grant agreement number 101004214 (EXPLORE project).

Appendix A Fitting procedure of the DTDs

As mentioned in Section 2.2, the DTDs whose functional form cannot be exactly described by Eq. 4 require a fitting approximation in order to be evaluated in the analytic solutions. In this section, we detail the procedure performed for the individual modelling of such DTDs and evaluate their associated errors.

Appendix A.1 Case of the MR01 DTD

According to MR01 and Fig. 1, the MR01 DTD is described by a piece-wise function of two components connected at t_knee ≈ 1.612 Gyr whose shape and slope depends on the γ parameter, where we assume γ = 0.5 as in Bonaparte et al. (2013). We model the leftmost component as follows: i) we perform an initial Gaussian Mixture fitting assuming N_G Gaussian distributions (Vanderplas et al. 2012; Ivezić et al. 2014). As a result, we obtain three N_G-dimensional arrays with the mean values (τ'), the widths (σ') and the amplitudes (A_G) of the Gaussian curves. We substitute the lowest value in τ' by 0.094 Gyr because we find a better fit of the MR01 DTD. ii) We propose a set of N_E decreasing exponential curves whose characteristic timescales (τ_D) are multiple of a fundamental timescale, τ_f (i.e. the n-th component of τ_D is n · τ_f, with 1 ≤ n ≤ N_E). The choice for τ_f is motivated by the naïve fitting of the modal and knee points (A and B of Fig. 1) with a single exponential curve, resulting in a timescale of 0.57 Gyr. We divide this value by a factor of two to account for shorter timescales; thus, τ_f ≈ 0.57 Gyr. iii) We add a ~ t⁻¹ term to the fitting function. iv) Fixing the non-linear parameters τ', σ' and τ_D, we search the values of the N_G + N_E + 1 amplitudes (A_G,i, A_E,i and A_I,_i, respectively) that minimise the discrepancy with the MR01 DTD. By using Lagrange multipliers, we impose four additional constraints on the least-squares fitting algorithm: we fix the values of the fit at t = 0.03 Gyr, ~ 0.094 Gyr (the maximum) and 1.612 Gyr (the “knee”), and impose zero derivative at the maximum. v) We repeat this procedure testing different combinations of N_G and N_E to find a good compromise between the complexity of the fitting and the similarity with the MR01 DTD. Based on these tests, we select the combination N_G = 3, N_E = 5 (see the first five rows in Table 1).

For the right-most part of the MR01 DTD (t > t_knee), we repeat the previous procedure with the following modifications: i) no imposed values are used for τ'. ii) The fundamental timescale, τ_f, is computed using the coordinates of the “knee” and the minimun (located at t = 13.8 Gyr), which leads to τ_f ≈ 1.79 Gyr. iii) Using the restricted least-squares fitting algorithm, we fix the values at the edges of the interval [t_knee, 13.8 Gyr]. iv) After testing different combinations, we consider the case with N_E = 4 and no Gaussian distributions the best choice for this part of the MR01 DTD.

The resulting set of parameters for this DTD are summarised in the upper part of Table 1. Figure A.1 illustrates the comparison of the original MR01 DTD with its fit, whose maximum discrepancy in absolute value is 0.018. We note the fitting function is defined on a shorter time interval than that of the MR01 because the former becomes negative when t ≲ 30 Myr, while the MR01 DTD is defined for t ≳ 27 Myr. This interval of ~ 3 Myr, however, has a negligible contribution to the total area (~ 0.02 %) of the MR01 DTD.

Appendix A.2 Case of the WIDE G05 DTD

As Fig. 1 illustrates, the WIDE G05 DTD increases asymptotically up to its maximum to decrease at later times following a power-law like relation, in which its slope is defined by the β_a parameter (see Section 4.3.1 of Greggio 2005). For illustrative purposes, in this work we consider the intermediate case β_a=0. The transition between the mentioned two regimes is determined by the nuclear timescale of the least massive secondary in Type Ia SN progenitor systems τ_n,x, whose value of 0.4 Gyr adopted in this work implies a mass of 3 M_⊙. Using τ_n,x as reference, we define two time intervals to perform the fitting.

Fig. A.1 Delay time distribution of Matteucci & Recchi (2001) normalised to the maximum (dashed black curve) and its fit using three Gaussian curves, nine exponentials, and two t−1 functions (solid blue line). Vertical dotted lines indicate the time domain (τ1,τ2)G,E,I of the fitting for each section of the DTD (see Table 2 for the adopted parameters).

For the t ≤ τ_n,x interval, we perform a restricted least-squares fitting procedure similar to those considered for the MR01 DTD. We model the increasing part of the G05 DTD with a Gaussian curve centred at τ' = 0.θ9 Gyr, an exponential function with a timescale of τ_D = 3.19 Gyr and a 1/t relation imposing continuity at τ_n,x. In contrast, the decreasing regime (t > τ_n,x) requires a more complex fitting function: i) we propose a set of five Gaussian curves with τ' = τ_n,x and widths σ' ranging from 0.1 to 0.5 Gyr (step 0.1 Gyr). ii) For the exponential term, we create a partition of the time interval [τ_n,x, 13.8 Gyr] into six sub-intervals and repeat the procedure described in Section A.1 in each subdivision. The resulting timescales, τ_D, are summarised in Table 1. iii) We introduce an offset τ₀ = 0.35 Gyr in the 1/t relation to improve the fitting in the t → τ_n,x regime, where the slope of the G05 DTD becomes steeper. iv) The amplitudes of the twelve fitting functions described above are optimised by the restricted least-squares method fixing the values at t = τ_n,x and 13.8 Gyr, and the first derivative at t = 13.8 Gyr.

The maximum discrepancy between the normalised G05 DTD and its fit is 0.041 (see right panel in Fig. A.2). Similarly to the case of the MR01 DTD, we find a difference of ~ 3 Myr between the time domains of the fit and the original DTD that excludes a negligible fraction of the WIDE G05 DTD area (~ 0.04%).

Fig. A.2 WIDE (left panel) and CLOSE (right panel) double degenerate delay time distributions of Greggio (2005) normalised to their maximum (dashed black curves). The form of the WIDE G05 DTD is defined by the parameters τn,x = 0.4 Gyr, βa = 0, while for the CLOSE G05 DTD, we use the combination τn,x = 0.4 Gyr, βg = −0.975. The WIDE DTD is modelled as six Gaussian curves, eight exponentials, and two t−1 functions (solid purple line); while the CLOSE DTD requires twelve Gaussian curves, nine exponentials, and four t−1 functions (solid pink line). Vertical dotted lines indicate the time domain (τ1, τ2)G,E,I of the fitting for each section of the DTD (see Table 2 for the adopted parameters).

Appendix A.3 Case of the CLOSE G05 DTD

As in the previous cases, we identify two regimes in the CLOSE double degenerate G05 DTD (pink curve in Fig. 1) connected at t = τ_n,x = 0.398 Gyr, where the slope in the power-law regime is determined by the parameter β_g (Greggio 2005, see Section 4.3.2 of). In this work, we illustrate the case β_g = −0.975 (i.e. a very steep DTD). The leftmost part of this DTD can be modelled using the following fitting functions: i) six Gaussian curves whose parameters are determined by a Gaussian Mixture process, imposing the value of the lowest offset τ' to 0.102 Gyr. ii) Four exponential curves whose timescales are estimated using the same partition procedure as for the WIDE DTD case. iii) One t⁻¹ function. The contribution of these eleven functions are optimised by a restricted least-squares algorithm fixing the values at t = 0.04 Gyr, 0.398 Gyr and the maximum at 0.102 Gyr. Similarly, the rightmost part of the CLOSE G05 DTD (t > 0.102 Gyr) is modelled imposing: i) five exponentials with τ_D determined by the partition of the [0.102 Gyr, 13.8 Gyr] interval. ii) Three (t − τ₀)⁻¹ functions with τ₀ = 0.25, 0.30 and 0.33 Gyr. iii) Six Gaussian curves whose offsets τ' and widths σ' are updated in an iterative process based on the mismatch between the original CLOSE G05 DTD and its fit. The resulting values are summarised in the third and fourth columns of Table 1. where the constraints fix the values at t = 0.102 Gyr, 13.8 Gyr as well as the slope at t = 13.8 Gyr. This fit results in a maximum discrepancy with the original CLOSE G05 DTD of 0.023 (left panel in Fig. A.2).

Appendix A.4 Case of the P08 DTD

Since the P08 DTD has the form ~ t^−1/2 we do not need to use piece-wise functions for its approximation as with the previous DTDs. We model the P08 DTD as a combination of the following functions: i) a Gaussian curve with τ' = 3.5 × 10⁻³ Gyr and σ' = 0. 1 Gyr. ii) Four exponential curves whose timescales, τ_D, are computed as in the G05 case. iii) Four (t − τ₀)⁻¹ functions, where τ₀ ∈ {0.025, 0.02, 0.015, 0.01}. The amplitudes of these nine basis functions are optimised imposing the exact values at the limits of the fitting interval (see Table 1). This results in a discrepancy between the P08 DTD and its fit lower than 0.017 (Fig. A.3).

Fig. A.3 Delay time distribution of Pritchet et al. (2008) normalised to the maximum (dashed black curve) and its fit using one Gaussian curve, four exponentials, and four (t − τ0)−1 functions (solid brown line). Vertical dotted lines indicate the time domain (τ1, τ2)G,Ε,I of the fitting for each section of the DTD (see Table 2 for the adopted parameters).

Appendix B Analytic solution

Appendix B.1 Case τ_j ≠ α⁻¹

In order to simplify the analytic expressions, it is useful to define the following parameters and functions: $${\beta }_j=\alpha -{\tau }_j^{-1},$$(B.1) $$\begin{array}{cc}{\lambda }_x\left(t\right)=\min \left(t,{\tau }_{1x}\right),& \text{with\hspace{0.17em}}x\in \left\{G,E,I\right\},\end{array}$$(B.2) $$\begin{array}{cc}{\text{Λ}}_x\left(t\right)=\min \left(t,{\tau }_{2x}\right),& \text{with\hspace{0.17em}}x\in \left\{G,E,I\right\},\end{array}$$(B.3) $${\eta }_j=\tau \prime +\sigma {\prime }^2/{\tau }_j,$$(B.4) $${\eta }_{\alpha }=\tau \prime +\sigma {\prime }^2\alpha ,$$(B.5) $$Q_D\left(x;a|c\right)=\{\begin{array}{ll}\frac{a\cdot c}{a-c}\cdot \text{exp\hspace{0.17em}}\left(\frac{\left(c-a\right)}{a\cdot c}\cdot x\right)\hfill & \text{\hspace{0.17em}\hspace{0.17em}if\hspace{1em}}a\ne c,\hfill \\ -x\hfill & \text{\hspace{0.17em}\hspace{0.17em}if\hspace{1em}}a=c,\hfill \end{array}$$(B.6) $$P_D\left(x;a,b|c\right)=c\times \{\begin{array}{ll}\frac{a\cdot b\cdot c}{\left(a-c\right)\left(b-c\right)}\cdot \text{exp\hspace{0.17em}}\left(\frac{\left(c-b\right)}{b\cdot c}\cdot x\right)\hfill & \text{if\hspace{1em}}c\ne b,\wedge c\ne b,\hfill \\ -\frac{a}{\left(a-c\right)}\cdot x\hfill & \text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{1em}}c\text{\hspace{0.17em}=\hspace{0.17em}}b,\hfill \\ \frac{b}{{\left(b-c\right)}^2}\cdot \left[\left(c-b\right)\cdot x-c\cdot b\right]\cdot \text{exp\hspace{0.17em}}\left(\frac{\left(c-b\right)}{b\cdot c}\cdot x\right),\hfill & \text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{1em}}c\text{\hspace{0.17em}=\hspace{0.17em}}a,\hfill \end{array}$$(B.7) $$S_D\left(x;a|c\right)=\{\begin{array}{l}\frac{a^2\cdot c^2}{{\left(a-c\right)}^2}\cdot \text{exp\hspace{0.17em}}\left(\frac{\left(c-a\right)}{a\cdot c}\cdot x\right)\text{\hspace{1em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{1em}}c\ne a,\hfill \\ \frac{x^2}{2}\text{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}\hspace{0.17em}if\hspace{1em}}c\text{\hspace{0.17em}=\hspace{0.17em}}a.\hfill \end{array}$$(B.8)

Using these definitions, the Type Ia SN rates are as follows: $$\begin{array}{l}{ℛ}_{Ia}^G\left(t\right)=C_{Ia}{A}_G\sigma \prime v_L\sqrt{\frac{\pi }{2}}\{{\displaystyle \sum _{\begin{array}{c}j=1\\ {\tau }_j\ne {\alpha }^{-1}\end{array}}^N\frac{A_j}{{\beta }_j}\cdot \text{exp\hspace{0.17em}}\left(-\frac{1}{{\tau }_j}\left[\text{Δ}{t}_j-\tau \prime -\frac{\sigma {\prime }^2}{2{\tau }_j}\right]\right)\cdot \left[\text{erf\hspace{0.17em}}\left(\frac{{\text{Λ}}_G\left(\text{Δ}{t}_j\right)-{\eta }_j}{\sqrt{2}\sigma \prime }\right)-\text{erf\hspace{0.17em}}\left(\frac{{\lambda }_G\left(\text{Δ}{t}_j\right)-{\eta }_j}{\sqrt{2}\sigma \prime }\right)\right]}\hfill \\ -{\displaystyle \sum _{\begin{array}{c}j=1\\ {\tau }_j\ne {\alpha }^{-1}\end{array}}^N\frac{A_j}{{\beta }_j}\cdot \text{exp\hspace{0.17em}}\left(-\alpha \left[\text{Δ}{t}_j-\tau \prime -\frac{\sigma {\prime }^2\alpha }{2}\right]\right)\cdot \left[\text{erf\hspace{0.17em}}\left(\frac{{\text{Λ}}_G\left(\text{Δ}{t}_j\right)-{\eta }_{\alpha }}{\sqrt{2}\sigma \prime }\right)-\text{erf\hspace{0.17em}}\left(\frac{{\lambda }_G\left(\text{Δ}{t}_j\right)-{\eta }_{\alpha }}{\sqrt{2}\sigma \prime }\right)\right]}\hfill \\ +{\sigma }_{gas}\left(0\right)\cdot \text{exp\hspace{0.17em}}\left(-\alpha \left[t-\tau \prime -\frac{\sigma {\prime }^2\alpha }{2}\right]\right)\cdot \left[\text{erf\hspace{0.17em}}\left(\frac{{\text{Λ}}_G\left(t\right)-{\eta }_{\alpha }}{\sqrt{2}\sigma \prime }\right)-\text{erf\hspace{0.17em}}\left(\frac{{\lambda }_G\left(t\right)-{\eta }_{\alpha }}{\sqrt{2}\sigma \prime }\right)\right]\},\hfill \end{array}$$(B.9a) $$\begin{array}{l}{ℛ}_{Ia}^E\left(t\right)=C_{Ia}{A}_E{v}_L\{{\displaystyle \sum _{\begin{array}{c}j=1\\ {\tau }_j\ne {\alpha }^{-1}\end{array}}^N\frac{A_j}{{\beta }_j}\text{exp\hspace{0.17em}}\left(-\frac{\text{Δ}{t}_j}{{\tau }_j}\right)\left[Q_D\left({\lambda }_E\left(\text{Δ}{t}_j\right);{\tau }_j|{\tau }_D\right)-Q_D\left({\text{Λ}}_E\left(\text{Δ}{t}_j\right);{\tau }_j|{\tau }_D\right)\right]}\hfill \\ -{\displaystyle \sum _{\begin{array}{c}j=1\\ {\tau }_j\ne {\alpha }^{-1}\end{array}}^N\frac{A_j}{{\beta }_j}\cdot \text{exp\hspace{0.17em}}\left(-\alpha \text{Δ}{t}_j\right)\left[Q_D\left({\lambda }_E\left(\text{Δ}{t}_j\right);{\alpha }^{-1}|{\tau }_D\right)-Q_D\left({\text{Λ}}_E\left(\text{Δ}{t}_j\right);{\alpha }^{-1}|{\tau }_D\right)\right]}\hfill \\ +{\sigma }_{gas}\left(0\right)\cdot \text{exp\hspace{0.17em}}\left(-\alpha t\right)\left[Q_D\left({\lambda }_E\left(t\right);{\alpha }^{-1}|{\tau }_D\right)-Q_D\left({\text{Λ}}_E\left(t\right);{\alpha }^{-1}|{\tau }_D\right)\right]\},\hfill \end{array}$$(B.9b) $$\begin{array}{l}{ℛ}_{Ia}^I\left(t\right)=C_{Ia}{A}_{\text{I}}{\tau }_I{v}_L{\displaystyle \sum _{\begin{array}{c}j=1\\ {\tau }_j\ne {\alpha }^{-1}\end{array}}^N\frac{A_j}{{\beta }_j}\{\text{exp\hspace{0.17em}}\left(-\frac{\text{Δ}{t}_j-{\tau }_0}{{\tau }_j}\right)\cdot \left[\text{Ei\hspace{0.17em}}\left(\frac{{\text{Λ}}_I\left(\text{Δ}{t}_j\right)-{\tau }_0}{{\tau }_j}\right)-\text{Ei}\left(\frac{{\lambda }_I\left(\text{Δ}{t}_j\right)-{\tau }_0}{{\tau }_j}\right)\right]+}\hfill \\ +\text{\hspace{0.17em}exp\hspace{0.17em}}\left(-\alpha \left[\text{Δ}{t}_j-{\tau }_0\right]\right)\cdot \left[\text{Ei}\left(\alpha \left[{\lambda }_I\left(\text{Δ}{t}_j\right)-{\tau }_0\right]\right)-\text{Ei\hspace{0.17em}}\left(\alpha \left[{\text{Λ}}_I\left(\text{Δ}{t}_j\right)-{\tau }_0\right]\right)\right]\}\hfill \\ +C_{Ia}{A}_{\text{I}}{\tau }_I{v}_L{\sigma }_{gas}\left(0\right)\text{\hspace{0.17em}exp\hspace{0.17em}}\left(-\alpha \left[t-{\tau }_0\right]\right)\cdot \left[\text{Ei\hspace{0.17em}}\left(\alpha \left[{\text{Λ}}_I\left(t\right)-{\tau }_0\right]\right)-\text{Ei\hspace{0.17em}}\left(\alpha \left[{\lambda }_I\left(t\right)-{\tau }_0\right]\right)\right],\hfill \end{array}$$(B.9c)

so that the global Type Ia SN rate is ${ℛ}_{Ia}={ℛ}_{Ia}^G+{ℛ}_{Ia}^E+{ℛ}_{Ia}^I$, where Δt_j = t − t_j. The Ei function presented in B.9c refers to the so-called exponential integral, though for computational purposes it is better to use a modified version $\tilde{\text{E}}\text{i}$ without divergences at the origin: $$\tilde{\text{E}}\text{i}\left(x\right)=\{\begin{array}{l}Ei\left(x\right)\equiv {\displaystyle {\int }_{-\infty }^x{y}^{-1}{e}^{y}dy,\text{\hspace{1em}\hspace{0.17em}\hspace{0.17em}if\hspace{1em}}x\ne 0,}\hfill \\ 0,\text{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{1em}}x=0.\hfill \end{array}$$(B.10)

We separate the individual contribution of each term in Eq. 10 to the global solution σ_X as σ_X = σ_X,IRA + σ_X,Ia. The term related to the IRA approximation is expressed as: $$\begin{array}{l}{\sigma }_{X,IRA}\left(t\right)={\sigma }_X\left(0\right)\text{exp\hspace{0.17em}}\left(-\alpha t\right)+\langle y_X\rangle \left(1-R\right)v_L{\displaystyle \sum _{\begin{array}{c}j=1\\ {\tau }_j\ne {\alpha }^{-1}\end{array}}^N\frac{A_j}{{\beta }_j}\text{exp\hspace{0.17em}}\left(-\alpha \text{Δ}{t}_j\right)\theta \left(\text{Δ}{t}_j\right)\left[\frac{\text{exp}\left({\beta }_j\text{Δ}{t}_j-1\right)}{{\beta }_j}-\text{Δ}{t}_j\right]}+\hfill \\ +\langle y_X\rangle \left(1-R\right)v_L{\sigma }_{gas}\left(0\right)\cdot t\cdot \text{exp\hspace{0.17em}}\left(-\alpha t\right)\theta \left(t\right),\hfill \end{array}$$(B.11)

We can write σ_X,Ia as the sum of the three components of the DTD as σ_X,Ia = σ_X,Ia,G + σ_X,Ia,E + σ_X,Ia,I. For sake of illustration, we consider the case in which σ_X,Ia,G, σ_X,Ia,E and σ_X,Ia,I are defined by a unique Gaussian, exponential, and (t − τ₀)⁻¹ distributions, respectively. This allows us to obviate the i-index in the parameters of the DTD (see Eq. 4) and simplify the notation. For more realistic DTDs, like the ones considered in this work, it is necessary to compute the associated σ_X,Ia,G, σ_X,Ia,E or σ_X,Ia,I of each individual component of the DTD. Thus, the terms σ_X,Ia,G, σ_X,Ia,E, and σ_X,Ia,I are expressed as: $$\begin{array}{l}{\sigma }_{X,Ia,G}\left(t\right)=\hfill \\ =\langle m_{X,Ia}\rangle C_{Ia}{A}_G{\sigma }^{\prime }{v}_L\sqrt{\frac{\pi }{2}}{\displaystyle \sum _{\underset{{\tau }_j\ne {\alpha }^{-1}}{j=1}}^N\frac{A_j}{{\beta }_j^2}\theta \left(\text{Δ}{t}_j-{\tau }_1\right)\exp \left(-\alpha \text{Δ}{t}_j+\frac{{\tau }^{\prime }}{{\tau }_j}+\frac{{{\sigma }^{\prime }}^2}{2{\tau }_j^2}\right)\cdot \{\left[\text{erf}\left(\frac{{\tau }_1-{\eta }_{\alpha }}{\sqrt{2}{\sigma }^{\prime }}\right)-\text{erf}\left(\frac{{\text{Λ}}_G\left(\text{Δ}{t}_j\right)-{\eta }_{\alpha }}{\sqrt{2}{\sigma }^{\prime }}\right)\right].}\hfill \\ \cdot \exp \left({\beta }_j{\eta }_j+\frac{1}{2}{\beta }_j^2{{\sigma }^{\prime }}^2\right)+\left[\text{erf}\left(\frac{{\text{Λ}}_G\left(\text{Δ}{t}_j\right)-{\eta }_j}{\sqrt{2}{\sigma }^{\prime }}\right)-\text{erf}\left(\frac{{\tau }_1-{\eta }_j}{\sqrt{2}{\sigma }^{\prime }}\right)\right]\exp \left({\beta }_j{\text{Λ}}_G\left(\text{Δ}{t}_j\right)\right)\hfill \\ +\theta \left(\text{Δ}{t}_j-{\tau }_2\right)\left[\exp \left({\beta }_j\text{Δ}{t}_j\right)-\exp \left({\beta }_j{\tau }_2\right)\right]\left[\text{erf}\left(\frac{{\tau }_2-{\eta }_j}{\sqrt{2}{\sigma }^{\prime }}\right)-\text{erf}\left(\frac{{\tau }_1-{\eta }_j}{\sqrt{2}{\sigma }^{\prime }}\right)\right]\}\hfill \\ -\langle m_{X,Ia}\rangle C_{Ia}{A}_G{\sigma }^{\prime }{v}_L{\displaystyle \sum _{\underset{{\tau }_j\ne {\alpha }^{-1}}{j=1}}^N\frac{A_j}{{\beta }_j}\theta \left(\text{Δ}{t}_j-{\tau }_1\right)\exp \left(-\alpha \text{Δ}{t}_j+{\tau }^{\prime }\alpha +\frac{{{\sigma }^{\prime }}^2{\alpha }^2}{2}\right)\cdot \{\sqrt{\frac{\pi }{2}}\cdot [\left({\text{Λ}}_G\left(\text{Δ}{t}_j\right)-{\eta }_{\alpha }\right)\cdot \text{erf}\left(\frac{{\text{Λ}}_G\left(\text{Δ}{t}_j\right)-{\eta }_{\alpha }}{\sqrt{2}{\sigma }^{\prime }}\right)}\hfill \\ -\left({\tau }_1-{\eta }_{\alpha }\right)\cdot \text{erf}\left(\frac{{\tau }_1-{\eta }_{\alpha }}{\sqrt{2}{\sigma }^{\prime }}\right)+\left({\tau }_1-{\text{Λ}}_G\left(\text{Δ}{t}_j\right)\right)\cdot \text{erf}\left(\frac{{\tau }_1-{\eta }_{\alpha }}{\sqrt{2}{\sigma }^{\prime }}\right)]+\hfill \\ {\sigma }^{\prime }\left[\exp \left(-\frac{{\left(\text{Λ}\left(\text{Δ}{t}_j\right)-{\eta }_{\alpha }\right)}^2}{2{{\sigma }^{\prime }}^2}\right)-\exp \left(-\frac{{\left({\tau }_1-{\eta }_{\alpha }\right)}^2}{2{{\sigma }^{\prime }}^2}\right)\right]+\sqrt{\frac{\pi }{2}}\theta \left(\text{Δ}{t}_j-{\tau }_2\right)\cdot \left(\text{Δ}{t}_j-{\tau }_2\right)\left[\text{erf}\left(\frac{{\tau }_2-{\eta }_{\alpha }}{\sqrt{2}{\sigma }^{\prime }}\right)-\text{erf}\left(\frac{{\tau }_1-{\eta }_{\alpha }}{\sqrt{2}{\sigma }^{\prime }}\right)\right]\}+\hfill \\ +\langle m_{X,Ia}\rangle C_{Ia}{A}_G{\sigma }^{\prime }{v}_L{\sigma }_{gas}\left(0\right)\cdot \theta \left(t-{\tau }_1\right)\exp \left(-\alpha t+{\tau }^{\prime }\alpha +\frac{{{\sigma }^{\prime }}^2{\alpha }^2}{2}\right)\cdot \{\sqrt{\frac{\pi }{2}}\cdot [\left({\text{Λ}}_G\left(t\right)-{\eta }_{\alpha }\right)\cdot \text{erf}\left(\frac{{\text{Λ}}_G\left(t\right)-{\eta }_{\alpha }}{\sqrt{2}{\sigma }^{\prime }}\right)\hfill \\ -\left({\tau }_1-{\eta }_{\alpha }\right)\cdot \text{erf}\left(\frac{{\tau }_1-{\eta }_{\alpha }}{\sqrt{2}{\sigma }^{\prime }}\right)+\left({\tau }_1-{\text{Λ}}_G\left(t\right)\right)\cdot \text{erf}\left(\frac{{\tau }_1-{\eta }_{\alpha }}{\sqrt{2}{{\sigma }^{\prime }}^2}\right)]+{\sigma }^{\prime }\left[\exp \left(-\frac{{\left({\text{Λ}}_G\left(t\right)-{\eta }_{\alpha }\right)}^2}{2{{\sigma }^{\prime }}^2}\right)-\exp \left(-\frac{{\left({\tau }_1-{\eta }_{\alpha }\right)}^2}{2{{\sigma }^{\prime }}^2}\right)\right]+\hfill \\ \theta \left(t-{\tau }_2\right)\cdot \left(t-{\tau }_2\right)\sqrt{\frac{\pi }{2}}\cdot \left[\text{erf}\left(\frac{{\tau }_2-{\eta }_{\alpha }}{\sqrt{2}{\sigma }^{\prime }}\right)-\text{erf}\left(\frac{{\tau }_1-{\eta }_{\alpha }}{\sqrt{2}{\sigma }^{\prime }}\right)\right]\}.\hfill \end{array}$$(B.12a) $$\begin{array}{l}{\sigma }_{X,Ia,E}\left(t\right)=\hfill \\ =\langle m_{X,Ia}\rangle C_{Ia}{A}_E{v}_L{\displaystyle \sum _{\underset{{\tau }_j\ne {\alpha }^{-1}}{j=1}}^N\frac{A_j}{{\beta }_j}\exp \left(-\alpha \text{Δ}{t}_j\right)\{Q_D\left({\tau }_2;{\tau }_j|{\tau }_D\right)\left[Q_D\left(\text{Δ}{t}_j;{\alpha }^{-1}|{\tau }_j\right)-Q_D\left({\text{Λ}}_E\left(\text{Δ}{t}_j\right);{\alpha }^{-1}|{\tau }_j\right)\right]}\hfill \\ +Q_D\left({\tau }_1;{\tau }_j|{\tau }_D\right)\left[Q_D\left({\lambda }_E\left(\text{Δ}{t}_j\right);{\alpha }^{-1}|{\tau }_j\right)-Q_D\left(\text{Δ}{t}_j;{\alpha }^{-1}|{\tau }_j\right)\right]+Q_D\left({\tau }_2;{\alpha }^{-1}|{\tau }_D\right)\left[Q_D\left({\text{Λ}}_E\left(\text{Δ}{t}_j\right);{\alpha }^{-1}|{\alpha }^{-1}\right)-Q_D\left(\text{Δ}{t}_j;{\alpha }^{-1}|{\alpha }^{-1}\right)\right]\hfill \\ +Q_D\left({\tau }_1;{\alpha }^{-1}|{\tau }_D\right)\left[Q_D\left(\text{Δ}{t}_j;{\alpha }^{-1}|{\alpha }^{-1}\right)-Q_D\left({\lambda }_E\left(\text{Δ}{t}_j\right);{\alpha }^{-1}|{\alpha }^{-1}\right)\right]\hfill \\ +P_D\left({\text{Λ}}_E\left(\text{Δ}{t}_j\right);{\tau }_j,{\alpha }^{-1}|{\tau }_D\right)-S_D\left({\text{Λ}}_E\left(\text{Δ}{t}_j\right);{\alpha }^{-1}|{\tau }_D\right)-P_D\left({\lambda }_E\left(\text{Δ}{t}_j\right);{\tau }_j,{\alpha }^{-1}|{\tau }_D\right)+S_D\left({\lambda }_E\left(\text{Δ}{t}_j\right);{\alpha }^{-1}|{\tau }_D\right)\}\hfill \\ +\langle m_{X,Ia}\rangle C_{Ia}{A}_E{v}_L{\sigma }_{gas}\left(0\right)\cdot \exp \left(-\alpha t\right)\{Q_D\left({\tau }_2;{\alpha }^{-1}|{\tau }_D\right)\left[Q_D\left(t;{\alpha }^{-1}|{\alpha }^{-1}\right)-Q_D\left({\text{Λ}}_E\left(t\right);{\alpha }^{-1}|{\alpha }^{-1}\right)\right]\hfill \\ -Q_D\left({\tau }_1;{\alpha }^{-1}|{\tau }_D\right)\left[Q_D\left(t;{\alpha }^{-1}|{\alpha }^{-1}\right)-Q_D\left({\lambda }_E\left(t\right);{\alpha }^{-1}|{\alpha }^{-1}\right)\right]+S_D\left({\text{Λ}}_E\left(t\right);{\alpha }^{-1}|{\tau }_D\right)-S_D\left({\lambda }_E\left(t\right);{\alpha }^{-1}|{\tau }_D\right)\}.\hfill \end{array}$$(B.12b) $$\begin{array}{l}{\sigma }_{X,Ia,I}\left(t\right)=\hfill \\ =\langle m_{X,Ia}\rangle C_{Ia}{A}_I{\tau }_I{v}_L{\displaystyle \sum _{\begin{array}{c}j=1\\ {\tau }_j\ne {\alpha }^{-1}\end{array}}^N\frac{A_j}{{\beta }_j^2}\exp \left(-\alpha \left[\text{Δ}{t}_j-{\tau }_0\right]\{\theta \left(\text{Δ}{t}_j-{\tau }_0\right)\cdot \left[\text{Ei}\left(\frac{{\text{Λ}}_I\left(\text{Δ}{t}_j\right)-{\tau }_0}{{\tau }_j}\right)\cdot \exp \left({\beta }_j\left[{\text{Λ}}_I\left(\text{Δ}{t}_j\right)-{\tau }_0\right]\right)\right]\right)+}\hfill \\ +\text{Ei}\left(\alpha \left[{\tau }_1-{\tau }_0\right]\right)-\text{Ei}\left(\alpha \left[{\text{Λ}}_I\left(\text{Δ}{t}_j\right)-{\tau }_0\right]\right)-\text{Ei}\left(\frac{{\tau }_1-{\tau }_0}{{\tau }_j}\right)\cdot \exp \left({\beta }_j\left[\text{Δ}{t}_j-{\tau }_0\right]\right)]+\hfill \\ +\theta \left(\text{Δ}{t}_j-{\tau }_2\right)\cdot \left[\exp \left({\beta }_j\left[\text{Δ}{t}_j-{\tau }_0\right]\right)-\exp \left({\beta }_j\left[{\tau }_2-{\tau }_0\right]\right)\right]\cdot \text{Ei}\left(\frac{{\tau }_2-{\tau }_0}{{\tau }_j}\right)\}\hfill \\ -\langle m_{X,Ia}\rangle C_{Ia}{A}_I{\tau }_I{v}_L{\displaystyle \sum _{\begin{array}{c}j=1\\ {\tau }_j\ne {\alpha }^{-1}\end{array}}^N\frac{A_j}{{\beta }_j}\exp \left(-\alpha \left[\text{Δ}{t}_j-{\tau }_0\right]\right)\{\theta \left(\text{Δ}{t}_j-{\tau }_1\right)\left[{\text{Λ}}_I\left(\text{Δ}{t}_j\right)-{\tau }_0\right]\cdot \text{Ei}\left(\alpha \left[{\text{Λ}}_I\left(\text{Δ}{t}_j\right)-{\tau }_0\right]\right)+}\hfill \\ +\frac{\exp \left(\alpha \left[{\tau }_1-{\tau }_0\right]\right)-\exp \left(\alpha \left[{\text{Λ}}_I\left(\text{Δ}{t}_j\right)-{\tau }_0\right]\right)}{\alpha }-\left(\text{Δ}{t}_j-{\tau }_0\right)\cdot \text{Ei}\left(\alpha \left[{\tau }_1-{\tau }_0\right]\right)]+\theta \left(\text{Δ}{t}_j-{\tau }_2\right)\cdot \text{Ei}\left(\alpha \left[{\tau }_2-{\tau }_0\right]\right)\cdot \left(\text{Δ}{t}_j-{\tau }_2\right)\}+\hfill \\ +\langle m_{X,Ia}\rangle C_{Ia}{A}_I{\tau }_I{v}_L{\sigma }_{gas}\left(0\right)\exp \left(-\alpha \left[t-{\tau }_0\right]\right)\{\theta \left(t-{\tau }_1\right)[\left({\text{Λ}}_I\left(t\right)-{\tau }_0\right)\cdot \text{Ei}\left(\alpha \left[{\text{Λ}}_I\left(t\right)-{\tau }_0\right]\right)-\frac{\exp \left(\alpha \left[{\text{Λ}}_I\left(t\right)-{\tau }_0\right]\right)}{\alpha }+\hfill \\ +\frac{\exp \left(\alpha \left[{\tau }_1-{\tau }_0\right]\right)}{\alpha }-\left(t-{\tau }_0\right)\cdot \text{Ei}\left(\alpha \left[{\tau }_1-{\tau }_0\right]\right)]+\theta \left(t-{\tau }_2\right)\cdot \text{Ei}\left(\alpha \left[{\tau }_2-{\tau }_0\right]\right)\cdot \left(t-{\tau }_2\right)\}.\hfill \end{array}$$(B.12c)