© ESO 2021

1. Introduction

The vast number of galaxies are rotationally-supported disk galaxies with quite constant star-formation histories (Schombert et al. 2019; Kroupa et al. 2020a; Hoffmann et al. 2020), and they comprise the standard outcome of galaxy formation. Elliptical (E) galaxies are rare pressure-supported stellar systems, comprising only a few per cent of the galaxy population (Delgado-Serrano et al. 2010), but they include the most massive galaxies with a dynamical mass (stars and stellar remnants) larger than M_dyn ≈ 10¹⁰ M_⊙. The formation of largely pressure-supported galaxies thus constitutes an interesting problem to solve in a model of cosmological structure formation, particularly since these galaxies are today known to have formed very soon after the Big Bang. E galaxies, furthermore, have a number of properties that are not yet well understood, including their chemical enrichment history.

Dynamical analysis (Poci et al. 2021), spectroscopic studies (Saracco et al. 2020), and the analysis of resolved galaxies (Barbosa et al. 2021; La Barbera et al. 2021) have confirmed that the central regions of massive early-type galaxies (ETGs) form in a monolithic collapse with a short star formation timescale (SFT), and therefore can be reasonably approximated by the ‘closed-box’ model (Matteucci 2012) with neither galactic inflow nor outflow (we return to this in Sect. 6.3). The central and satellite ETGs show no significant differences in their age and α-element gradients, suggesting a similar formation history (Santucci et al. 2020). With a well-constrained galaxy-wide stellar initial mass function (gwIMF), stellar yields, type Ia supernovae (SNIa) normalisation, and delay-time-distribution (DTD), the stellar-element abundance ratio of a galaxy with its present-day mass is determined in the closed-box model solely by the initial gas supply and the star formation history (SFH). One can then estimate the SFT of a monolithically collapsed E galaxy through its stellar mean α-element-to-iron peak element abundance ratio (e.g., [Mg/Fe]). This was achieved by Thomas et al. (2005, 2010), who showed that (i) more massive E galaxies started to form sooner after the Big Bang; (ii) the resulting SFTs constrained by [Mg/Fe], τ_SF, Mg/Fe(M_dyn) become shorter with increasing M_dyn (comprising the downsizing problem as it is contrary to the expectation that more massive galaxies need longer to form through mergers in the standard dark-matter-based cosmological models); and (iii) the τ_SF, Mg/Fe(M_dyn) values are shorter than the τ_SF, SPS(M_dyn) values obtained by stellar population synthesis (SPS) studies by, for example, McDermid et al. (2015) and Lacerna et al. (2020). On the other hand, dark-matter cosmology-motivated hydrodynamical simulations of the hierarchical formation of the most massive ellipticals lead to the result that the empirical τ_SF, Mg/Fe values are too short (less than a Gyr); there not being enough time for the synthesised and released elements to recycle and increase the mean stellar metallicity to the observed level (Colavitti et al. 2009; Pipino et al. 2009; De Lucia et al. 2017; Okamoto et al. 2017; Jafariyazani et al. 2020).

The situation worsens when one considers the galactic stellar metallicity in addition to [Mg/Fe] in a chemical evolution model due to the metal-rich stars having a lower [Mg/Fe] yield (see e.g., Matteucci 2012, Sect. 2.1.3). This was studied in Pipino & Matteucci (2004), where three different masses of galaxies were modelled (shown in Fig. 4 below), although the fits of observed E galaxies are not ideal (cf. Pipino & Matteucci 2011). With the standard assumption of an invariant and canonical gwIMF (Salpeter 1955; Kroupa 2001) and using closed-box modelling, Yan et al. (2019a) obtained a good fit to the metallicity and α element abundances, but the implied τ_SF–M_dyn relation is steeper than the relation suggested originally by Thomas et al. (2005) regardless of a potential systematic bias on the observed metal abundances or on the stellar magnesium yield uncertainty. This revised τ_SF–M_dyn relation suggests an even shorter τ_SF for massive galaxies than that implied by the SPS studies, and it poses an even more severe downsizing problem that cannot be resolved in the standard cosmological hydrodynamical simulations. Different gwIMF and/or stellar yields, SNIa delay-time distributions (DTDs), gas mixing, and expulsion physics, or a combination of these, are probably needed to solve this problem.

The most promising solution is to apply a systematically varying gwIMF. With an increasing number of observational studies supporting this idea, it has become a crucial aspect in a galaxy chemical evolution model and needs to be investigated, given the above problems between theory and observation. Other possible solutions assuming an invariant canonical stellar initial mass function (IMF) are discussed in Sect. 6.4. The gwIMF appears to be top-heavy (containing more massive stars) when the star formation rate (SFR) is high (Gunawardhana et al. 2011; Zhang et al. 2018) and top-light for low-SFR dwarf galaxies (Jeřábková et al. 2018). The gwIMF of low-mass stars is also variable. It becomes bottom-light (containing fewer low-mass stars) in metal-poor environments (Gennaro et al. 2018; Yan et al. 2020) and bottom-heavy for the super-solar-metallicity E galaxies (e.g., Salvador-Rusiñol et al. 2021). This means the gwIMF would change over time depending on the SFH and metal enrichment history, which helps to build up the high-metallicity and [Mg/Fe] ratio for the massive E galaxies. The need for a time-variable gwIMF to explain the galactic abundance has been argued in, for example, Vazdekis et al. (1997), Larson (1998), Weidner et al. (2013a), Narayanan & Davé (2013), Bekki (2013), Ferreras et al. (2015), and Martín-Navarro et al. (2018). We invite the reader to also consult the reviews given by Hopkins (2018) and Smith (2020) giving other reasons for considering IMF variations.

The applications of environment-dependent gwIMFs do show promising results (Recchi et al. 2009; Fontanot et al. 2017; Barber et al. 2018; De Masi et al. 2019; Gutcke & Springel 2019; Palla et al. 2020; Yan et al. 2020), but it is important to use formulations of the gwIMF that are consistent with the extragalactic data and at the same time consistent with observed resolved stellar populations for both massive and low-mass stars (detailed in Sect. 2.1). With this in mind, for the present work we focused on the advanced variable gwIMF theory with various verified predictions, i.e., we applied the integrated galaxy-wide IMF (IGIMF) theory (Kroupa & Weidner 2003; Weidner et al. 2011; Yan et al. 2017, 2020; Jeřábková et al. 2018) to compute the gwIMF. The IGIMF theory is a framework that accounts for the gwIMF being made up of the IMFs of all embedded clusters forming in the galaxy (Kroupa & Weidner 2003) and applies the empirical rules that modify the stellar IMF in embedded clusters according to the metallicity and density of an embedded cluster (Marks et al. 2012; Jeřábková et al. 2018); this is supported by the recent studies of Villaume et al. (2017) and Martín-Navarro et al. (2019a). Thus, here we apply, for the first time, the IGIMF theory to the monolithic E galaxy chemical-evolution model. We study how it affects the SFTs of E galaxies, whether the result is consistent with SPS studies, and what the implications might be. The open-source chemical evolution model developed by Yan et al. (2019b) is applied, accounting specifically for the strong gwIMF variation and calculating the element enrichment self-consistently with the number of SNIa affected by the IMF. The long-term aim of this research project is to ultimately understand how hydrodynamical models of galaxy formation could be made consistent with the observational constraints (mass, α-element abundances, and metallicities of gas and stars) of E galaxies. The closed-box models studied here are to be seen as a comprehensive parameter study, which will allow future computationally expensive fully self-consistent galaxy formation computations to be performed using this knowledge gain.

Section 2 summarises the method used here to calculate the gwIMF and the chemical evolution of a galaxy. In particular, we emphasise how the number of SNIa is affected by the IMF shape. The observed chemical abundances of elliptical galaxies are briefly introduced in Sect. 3. Then, the likelihood of certain SFTs for a galaxy can be determined by comparing the observed and modelled element abundances of galaxies as is detailed in Sect. 4. The resulting most-likely τ_SF–M_dyn relations adopting different IMF and SNIa formulations are shown in Sect. 5. The reliability of the assumptions applied in our model is discussed in Sect. 6. Finally, Sect. 7 contains our conclusions.

2. Galaxy model

This section introduces the GalIMF code, which combines the IGIMF theory (Sect. 2.1) with the galaxy chemical-evolution model (Sect. 2.3). Following Yan et al. (2019b), we introduce the mathematical formulation that calculates the number of SNIa in Sect. 2.2 and demonstrate how the number of SNIa is affected by the IMF variation.

2.1. The systematically varying gwIMF

The IGIMF theory describes how the gwIMF should vary as a function of galactic properties and is computed by assuming that in each star formation epoch of duration δt in which the average SFR¹ is ${\overline{\psi }}_{\delta t}$, the galaxy forms a total mass in stars, $M_{\mathrm{tot}}={\overline{\psi }}_{\delta t}\delta t$, and that this mass is distributed over a fully populated mass function of embedded star clusters, the ECMF². The timescale of δt = 10 Myr, discussed in Schulz et al. (2015), is essentially the lifetime of molecular clouds and is also the timescale during which the interstellar medium of a galaxy churns out a full population of freshly formed embedded clusters, each of which dissolves into the galactic field through gas expulsion; stellar evolution mass loss; and two-body, relaxation-driven evaporation.

The IGIMF theory is based on a set of axioms that are formulated based on observational constraints (e.g., Recchi & Kroupa 2015). These include how the stellar IMF and the ECMF change (Kroupa & Weidner 2003; Weidner et al. 2011; Jeřábková et al. 2018; Yan et al. 2020). The IGIMF theory has solved a number of previously outstanding extragalactic problems, such as explaining the UV extended galactic discs (Pflamm-Altenburg & Kroupa 2008); predicting the diverging Hα- versus UV-fluxes of dwarf galaxies (Pflamm-Altenburg et al. 2009), verified by Lee et al. (2009), a lower α element-to-iron peak element ratio in dwarf galaxies (Yan et al. 2020; Theler et al. 2020; Minelli et al. 2021); and naturally accounting for the timescale problem for building up a sufficient stellar population in dwarf galaxies, given their low SFRs (Pflamm-Altenburg & Kroupa 2009). In the following, we specify the exact formula for calculating the gwIMF that is applied for this work.

The IMF for stars with a mass higher than 1 M_⊙ (Eq. (5), below) follows the prescription given in Yan et al. (2017), while the IMF of low-mass stars (i.e. α₁ and α₂ in Eq. (1)) follows Eq. (9) of Yan et al. (2020). For the stellar IMF (m is in unit of solar mass), Assuming that the stellar mass upper limit, m_max, is higher than 1 M_⊙, we have

$$\begin{array}{c}\hfill {\xi }_{\star }(m)=\mathrm{d}{N}_{\star }/\mathrm{d}m=\{\begin{array}{c}2k_{\star }{m}^{-{\alpha }_1},0.08\le m/M_{\odot }<0.50,\hfill \\ k_{\star }{m}^{-{\alpha }_2},0.50\le m/M_{\odot }<1.00,\hfill \\ k_{\star }{m}^{-{\alpha }_3},1.00\le m/M_{\odot }<m_{\max },\hfill \end{array}\end{array}$$(1)

where the number of stars in the mass interval m to m + dm is dN_⋆. We note that the normalisation parameter for stars in the smallest mass interval, 2k_⋆, is two times the normalisation parameter for more massive stars, because we presume that IMF is a continuous function and that α₂ − α₁ = 1 always holds (see Eq. (4), below). The embedded-cluster-mass-dependent stellar mass upper limit, m_max(M_ecl)≤m_max* = 150 M_⊙, and the normalisation parameter, k_⋆, are determined by simultaneously solving for the mass in stars formed in the embedded cluster,

$$\begin{array}{c}\hfill M_{\mathrm{ecl}}={\displaystyle {\int }_{0.08{\mathrm{M}}_{\odot }}^{m_{\max }}}m{\xi }_{\star }(m)\mathrm{d}m,\end{array}$$(2)

and the mass limit for the stars, m_max, is given by the optimal sampling law (Kroupa et al. 2013; Schulz et al. 2015; Yan et al. 2017),

$$\begin{array}{c}\hfill 1={\displaystyle {\int }_{m_{\max }}^{m_{\max \ast }}}{k}_{\star }{m}^{-{\alpha }_3}\mathrm{d}m;\\ \hfill \end{array}$$(3)

m_max* is the adopted physical upper mass limit, 150 M_⊙ (Yan et al. 2017), given m_max >  1.

The power-law indices or slopes of the IMF for low- and intermediate-mass stars is determined empirically by Yan et al. (2020) as

$$\begin{array}{cc}& {\alpha }_1=1.3+\mathrm{\Delta }\alpha ·(Z-Z_{\odot }),\hfill \\ \hfill & {\alpha }_2=2.3+\mathrm{\Delta }\alpha ·(Z-Z_{\odot }),\hfill \end{array}$$(4)

where Z and Z_⊙ = 0.0142 are, respectively, the mean stellar metal-mass fraction for the target system and the Sun. The value of Δα = 35 suggested by Yan et al. (2020) in their Eq. (9) is now updated to Δα = 63 due to a different Z_⊙ assumption (Yan et al. 2020 should also have applied 0.0142 but actually applied 0.02 in the code; this change negligibly affected the results). We note that this latest IGIMF formulation is different from the one applied in Fontanot et al. (2017), where the IMFs for low- and intermediate-mass stars are invariant. Evidence of a systematic change of α₁ and α₂ with metallicity was already noted by Kroupa (2001) and formulated in Kroupa (2002), Marks et al. (2012), and Recchi et al. (2014). Recent observations of massive ETGs also support a bottom-heavy IMF for galaxies with super-solar metallicity (Zhou et al. 2019; Smith 2020), in agreement with Eq. (4). This is an important factor determining the variable ‘overall SNIa realisation parameter’ introduced in Sect. 5.2.2, below.

The variation of the power-law indices of the IMF for massive stars is determined empirically by Marks et al. (2012) as

$$\begin{array}{c}\hfill {\alpha }_3=\{\begin{array}{cc}2.3,\hfill & x<-0.87,\hfill \\ -0.41x+1.94,\hfill & x>-0.87.\hfill \end{array}\end{array}$$(5)

The parameter x in Eq. (5) depends on the metallicity, [Z/X] = log₁₀(Z/X)−log₁₀(Z_⊙/X_⊙), where Z and X are the initial metal and hydrogen mass fractions in stars, respectively. We assume the initial hydrogen mass fraction of any star, X, is the same as the Sun, i.e., X ≈ X_⊙; therefore, [Z/X]≈log₁₀(Z/Z_⊙) = [Z]. In addition, the parameter x in Eq. (5) depends on the core density (final stars plus residual gas assuming a star formation efficiency of 1/3), ρ_cl, of the molecular cloud core that forms the embedded star cluster,

$$\begin{array}{c}\hfill x=-0.14[Z]+0.99{log}_{10}({\rho }_{\mathrm{cl}}/{10}^6),\end{array}$$(6)

with

$$\begin{array}{c}\hfill {log}_{10}{\rho }_{\mathrm{cl}}=0.61{log}_{10}{M}_{\mathrm{ecl}}+2.85,\end{array}$$(7)

as is explained and applied in Yan et al. (2017, their Eq. (7)), with ρ_cl and M_ecl being in astronomical units. For example, an embedded star cluster weighing M_ecl = 10³ M_⊙ in stars has ρ_cl = 4.79 × 10⁴ M_⊙ pc⁻³.

The distribution of all embedded star clusters formed in a galaxy within the δt = 10 Myr star-formation epoch is approximated by a power-law ECMF (Gieles et al. 2006; Lieberz & Kroupa 2017):

$$\begin{array}{c}\hfill {\xi }_{\mathrm{ecl}}=\mathrm{d}{N}_{\mathrm{ecl}}/\mathrm{d}{M}_{\mathrm{ecl}}=k_{\mathrm{ecl}}{M}_{\mathrm{ecl}}^{-\beta },5M_{\odot }⩽M_{\mathrm{ecl}}<M_{\mathrm{ecl},\max }.\end{array}$$(8)

The total mass formed in stars during the δt epoch is given by

$$\begin{array}{c}\hfill M_{\mathrm{tot}}={\displaystyle {\int }_{M_{\mathrm{ecl},\min }}^{M_{\mathrm{ecl},\max }}}{M}_{\mathrm{ecl}}{\xi }_{\mathrm{ecl}}(M_{\mathrm{ecl}})\mathrm{d}{M}_{\mathrm{ecl}}={\overline{\psi }}_{\delta t}\delta t,\end{array}$$(9)

and the mass limit for the embedded clusters, M_ecl, max, is given by the optimal sampling law

$$\begin{array}{c}\hfill 1={\displaystyle {\int }_{M_{\mathrm{ecl},\max }}^{{10}^9{\mathrm{M}}_{\odot }}}{k}_{\mathrm{ecl}}{M}_{\mathrm{ecl}}^{-\beta }\mathrm{d}{M}_{\mathrm{ecl}},\end{array}$$(10)

with 10⁹ M_⊙ being about the mass of the most-massive ultra-compact-dwarf galaxy. Solving the two equations above yields k_ecl and M_ecl, max. We note that $M_{\mathrm{ecl},\max }=M_{\mathrm{ecl},\max }({\overline{\psi }}_{\delta t}),$ which is a relation that is consistent with the extragalactic most-massive-very-young star clusters versus galaxy-wide SFR data (Weidner et al. 2004; Schulz et al. 2015; Randriamanakoto et al. 2013). The slope of the ECMF, β, depends on the mean galaxy-wide SFR over the δt time period, ${\overline{\psi }}_{\delta t}$ (Weidner et al. 2013b):

$$\begin{array}{c}\hfill \beta =-0.106{log}_{10}{\overline{\psi }}_{\delta t}+2.\end{array}$$(11)

The environment-dependent gwIMF (i.e. the IGIMF) follows by adding up all the stars formed in the δt epoch:

$$\begin{array}{cc}& \xi (t)=\mathrm{d}{N}_{\star }/\mathrm{d}m={\displaystyle {\int }_{5M_{\odot }}^{M_{\mathrm{ecl},\max }}}{\xi }_{\star }(m,M_{\mathrm{ecl}},[Z/X])\hfill \\ \hfill & {\xi }_{\mathrm{ecl}}(M_{\mathrm{ecl}},{\overline{\psi }}_{\delta t})\mathrm{d}{M}_{\mathrm{ecl}}.\hfill \end{array}$$(12)

The environment dependence comes in through the SFR, which is a function of time, ${\overline{\psi }}_{\delta t}={\overline{\psi }}_{\delta t}(t)$, and through the stellar IMF and the ECMF being functions of the time-changing metallicity and SFR.

In summary, Eqs. (1), (4), (5), (6), (8), and (11) are empirical formulations. The integral of Eqs. (3) and (10) being 1 is an ansatz of the optimal sampling theory (Kroupa et al. 2013, their Sect. 2.2 and Schulz et al. 2015). The existence of the m_max = m_max(M_ecl) relation is supported by observed m_max–M_ecl data as is demonstrated in Yan et al. (2017) and Oh & Kroupa (2018). The above formulation is thus consistent with the observed constraints on the IMF in resolved star clusters and the Galactic field, as well as with the above-mentioned extragalactic data, and thus it fulfils this necessary condition.

2.2. Number of type Ia supernovae

The number of SNIa depends on the number of possible SNIa progenitors, i.e. on the number of stars with a mass between about 3 and 8 M_⊙ (see below). Therefore, the number of SNIa per unit mass of stars formed is a function of the IMF or the gwIMF. In previous research, the SNIa incidence has been taken to correlate with the total number of stars. Such an analytical formulation, first developed by Greggio (2005), assumes the canonical universal IMF and cannot be directly applied when assuming a different or evolving gwIMF.

Here, we adjust the Greggio (2005) formulation in order to account for the variable gwIMF. The total number of SNIa explosions for a simple stellar population (SSP, stars formed at the same time) after t years of the birth of the SSP per unit stellar mass of the SSP (i.e. the time-integrated number of SNIa per stellar mass formed until time t) is

$$\begin{array}{c}\hfill n_{\mathrm{Ia}}(t,\xi ,{\overline{\psi }}_{\delta t})=N_{\mathrm{Ia}}(\xi ,{\overline{\psi }}_{\delta t}){\displaystyle {\int }_0^t}{f}_{\mathrm{delay}}(t)\mathrm{d}t,\end{array}$$(13)

where N_Ia is the total number of SNIa per unit mass of stars formed in the SSP. The integral from t = 0 to ∞ is 1; therefore, N_Ia = n_Ia(t = ∞)≈n_Ia(t = 10 Gyr). It depends on the gwIMF of the stellar population, ξ, and on the star formation environment represented by the SFR, ${\overline{\psi }}_{\delta t}$. Here, f_delay is the DTD of the SNIa events for which we adopt the empirical power-law function given by Maoz & Mannucci (2012):

$$\begin{array}{c}\hfill f_{\mathrm{delay}}(t)=\{\begin{array}{cc}0,\hfill & t⩽40\mathrm{Myr},\hfill \\ k·t^{-1},\hfill & t>40\mathrm{Myr},\hfill \end{array}\end{array}$$(14)

where f_delay(t) is the fraction of exploded SNIa for a SSP with age t. The normalisation factor k is determined by the condition ${\displaystyle {\int }_{t=0}^{\infty }}{f}_{\mathrm{delay}}(t)\mathrm{d}t=1$.

In the present context, in which a systematically varying IMF/gwIMF is applied, the total number of SNIa per unit stellar mass of the SSP, N_Ia, (or SNIa production efficiency), is calculated as

$$\begin{array}{c}\hfill N_{\mathrm{Ia}}(\xi ,{\overline{\psi }}_{\delta t})=\frac{n_{3,8}(\xi )}{M_{0.08,150}(\xi )}·B_{\mathrm{bin}}({\overline{\psi }}_{\delta t})·\frac{n_{3,8}(\xi )}{n_{0.08,150}(\xi )}·C_{\mathrm{Ia}}({\overline{\psi }}_{\delta t}),\end{array}$$(15)

where n is the number of stars within the mass range given by its subscript (in the unit of M_⊙). Similarly, M is the mass of stars in a given mass range indicated by its subscript. The first term represents the number of stars that are possible to become SNIa progenitors (the primary star in a binary system) per unit stellar mass of the SSP. The fraction of stars in the initial mass range 3 to 8 M_⊙ that eventually explode as SNe Ia consists of the following three terms. The second term, B_bin, denotes the binary fraction of the SSP which can be assumed to be a constant or a function of the environment represented by the galaxy-wide SFR, ${\overline{\psi }}_{\delta t}$. The third term then represents the likelihood of a binary system to have the companion star also in the suitable stellar mass range thus being a potential SNIa progenitor system (a binary white dwarf, WD). Both the first and third terms depend on the IMF of each star cluster since it is extremely unlikely for the binaries to form outside a star cluster. However, since the gwIMF is the mass-weighted sum of all IMFs of individual star clusters in a galaxy, we apply gwIMF to Eq. (15) to simplify the calculation when considering a galaxy. Finally, C_Ia is the realisation probability of an SNIa explosion for the potential SNIa progenitor system; that is, the fraction of the above-selected binary systems that are able to give rise to an SNIa explosion. Again, it was assumed to be a constant in previous chemical evolution studies but should, in principle, be a function of the star formation environment, because it depends on the distribution of semi-major axes and the mass ratios of the binary WD systems. When the galactic SFR is high, more massive clusters form, as calculated by the IGIMF theory (Eqs. (8)–(11)). Massive clusters lead to substantial long-term dynamical processing of the binary-star population (Heggie 1975; Marks et al. 2011). The soft initial binaries, which have binary-orbital velocities slower than the velocity dispersion in the cluster, are disrupted, while hard binaries (with binary-orbital velocities that are faster than the velocity dispersion in the cluster) are not disrupted but sink to the cluster centre through dynamical friction. This is also true for WD–stellar or WD–WD binaries. There, they harden through encounters that also increase their binary-orbital eccentricities. Few would be expected to be ejected in dynamical strong encounters, and these are likely to merge and explode as SNIa. Many will not be able to escape from their massive clusters and will continue to encounter other stars and remnants to also ultimately explode as SNIa (Shara & Hurley 2002). Massive star clusters (and thus high SFRs) are thus expected to correlate with enhanced SNIa events. This appears to be born out, not only by the here-derived Eq. (15), but also by Friedmann & Maoz (2018), Freundlich & Maoz (2021) finding evidence for a significantly enhanced occurrence of SNIa in galaxy clusters.

Regarding the IMF-dependent terms in Eq. (15), we assume that the stars with a mass between 3 and 8 M_⊙ have an equal probability of leading to an SNIa event, while stars with other masses have zero probability. Stars of 1.5 M_⊙, for example, may also become WDs in reality, but, compared to the more massive stars, they have a significantly longer lifetime and almost no chance of being primary stars in binary systems with a total mass above 3 M_⊙, which is considered as the lower binary system mass limit for SNIa events (Matteucci & Recchi 2001). As a simplification, we apply the same mass limit (3–8 M_⊙) for the primary and the secondary star. These assumed mass limits insignificantly affect the dependency of the SNIa number on the IMF. Changing both SNIa progenitor mass limits for the primary and the secondary stars to 2–8 M_⊙, for example, would slightly increase the SNIa production efficiency for the low-SFR galaxies (shown by Fig. A.1), but it would not affect our conclusions. Having different mass limits for the primary and the secondary stars would give similar results (Fig. A.2).

Therefore, for the present study, the only free parameter in Eqs. (13)–(15) is the product $B_{\mathrm{bin}}({\overline{\psi }}_{\delta t})·C_{\mathrm{Ia}}({\overline{\psi }}_{\delta t})$, i.e., the overall SNIa realisation parameter. We set this parameter as a constant in our ‘fiducial SNIa model’. It could become higher for more massive elliptical galaxies because the stars are formed, on average, in a more crowded and metal-rich environment. This leads to a higher fraction of hard binaries, larger asymptotic giant branch (AGB) star radii, and a more efficient accretion onto the pre-supernova white dwarf (Hachisu et al. 1996). We note that the possible variation of the overall SNIa realisation parameter could be due to a variation of both SFR and metallicity but these two dependencies degenerate in the presented description (Eqs. (15) and (20), below), which is only a function of SFR because both the peak SFR and mean stellar metallicity of E galaxies depend on their mass monotonically.

In order to calibrate the overall SNIa realisation parameter $B_{\mathrm{bin}}({\overline{\psi }}_{\delta t})·C_{\mathrm{Ia}}({\overline{\psi }}_{\delta t})$, we start from a given galaxy-wide 10 Myr averaged SFR, ψ₀ = 1 M_⊙/yr. The value of B_bin(ψ₀)⋅C_Ia(ψ₀) can be determined by the observed number of SNIa events in nearby stellar systems assuming they have the canonical gwIMF and a Galactic SFR of ψ₀. The local SNIa incidence per unit mass of stars formed is determined by the observation of different systems with different SFR. Therefore, the calibration assuming a mean SFR of the local universe of ψ₀ = 1 M_⊙/yr may not be exactly correct, but it has led to successful reproduction of the work of other groups (Yan et al. 2019a, 2020). Following Maoz & Mannucci (2012), we set

$$\begin{array}{c}\hfill n_{\mathrm{Ia}}(t=10\mathrm{Gyr},{\xi }_{\mathrm{canonical}},{\psi }_0)=0.0022/M_{\odot },\end{array}$$(16)

where ξ_canonical denotes that the relation holds when adopting the canonical Kroupa IMF (Kroupa 2001, i.e. Eq. (1) with α₁ = 1.3 and α₂ = α₃ = 2.3).

Under the fiducial SNIa model with a constant overall SNIa realisation parameter, $B_{\mathrm{bin}}({\overline{\psi }}_{\delta t})·C_{\mathrm{Ia}}({\overline{\psi }}_{\delta t})=B_{\mathrm{bin}}({\psi }_0)·C_{\mathrm{Ia}}({\psi }_0)$, the above formulation in Eq. (15) is completely defined and quantifies how n_Ia(t = 10 Gyr, ξ) varies given different gwIMFs, ξ(m), and SFRs, ${\overline{\psi }}_{\delta t}$. We note that we drop ${\overline{\psi }}_{\delta t}$ in $n_{\mathrm{Ia}}(t,\xi ,{\overline{\psi }}_{\delta t})$ because this dependency is contained in the gwIMF, ξ, by the latter being a function of ${\overline{\psi }}_{\delta t}$. As an example, Fig. 1 demonstrates how n_Ia(t = 10 Gyr, ξ) varies with the IMF power-law index, α₃. As can be seen, the SNIa production efficiency is maximised for α₃ ≈ 1.8 and suppressed for extremely top-heavy IMF cases. Such stellar populations may not be impossible: α₃ <  1 for young globular clusters with masses larger than about 10⁸ M_⊙ according to Marks et al. (2012, Fig. 2), while α₃ >  4 for the gwIMF in dwarf galaxies with ${\overline{\psi }}_{\delta t}<{10}^{-4}{M}_{\odot }$/yr according to Yan et al. (2017, Fig. 6).

Fig. 1. Total number of SNIa explosions for a SSP after 10 Gyr of the birth of the SSP per unit stellar mass of the SSP predicted by Eqs. (13)–(16), assuming constant parameter Bbin and CIa and an IMF defined in Eq. (1), with α1 = 1.3, α2 = 2.3, as a function of α3. A smaller α3 indicates a flatter IMF/gwIMF for the massive stars, that is, a top-heavy IMF/gwIMF. Very large α3 mean that essentially no stars above 1 M⊙ are formed. The function is calibrated by the vertical and horizontal dashed lines indicating the canonical IMF/gwIMF slope, α3 = 2.3, and the value nIa(t = 10 Gyr, ξcanonical) = 0.0022/M⊙ suggested by Maoz & Mannucci (2012) for the local Universe, respectively.

Since the gwIMF predicted by the IGIMF theory is a function of the galaxy-wide SFR and metallicity (see Eq. (12)), Fig. 2 shows how n_Ia(t = 10 Gyr) changes as a function of these two parameters for our fiducial SNIa model. The resulting n_Ia(t = 10 Gyr) for the metal-poor galaxies can be much higher than the values shown in Fig. 1 because they have a bottom-light gwIMF according to Eq. (4). The galaxy-wide SFR also affects the number of SNIa per stellar mass formed because of Eq. (5). The calculated n_Ia(t = 10 Gyr) variation due to different gwIMF agrees with the observational constraints for SNIa surveys (listed in Table 1) and could be a natural explanation for the difference in these observations. This is not the case for the older version of the IGIMF formulations summarised in Jeřábková et al. (2018, Table 3) as is demonstrated by Figs. A.3 to A.5 of this paper.

Fig. 2. Ten Gyr time-integrated number of SNIa per unit stellar mass formed for the fiducial SNIa model calculated according to Eqs. (13)–(16). The fiducial SNIa model assumes a constant overall SNIa realisation parameter and the gwIMF as given by the IGIMF theory (Eq. (12)) as a function of the galaxy-wide SFR, ${\overline{\psi }}_{\delta t}$, and metallicity, [Z/X]. The black line is the relation for [Z/X] = 0. Other lines with different colours represent different values of [Z/X], as indicated by the white stripes on the colour map on the right: [Z/X] = 0.3, 0.2, 0.1, −0.1, −0.2, −0.3, −0.4, −0.5, −1, −2, −4, and −6. The black horizontal dotted line represents the canonical nIa(t = 10 Gyr, ξcanonical) value of 0.0022/M⊙ (Maoz & Mannucci 2012). The green and red horizontal dashed lines indicate observational constraints on the nIa(t = 10 Gyr) for SNIa surveys up to a certain redshift and in galaxy clusters, respectively. The shaded regions represent the uncertainty ranges of the horizontal dashed lines. References are given in Table 1.

In addition, we divide the estimations of the SNIa production efficiency in Table 1 into two groups, those targeted at galaxy clusters and volume-limited and field galaxies (plotted as the red and green shaded regions, respectively, in Fig. 2). These results suggest that galaxy clusters have a higher SNIa production efficiency. The n_Ia(t = 10 Gyr) of our predicted high-SFR metal-rich galaxies is lower than these estimates. As mentioned above, it is reasonable to consider that the overall SNIa realisation parameter is higher for the galaxies that have experienced more intense star formation in a galaxy cluster environment. We discuss such scenarios in Sect. 5.2.2.

2.3. The galaxy chemical evolution model

Motivated by the monolithic collapse formation scenario of E galaxies (see Sect. 1) and following Thomas et al. (1999), the single-zone, closed-box, galaxy chemical evolution modelling computer programme, GalIMF, described in Yan et al. (2019b), is applied. GalIMF assumes that all the gas is always well mixed (i.e. single-zone) and that the stars inject a certain amount of mass of each element into the gas when they exhaust their lifetime (non-instantaneous recycling of metals). The smallest time step in our galaxy-wide calculation is δt = 10 Myr because this is the timescale required for a galaxy to fully populate the ECMF (Yan et al. 2017), and the gwIMF cannot be defined over shorter timescales (see also Sect. 2.1).

The stellar yield table assumed here is the same as in Yan et al. (2020). That is, the stellar lifetime, remnant mass, and metal yield for low-mass stars are adopted from Marigo (2001). For massive stars, they are adopted from Kobayashi et al. (2006). The stellar yield is given according to the initial mass and metallicity of a star (cf. Yan et al. 2019b, Figs. 3 and 4). The yield table is first interpolated in the dimension of stellar mass, while for stars with a mass higher or lower than the mass range provided in the table, a fixed value is applied; this is respectively, the value for the most massive or least massive star (cf. Yan et al. 2020, Fig. 8). Then, the table is interpolated again in the dimension of the stellar initial metallicity. For metallicities higher or lower than the range provided in the table, the value from the largest or lowest metallicity in the table is applied, respectively. We note that the yield table for massive stars is different from our previous work on constraining the SFT of elliptical galaxies (Yan et al. 2019a), which follows the prescription of Thomas et al. (1999). This introduces a systematic difference but preserves the general trend of the τ_SF–M_dyn relation and does not change our conclusions, as is explained in Sect. 6.2.

3. Observational constrains

The galaxy models need to be compared with the observations to constrain the parameters of the galaxy models such as the SFTs. Here, we apply the observational [Z/X] and [Mg/Fe] constraints shown by the thick solid and dashed lines in Fig. 3. The [Z/X]–M_dyn relation is constrained by Arrigoni et al. (2010), and the mean and the standard deviations of the relation are calculated using the procedure detailed in Yan et al. (2019a, Sect. 2.1). The comparison between the galaxy dynamical mass derived from observation and the dynamical mass from our model is detailed in Yan et al. (2019a, Sect. 2.3).

Fig. 3. [Z/X], [Mg/Fe], and age of the E galaxies as a function of their dynamical masses. The grey circles are observed galaxies from Arrigoni et al. (2010), with the typical uncertainty shown at the lower right corner of each panel. The solid and dashed lines are, respectively, the mean and standard deviations of the observational constraints. The lines in the upper panel are given by Arrigoni et al. (2010) and smoothed by Yan et al. (2019a). The mean relation for the middle panel follows Thomas et al. (2005), while the standard deviation of the relation follows Liu et al. (2016). The points and triangles are the resulting Mdyn, [Z/X], and [Mg/Fe] values of our fiducial galaxy models adopting the IGIMF theory and a constant overall SNIa realisation parameter defined in Sect. 2.2. The triangles are models with ${\log }_{10}({\overline{\psi }}_{\delta t})=0,$ while the points are models with other SFRs. The red and the blue dots are models with a SFT τSF ≥ 1 Gyr and τSF <  1 Gyr, respectively. The darker enhanced colours represent models with a higher gconv (defined in Sect. 4).

We apply the [Mg/Fe]–M_dyn relation given by Thomas et al. (2005, Eq. (3)) but a larger [Mg/Fe] scatter for galaxies with a mass smaller than 10^9.218 M_⊙, as is suggested by Liu et al. (2016, Table 1), which was calculated with the galaxy-mass–central-velocity-dispersion relation given using Thomas et al. (2005, Eq. (2)). The [Mg/Fe]–M_dyn relation given by Thomas et al. (2005) has a higher [Mg/Fe] value of about 0.05 dex than the relation given by the galaxies in Arrigoni et al. (2010), Table B.1 (see the solid line and the grey circles in Fig. 3). We adopt the Thomas et al. (2005) relation, because it is supported by the more recent and complete study by Liu et al. (2016), and this choice is consistent with our previous study (Yan et al. 2019a). The standard deviation of [Mg/Fe], σ_Mg/Fe, is 0.126 for galaxies with M_dyn <  10^9.018 M_⊙ and 0.064 for galaxies with M_dyn >  10^9.418 M_⊙ and a linear transition in between.

The age of the observed galaxies given in Arrigoni et al. (2010, Table B.1) has a large scatter: $6.3_{-2.7}^{+4.7}$ Gyr. There is no strong correlation between the age and mass or abundance of the galaxies as is shown by the lower panel of Fig. 3, although Thomas et al. (2005) calculated that massive ETGs and ETGs in high-density environments form earlier statistically, which is known as ‘archaeological downsizing’. We neglect the age difference of these galaxies and evolve all our modelled galaxies to 14 Gyr for simplicity. Assuming a different galaxy age for each galaxy or all the galaxies would not affect our conclusions because the modelled galaxies are assumed to have quenched their star formation completely. The properties of their stellar population would not have any significant change after the quenching.

These observational relations are derived from the ETG database, but since elliptical and lenticular galaxies follow similar [Z/X]–M_dyn and [Mg/Fe]–M_dyn relations (cf. Gallazzi et al. 2005; Zhu et al. 2010; Camps-Fariña et al. 2021), the difference in galaxy type is unlikely to affect our conclusions, although a follow-up study distinguishing galaxy types would be valuable.

4. Method

Here, we explain how the GalIMF code is used to constrain the SFTs of galaxies following the same method as in Yan et al. (2019a). In short, there are three input parameters that describe the gas supply and SFH of a galaxy and there are three outputs after running the code that describe the mass and metal abundance of this galaxy after 14 Gyr. The output is compared with the observational values to evaluate the likelihood of the input parameters.

The chemical evolution of a set of different galaxies with different combinations of input parameters is calculated. Following Yan et al. (2019b) and Yan et al. (2019a), the SFHs of the E galaxies are assumed to have a constant SFR, ${\overline{\psi }}_{\delta t}$, over a certain SFT, τ_SF. A low level of star formation has been found as a commonality in ‘quenched’ galaxies (de Lorenzo-Cáceres et al. 2020). This could be the cause of the large scatter of the observed abundance, but it does not have a significant effect on our results because we focus on the smoothed mean relation of galaxies. In addition, Salvador-Rusiñol et al. (2020) find only a percent of stellar mass to have been added to ETGs during the past 2 Gyr (see also Kroupa et al. 2020b for an in-depth discussion of ETG assembly in connection to the formation of supermassive black holes). The third input parameter is the gas-conversion parameter, g_conv, defined as the ratio between the accumulated stellar mass (total mass of stars ever formed in a galaxy) and initial-gas mass (total gas supply since the model assumes no galactic gas inflow or outflow). By adjusting the amount of initial gas supply for the same SFH, different g_conv values are tested. The input parameters for our fiducial SNIa model take the values of: g_conv= 0.05, 0.1, 0.2, 0.4, or 0.8, ${\log }_{10}({\overline{\psi }}_{\delta t}[M_{\odot }/\mathrm{yr}])=-1.0$, 0, 1.0, …, or 4.0, and τ_SF = 0.05, 0.2, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.25, 2.5, 2.75, 3, 3.5, or 4 [Gyr]. In total, 5 × 6 × 15 = 450 galaxy models are calculated. The other models discussed in Sect. 5.2.2 are similar but focus on different parameter ranges.

All galaxies are evolved to 14 Gyr for comparison with the observations, because the evolution of the mean stellar properties later than 1 Gyr after cessation of star formation for a galaxy is negligible. Therefore, the applied time steps for each model are 0.01, 0.02, …, τ_SF − 0.01, τ_SF, and 14 [Gyr]. The resulting galactic stellar plus remnant mass, M_dyn, mean stellar metallicity, and mean stellar [Mg/Fe] for the model galaxies with every combination of the input parameters are calculated and plotted in Fig. 3. Then, the results are interpolated in 3D for the three parameters (g_conv, ${\overline{\psi }}_{\delta t}$, and τ_SF) and compared with the mean observational values of galaxies with the same mass.

The goodness of the fit for each galaxy model is calculated as

$$\begin{array}{cc}& p(g_{\mathrm{conv}},{\tau }_{\mathrm{SF}})=p_{\mathrm{Z}/\mathrm{X}}\times p_{\mathrm{Mg}/\mathrm{Fe}},\hfill \\ \hfill & p_{\mathrm{Z}/\mathrm{X}}=\hfill \\ \hfill & 1-\mathrm{erf}[({[\mathrm{Z}/\mathrm{X}]}_{\mathrm{obs}}-{[\mathrm{Z}/\mathrm{X}]}_{\mathrm{mod}}(g_{\mathrm{conv}},{\tau }_{\mathrm{SF}}))/{\sigma }_{\mathrm{Z}/\mathrm{X}}/\sqrt{2}],\hfill \\ \hfill & p_{\mathrm{Mg}/\mathrm{Fe}}=\hfill \\ \hfill & 1-\mathrm{erf}[({[\mathrm{Mg}/\mathrm{Fe}]}_{\mathrm{obs}}-{[\mathrm{Mg}/\mathrm{Fe}]}_{\mathrm{mod}}(g_{\mathrm{conv}},{\tau }_{\mathrm{SF}}))/{\sigma }_{\mathrm{Mg}/\mathrm{Fe}}/\sqrt{2}],\hfill \end{array}$$(17)

where [Z/X]_obs and [Mg/Fe]_obs are the mean observational values and the σ stands for the standard deviation of the element abundance ratios for different galaxies for a given mass (see Yan et al. 2019a). For each M_dyn and τ_SF, the highest p value with any g_conv is calculated and shown in Figs. 4 and 7. Then, the τ_SF with the highest possible p for a given M_dyn indicates the mean SFT for the elliptical galaxies with the mass M_dyn.

Fig. 4. τSF as a function of Mdyn. The colour bar indicates the goodness of the fit, p(gconv, τSF) (Eq. (17)), for our fiducial SNIa model. For a given mass and SFT, τSF, any star-to-gas-mass fraction, gconv, is allowed. The most likely SFT for a given galaxy mass is indicated by the yellow ridge line. The black dotted curve (Eq. (18)) is a power-law fit of the yellow ridge line. The red dashed curve and the red dash-dotted curve (Eq. (19)) represent the SFT constrained by SPS studies by de La Rosa et al. (2011) and McDermid et al. (2015), respectively. As is introduced in Sect. 1, The SFT suggested by chemical evolution studies assuming the invariant canonical gwIMF performed by Pipino & Matteucci (2004) and Thomas et al. (2005) are shown by the blue triangles and blue dotted curve, respectively. The lower and upper yellow dotted curves correspond, respectively, to polynomial regressions of the yellow ridge lines in Figs. 5 and 7 of Yan et al. (2019a), derived using invariant canonical gwIMF. The mass ranges of the coloured curves are limited by the mass ranges of the galaxy data sets they are based on.

Finally, we modify the default assumptions in the GalIMF code, including the IMF and SNIa formulation to study their effects on the results. That is, we calculate a galaxy model grid for each set of IMF/SNIa assumptions.

5. Results

In Sect. 5.1, the results obtained by applying an invariant canonical IMF for all galaxies are documented. This serves as a benchmark model set for the more realistic case documented in Sect. 5.2, where the galaxy-wide IMF changes with varying conditions, as is calculable self-consistently (from the time-dependent SFR and metallicity) using the IGIMF theory.

5.1. Canonical IMF

5.1.1. Canonical IMF and invariant SNIa production efficiency

We explored the possible τ_SF–M_dyn relations in Yan et al. (2019a) assuming the invariant canonical gwIMF. The results are plotted in Fig. 4 as dotted yellow lines along with the estimates given by Pipino & Matteucci (2004) and Thomas et al. (2005). The blue dotted line is from Thomas et al. (2005). It assumes a constant total-star-to-gas-mass fraction, g_conv, in their galaxy chemical evolution model and takes into account only the [Mg/Fe] constraints. On the other hand, Pipino & Matteucci (2004, blue triangles) and Yan et al. (2019a, yellow dotted lines) adopt a flexible total-star-to-gas-mass fraction and simultaneously fit the stellar [Mg/Fe] and the metallicity, resulting in steeper relations.

These results, constrained by the chemical abundance of the galaxies, suggest shorter SFTs than the timescale estimated by the SPS studies (de La Rosa et al. 2011; McDermid et al. 2015). As is mentioned in Sect. 1, the cosmological hydrodynamical simulations of galaxy formation taking into account a realistic gas recycling time also have difficulty in reproducing the stellar [Mg/Fe] and metallicity of the most massive galaxies when assuming the canonical universal IMF. The short SFT for the most massive E galaxies suggested by their high [Mg/Fe] value may be relaxed if they have a top-heavy IMF or there are fewer SNIa in these galaxies.

5.1.2. Canonical IMF and variable SNIa production efficiency

To have a longer SFT (more SNIa explosion and iron produced during the star formation) for the same observed [Mg/Fe] of massive ETGs, a lower SNIa production efficiency would be required assuming an invariant IMF. This scenario contradicts the expectation that massive galaxies that formed more massive clusters should have formed more hard WD–WD and WD–stellar binaries and, therefore, should have more SNIa explosions (Shara & Hurley 2002). This led us to consider the observationally motivated variation of the gwIMF, as described in Sect. 5.2.

5.2. IGIMF theory

It is always possible to fit the observed galactic abundance if one is allowed to fine-tune the gwIMF formulation with no constraints. However, one needs to follow the observed IMF variation of the resolved stellar populations; this means that the gwIMF is not arbitrary. The here-applied IGIMF formulation (Sect. 2.1) is determined by the empirically deduced IMF and ECMF variations, and the calculated gwIMF automatically emerges and happens to be consistent with the observed gwIMF of vastly different systems.

5.2.1. IGIMF theory and the fiducial SNIa model

For the fiducial SNIa model, we assume the parameters B_bin and C_Ia in Eq. (15) to be invariant, i.e., SNIa production efficiency depends only on the gwIMF but not the star formation environment. Then, there is only one unknown overall SNIa realisation parameter in Eq. (13), i.e., B_bin(ψ₀)⋅C_Ia(ψ₀), which can be determined by Eq. (16). Under this assumption, the number of SNIa is only affected by the shape of the gwIMF (through the ξ related parameters in Eq. (15)), but it is not directly impacted by the physical condition (density and metallicity) of the star formation cloud.

The resulting [Z/X]–M_dyn and [Mg/Fe]–M_dyn relations of our galaxy models adopting the IGIMF theory with different combinations of input ${\overline{\psi }}_{\delta t}$, τ_SF, and g_conv (detailed in Sect. 4) are plotted in Fig. 3. Comparing the results with different input parameters shown by the different symbols (see e.g., the middle panel, the six nearly vertical sequences of galaxy model results, from left to right, are models with ${\log }_{10}({\overline{\psi }}_{\delta t})=-1.0$, 0, 1.0, …, or 4.0), it appears that the SFR is the main factor for establishing the resulting galaxy mass. Darker points have a higher g_conv, i.e., more stellar mass formed per unit of primordial gas provided. Thus, the g_conv value dominates the [Z/X] variation as is shown by the top panel. The change of SFT from 0.05 to 4 Gyr results in the galaxy model groups from the top to the bottom in the middle panel of Fig. 3, where the models with τ_SF ≥ 1 Gyr are represented by red symbols, and this change dominates the [Mg/Fe] variation.

The likelihoods of different SFTs for each galaxy mass are calculated according to Eq. (17) and shown in Fig. 4 by the colour map. The best-fit SFT values (shown as the yellow ridge line) are longer for more massive galaxies because they have more top-heavy gwIMFs, leading to a larger production of magnesium. A power-law function representing the yellow ridge line in Fig. 4, shown by the black dotted line, is given by

$$\begin{array}{c}\hfill {\tau }_{\mathrm{SF}}[\mathrm{Gyr}]=0.003·{(M_{\mathrm{dyn}}/M_{\odot })}^{0.3}\end{array}$$(18)

for a set of galaxies with a M_dyn above 10^8.3 M_⊙, although there is only one galaxy in the data set with a mass below 10⁹ M_⊙ (see Fig. 3). This resulting relation is not consistent with the relations suggested by the independent SPS studies shown by the red dashed and red dash-dotted lines given by de La Rosa et al. (2011) and McDermid et al. (2015), respectively. For example, in combination with Eq. (3) of Thomas et al. (2005), McDermid et al. (2015, Eq. (3)) suggests

$$\begin{array}{c}\hfill {\tau }_{\mathrm{SF}}[\mathrm{Gyr}]=49·{(M_{\mathrm{dyn}}/M_{\odot })}^{-0.14}\end{array}$$(19)

for a set of galaxies with a M_dyn between about 10^9.8 and 10¹² M_⊙. This relation is shown as the red dash-dotted curve in Figs. 4, 7, and 8. However, as is shown next, taking into account that more massive E galaxies are likely to have formed more massive and dense starburst clusters, which produce a larger fraction of SNIa-progenitor binaries per star (Shara & Hurley 2002), it is possible to obtain consistent results.

5.2.2. IGIMF theory and a boosted SNIa production efficiency for high-SFR galaxies

We define an SNIa realisation re-normalisation parameter, ${\kappa }_{\mathrm{Ia}}({\overline{\psi }}_{\delta t})$, which represents the variation of the overall SNIa realisation parameter (the terms independent of the IMF in Eq. (15)) as a function of the galaxy-wide SFR:

$$\begin{array}{c}\hfill {\kappa }_{\mathrm{Ia}}({\overline{\psi }}_{\delta t})=\frac{B_{\mathrm{bin}}({\overline{\psi }}_{\delta t})·C_{\mathrm{Ia}}({\overline{\psi }}_{\delta t})}{B_{\mathrm{bin}}({\psi }_0)·C_{\mathrm{Ia}}({\psi }_0)},\end{array}$$(20)

where the constant B_bin(ψ₀)⋅C_Ia(ψ₀) is calibrated by reproducing this value for a single stellar population with the canonical Kroupa (2001) IMF, while the variable $B_{\mathrm{bin}}({\overline{\psi }}_{\delta t})·C_{\mathrm{Ia}}({\overline{\psi }}_{\delta t})$ can become larger in massive galaxies, as mentioned in Sect. 2.2.

Through trial and error, we find that assuming the error function

$$\begin{array}{c}\hfill {\kappa }_{\mathrm{Ia}}({\overline{\psi }}_{\delta t})=1.75+0.75·\mathrm{erf}[{\log }_{10}({\overline{\psi }}_{\delta t})·1.25-2],\end{array}$$(21)

where erf stands for the Gauss error function (the solid curve in Fig. 5) and leads to a result (the yellow ridge line in Fig. 7) that roughly fits the τ_SF, SPS–M_dyn relation suggested by McDermid et al. (2015). The SNIa production efficiency is boosted by κ = 2.5 for the most massive ellipticals compared to the unchanged value of κ = 1 for the low-mass galaxies. This increases the number of SNIa per stellar mass formed for high-SFR galaxies, as is shown in Fig. 6 in comparison with Fig. 2, and therefore, it decreases their SFT estimation. The significantly increased SNIa production efficiency of the high-SFR galaxies shown in Fig. 6 is consistent with the observational findings that the SNIa production efficiency in galaxy clusters (red regions in Fig. 6) that host massive ETGs is only two to three times higher than for field galaxies (green regions in Fig. 6, see e.g., Friedmann & Maoz 2018; Freundlich & Maoz 2021). This is because the high-SFR galaxies have an overall metal-rich stellar population, leveraging down the n_Ia(t = 10 Gyr) value.

Fig. 5. SNIa realisation re-normalisation parameter, κIa defined in Eq. (20), as a function of the galaxy-wide SFR, ${\overline{\psi }}_{\delta t}$. Three different models are tested, and the corresponding best-fit SFT–galaxy-mass relations are shown in Figs. 4, 7, and B.1 for the fiducial model (dotted line), error function model (solid line), and power-law function model (dashed line), respectively. The error function model, formulated in Eq. (21), best reproduces the observed ETGs.

Fig. 6. Ten Gyr time-integrated number of SNIa per unit of stellar mass formed for the error function model. Same as Fig. 2, but here the calculation assumes the variable κIa–${\overline{\psi }}_{\delta t}$ relation defined by Eq. (21) and is shown by the solid curve in Fig. 5 instead of a constant overall SNIa realisation parameter (dotted line in Fig. 5).

Figure 7 shows the resulting τ_SF–M_dyn relation for the error function model. The E galaxies with masses below about 10^9.5 M_⊙ most likely have an SFT that is shorter for lower mass galaxies, following the same black dotted curve as in Fig. 4, i.e., Eq. (18). This could arise because low-mass galaxies have a shallower potential and lose their gas supply more easily. Although the SFT for individual dwarf ellipticals can have a large scatter (cf. Fig. 3), our result is in good agreement with the independent SFT estimations for the ultra-diffuse galaxy Dragonfly 44 (Haghi et al. 2019) when using the SPS method and for the ultra-faint dwarf galaxy Boötes I (Yan et al. 2020) and analysing the chemical enrichment history of the galaxy. Both find an SFT less than 1 Gyr. On the other hand, E galaxies with masses larger than about 10^9.5 M_⊙ have SFT values that are shorter with increasing M_dyn, following the relation suggested by McDermid et al. (2015), i.e., Eq. (19). These massive galaxies are likely regulated by the gas collapsing timescale. As a result, the highest SFT is reached for galaxies with about 10^9.5 M_⊙. We note that the galaxy mass reaching the peak SFT can be smaller if κ is lower than 1 for low-SFR galaxies. It is possible to further fine-tune the ${\kappa }_{\mathrm{Ia}}({\overline{\psi }}_{\delta t})$ function to adjust the resulting τ_SF–M_dyn relation, which we do not do considering the substantial computational cost. The short SFT for low-mass galaxies is concluded here due to a top-light gwIMF for the low-SFR galaxies. This is a distinct difference compared to the conclusions reached by Thomas et al. (2005) and Yan et al. (2019a) assuming the invariant canonical IMF and the closed-box chemical evolution model. However, the very top-light gwIMF of the present-day low-SFR galaxies (Lee et al. 2009; Pflamm-Altenburg et al. 2009; Yan et al. 2017, Fig. 6) may not represent the gwIMF of dwarf galaxies during their formation. Dwarf galaxies can be perturbed and have a more complex SFH than the monolithic collapse scenario assumed here (see discussions in Sect. 6.3). With a fluctuated SFR, the stellar populations of the dwarf galaxies are formed with a higher instantaneous SFR, ${\overline{\psi }}_{\delta t}$, than the mean SFR over the same SFT, leading to a less top-light gwIMF and a longer derived SFT (cf. Sect. 6.5). In other words, for an isolated E galaxy, the difference between its real SFT and the calculated SFT using its abundance constraints and a constant SFR assumption indicates how fluctuated its SFR is.

Fig. 7. τSF as a function of Mdyn for the error function model. Same as Fig. 4, but the results shown by the colour map adopt the variable ${\kappa }_{\mathrm{Ia}}({\overline{\psi }}_{\delta t})$ defined by Eq. (21) and shown by the solid curve in Fig. 5. The best-fit solutions (the yellow ridge line) for different Mdyn can be described by two relations, i.e., the black dotted curve (Eq. (18)) for Mdyn <  109.5 M⊙, and the red dash-dotted curve (Eq. (19)) for Mdyn >  109.8 M⊙.

In addition, we explore the possibility that the SNIa production efficiency for high SFR galaxies continues to increase following a power-law relation. In the case of the power-law function model shown by the dashed curve in Fig. 5, we adopt ${\kappa }_{\mathrm{Ia}}({\overline{\psi }}_{\delta t})=0.8+0.2·{\overline{\psi }}_{\delta t}^{0.38}$, which is similar to Eq. (21) out to ${\log }_{10}({\overline{\psi }}_{\delta t})\approx 2$. The best-fit SFTs for the massive galaxies are much shorter and do not fit with the SPS constraints. This is shown by Fig B.1.

We note again that applying the same method of tuning the ${\kappa }_{\mathrm{Ia}}({\overline{\psi }}_{\delta t})$ function under the assumption that the gwIMF is the invariant canonical IMF does not result in a meaningful solution. Assuming an invariant canonical IMF suggests shorter SFT values for more massive galaxies and/or longer SFTs for low-mass ellipticals; therefore, they require a lower κ_Ia for massive galaxies and/or larger κ_Ia for low-mass galaxies to fit the SFT estimated by the SPS studies. This contradicts the expectation that massive galaxies with a more intense star formation activity should have a greater chance of forming hard WD–WD and WD–stellar binary systems and SNIa events through the stellar-dynamical processing of binaries in massive star clusters, which cannot form in galaxies with low SFRs.

The SFT, SFR, gas-conversion parameter, and dynamical-to-stellar mass ratio of our best-fit models are shown in Fig. 8. We find that roughly the same fraction of gas is consumed in our error function SNIa model (Eq. (21)) for galaxies with all different masses (the panel for g_conv), being similar to the consumption fraction in embedded star clusters (Sect. 2.1). This indicates a universal star formation efficiency of about 1/3 for monolithically (free-fall collapse) formed systems. We note that the wiggles are caused by the interpolation of the results from a limited number of galaxy models. Other than the wiggles, the ratios of dynamical mass and stellar mass, M_dyn/M_*, are higher for more massive galaxies because they have a more top-heavy gwIMF and more stellar remnants. However, this result still needs to be translated to mass-to-light ratio, M_dyn/L_V, to compare with the observation. The difference between the M_dyn/L_V values of massive and low-mass galaxies is even larger than the difference between the M_dyn/M_* values. This is because massive galaxies are also more metal-rich and subsequently have a more bottom-heavy gwIMF, while the additional low-mass stars barely contribute to the V-band galaxy luminosity. A detailed discussion of this is given in Sect. 6.1.

Fig. 8. SFT, SFR, gas-conversion parameter, and dynamical-mass-to-stellar mass ratio of our best-fit model, from the top to the bottom panel, for the error function model. The solid black line in the top panel is the same as the yellow ridge line of Fig. 7. The black dotted and the red dash-dotted curves in the top panel follow the same relation as in Fig. 7, described by Eqs. (18) and (19), respectively. The thin horizontal lines indicate gconv values applied in our galaxy model (see Sect. 4), where gconv ≈ 1/3 is about the star formation efficiency found for embedded star clusters (Sect. 2.1).

6. Discussion

6.1. Mass-to-light ratio

We explore the V-band dynamical mass-to-light values, M_dyn/L_V, of the best-fit models described in Sect. 5.2.2. To simplify the calculation, we apply a simple-population approximation that assumes that all the stars in a galaxy have the same age and metallicity, where the gwIMF of the galaxy is calculated with the best-fit SFR (second panel of Fig. 8) and the observed mean stellar metallicity (upper panel of Fig. 3). The simple-population gwIMFs of a few galaxies with different masses are shown in Fig 9. The M_dyn/L_V value of a galaxy can be calculated for these gwIMFs for a given age. We apply t = 3.6, 6.3, and 11 Gyr corresponding to the mean and standard deviation of the estimated galaxy ages of the observed galaxy population (see Sect. 3). Since in the GalIMF code the luminosity is calculated using the main-sequence stellar luminosities only, as described by Yan et al. (2019b, Eq. (1)), the real galactic luminosity would be systematically higher due to the AGB stars such that the M_dyn/L_V ratios would be lower than the calculation described here. Therefore, we divide the M_dyn/L_V ratios by the same factor of 1.425 for all galaxies to normalise the M_dyn/L_V values such that the galaxy with a mass of M_dyn = 10^9.5 M_⊙ has M_dyn/L_V = 3 in solar units. The E-MILES stellar population model (Vazdekis et al. 2016) combined with the GalIMF code will be applied to give more accurate M_dyn/L_V values in the near future. The results are shown and compared with the mass-to-light ratios from direct observations in Fig. 10. The agreement is remarkable: the here-computed best-fit models follow the data trend very well, with the implication that the masses of the most massive E galaxies are significantly dominated by faint M dwarfs, and stellar remnants are in agreement with the spectroscopic studies (see discussion in Loubser et al. 2021 and references therein). These massive galaxies have higher M_dyn/L_V values due to the fact that their bulk super-solar metallicity leads to a bulk bottom-heavy gwIMF, as is shown by the red line in Fig 9.

Fig. 9. gwIMFs, defined in Eq. (12), for the best-fit galaxies of our error function model (Eq. (21)) with masses log10(Mdyn/M⊙) = 8.5, 9.0, 9.5, …, 11.5, and 11.9 from the bottom to the top line. The galaxies with log10(Mdyn/M⊙) = 8.5, 9.5, and 11.9 are highlighted in the green, blue, and red, respectively. The galactic SFRs ($\lg \psi \equiv {\log }_{10}{\overline{\psi }}_{\delta t}[M_{\odot }/\mathrm{yr}]$) of these models are given by the legend in units of M⊙/yr. The red gwIMF has an approximated α3 of about 1.2 although the top part of the IGIMF no longer preserves a power-law form.

Fig. 10. Normalised dynamical-mass-to-light ratio as a function of the dynamical mass of the galaxies (the grey circles are observational data given by Dabringhausen & Fellhauer 2016, DF16) at different galactic ages (denoted by the legend). The green, blue, and red points, i.e., galaxies with a final mass of Mdyn = 108.5, 109.5, and 1011.9 M⊙ are calculated with the gwIMF of the same colour given in Fig. 9.

Unlike the gwIMF calculated here using the simple-population approximation and the observed mean stellar metallicity, the real gwIMF evolves through each star formation epoch. At the onset of the formation of all galaxies, the gwIMFs were all similarly bottom-light, with a top-heaviness depending on the SFR and the gwIMF becoming more bottom-heavy as the metallicity of the galaxy increases. For example, Fig. 11 plots the gwIMF for each 10 Myr star formation epoch and also for the integrated gwIMF of all star formation epochs, i.e., a time-integrated gwIMF (TIgwIMF) for a galaxy with ${\log }_{10}({\overline{\psi }}_{\delta t}[M_{\odot }/\mathrm{yr}])=3.4468$, τ_SF = 1.19 Gyr, g_conv = 0.1898, [Z/X](t = 0)= − 6, and a final mass of 10^11.9 M_⊙, corresponding to the high-SFR red gwIMF in Fig. 9. The TIgwIMF for an evolved stellar population with different ages reproduces the gwIMF calculated using the observed mean stellar metallicity well, indicating that it is valid to calculate M_dyn/L_V with the simple-population approximation. Similarly, the TIgwIMFs for galaxies with a mass of 10^8.5 or 10^9.5 M_⊙ are shown in Figs. C.1 and C.2, respectively. Their TIgwIMFs are similar to the canonical IMF, leading to a standard M_dyn/L_V value. We refer the reader to Dabringhausen et al. (2016) and Dabringhausen (2019) for further analysis of the dynamical mass to light ratios of ETGs within the IGIMF framework.

Fig. 11. gwIMF, as defined in Eq. (12), of each 10 Myr star formation epoch (thin lines) and the time-integrated gwIMF for all formation epochs (TIgwIMF) for a galaxy with a final mass of Mdyn = 1011.9 M⊙. This TIgwIMF corresponds to the red model in Fig. 10. The corresponding approximating simple-population IGIMF model is shown in Fig. 9 as a red line. The gwIMF for the first and the last star formation epochs are highlighted by slightly thicker lines. The canonical IMF as given by Kroupa (2001) and the power-law IMF given by Salpeter (1955) are shown via the red dotted and blue dashed lines, respectively.

6.2. Stellar yield tables

The estimated SFT for a galaxy with a given mass is highly dependant on the adopted stellar yield table. However, since the yield is applied to all galaxies uniformly, the general trend (the slope) of the τ_SF–M_dyn relation and the scientific conclusions remain the same. We verified that the effect of applying a different yield table on the resulting τ_SF–M_dyn relation is similar to having a uniform bias of the observed galactic metal abundance. This was demonstrated in our previous work (Yan et al. 2019a, Sect. 5.1).

6.3. Application of closed-box models

The closed-box modelling approach has the advantage that a large range of parameters (e.g., the galaxy mass, formation timescale, gwIMF properties) are computationally reachable, while cosmological hydrodynamical galaxy-formation models with star formation and feedback processes are computationally too costly to allow a comprehensive investigation of parameter space.

The observations of E galaxies have shown that in situ star formation contributes a large fraction of the stellar population, while galaxy mergers are invoked to gradually increase the size and the mass of the most massive E galaxies over a long period of time, building up the outer regions of the galaxies (Renzini 2006; Oser et al. 2010; Navarro-González et al. 2013; McDermid et al. 2015; Ferré-Mateu et al. 2015; Martín-Navarro et al. 2018, 2019b; Zibetti et al. 2020; Lacerna et al. 2020; for a discussion, see Sect. 8.2 in Kroupa et al. 2020b). Since the timescale for most of the stars to form is short, the amount of galactic gas inflow and/or outflow cannot be significant enough to strongly affect the abundances of these stars (Sect. 5.1 in Yan et al. 2019b). For this reason, the central region of giant ellipticals can be reasonably approximated by the closed-box model with neither galactic inflow nor outflow (Matteucci 2012).

On the other hand, the theory of standard-dark-matter-based cosmology predicts that more massive galaxies form from more mergers instead of the in-situ star formation described by the closed-box model (e.g., Moster et al. 2020), in conflict with the observational studies. Reproducing the abundance scaling relations within the hierarchical galaxy formation framework has long been a problem. By tuning the star formation and feedback parameters, it has been shown possible for the simulated galaxies to broadly fulfil the observed mass–metallicity relation for galaxies at different redshifts and also the radial metallicity gradient of resolved galaxies (Hirschmann et al. 2016; Torrey et al. 2018, 2019; Hemler et al. 2021). However, the reproduction for the most massive galaxies is not entirely satisfactory. The observed M_dyn–stellar-metallicity relation has a higher metallicity and a smaller scatter of the metallicity value for the more massive galaxies, which is not seen in the simulations (e.g., the comparison between observed and predicted relation in Arrigoni et al. 2010, Fig. 7, and Barber et al. 2018, top panel of Fig. 11). In addition, the α-enhancement abundances of the massive galaxies are not simultaneously reproduced with the metal-enhancement (Okamoto et al. 2017; De Lucia et al. 2017; Rosito et al. 2019). This turns out to be a tricky problem requiring non-canonical solutions including a non-canonical IMF (Nagashima et al. 2005; Calura & Menci 2009; Arrigoni et al. 2010; Gargiulo et al. 2015; Fontanot et al. 2017; Barber et al. 2018; Gutcke & Springel 2019), non-canonical SNIa incidence (Arrigoni et al. 2010), non-canonical stellar yields (see Sect. 6.4), or differential galaxy winds (Yates et al. 2013). Hence, the hierarchical formation scenario appears to need further tuning and parameter-addition of the baryonic physics parameters.

Even if the merger scenario for more massive galaxies is correct, it does not necessarily disqualify the closed-box galaxy chemical evolution models. The satellite galaxies that merged into and contributed a large mass fraction to the massive galaxies should not have the same property as the low-mass galaxies surviving to date that have undergone processing over a Hubble time. Otherwise, the massive galaxies would not preserve their unique higher metallicity. The merged satellite galaxies are thought to have been more metal-rich than the central galaxies of the same mass (Pasquali et al. 2010) since they have lived in the vicinity of a massive galaxy or in a group of galaxies with a higher total mass. Thus, the satellites are not isolated but subject to the chemical enrichment of the entire gravitationally bound gas cloud and/or galaxy group (Bahé et al. 2017) that later merged into a single massive galaxy where the enriched galactic outflows of each galaxy in the group are recycled back within the group. In this sense, the closed-box approximation can broadly describe the chemical evolution of the massive E galaxies when considering the entire progenitor galaxy group as a whole, while the exact dynamical evolution and assembly history are not a concern here. This description is certainly only an approximation, and not all galaxies can be well represented by the closed-box model. There have been observations of outliers with unexpected abnormal abundance ratios by the closed-box model (e.g., the iron-poor population documented in McDermid et al. 2015 and Jafariyazani et al. 2020). Such galaxies can only be explained by a more complicated formation scenario; for example, through a larger amount of primordial gas inflow or other mechanisms that modify the chemical abundance. Nevertheless, taken at face value, the properties of ETGs are simplest and most concisely understood through the monolithic post-Big Bang gas-cloud collapse theory for their origin (Kroupa et al. 2020b).

6.4. Alternative solutions assuming a canonical IMF

We asked ourselves if it was possible for different gwIMF and/or different stellar yields, different SNIa DTDs, or gas mixing and expulsion physics, or even a combination of these to reproduce the observed chemical abundance ratios of E galaxies introduced in Sect. 1. Moreover, we wanted to know if it was possible for the most massive E galaxies to reach their high mean stellar [Mg/Fe] value with an SFT no shorter than about 1 Gyr, assuming that the gwIMF is invariant and canonical.

We first considered the possibility of a variable κ_Ia as a function of galaxy mass or the SFR when the IMF is canonical and invariant. In order to reproduce the high [Mg/Fe] ratios of massive E galaxies, smaller κ_Ia values would be required, and, therefore, this would contradict the observation. For example, Arrigoni et al. (2010) found that, with the canonical and a slightly top-heavy gwIMF, the [α/Fe] values of massive galaxies require too low an SNIa production efficiency, which is hard to reconcile with the observed SNIa rate of the star-forming galaxies. Also, the expectation that high-density, metal-rich star-forming regions in a massive E galaxy are likely to have a higher fraction of short-period binary stars and thus higher κ_Ia values (Shara & Hurley 2002; Greggio 2005, Sect. 2 and references therein, and Friedmann & Maoz 2018). Furthermore, if the bulk stellar population of E galaxies forms through a monolithic collapse of a post-Big Bang gas cloud, as is suggested by much observational evidence and theoretical modelling (Kroupa et al. 2020b and references therein), then these galaxies would be forming extremely massive star-burst clusters in-line with the observed extragalactic correlation between the SFR and cluster masses (Weidner et al. 2004; Randriamanakoto et al. 2013; Yan et al. 2017). As massive star-burst clusters are SNIa factories (Shara & Hurley 2002) and given the above arguments by Greggio (2005), it appears likely that κ_Ia indeed increases with increasing E galaxy mass. We show in Sect. 5.2.2 that κ_Ia indeed needs to increase for the more massive E galaxies under the framework of the IGIMF theory.

The other solutions include a more efficient loss of the elements produced by SNIa than those produced by type II supernovae. This could help to increase the galactic [Mg/Fe] values, allowing for longer SFT values. However, such a supernova-type-dependent wind would be more prominent in dwarf galaxies that have shallower potential wells (Recchi et al. 2001) than for the massive galaxies. Given the large E galaxy masses, here we apply the closed-box model and neglect any galactic winds. Also, the stellar yields may change if the average stellar rotation speed differs in different galaxies due to, for example, a difference in the stellar metallicity (Romano et al. 2019, 2020). This possibility requires further investigation.

6.5. Comparison with the SFTs estimated with other methods

The downsizing behaviour of galaxy formation is supported independently by dynamical simulations (Wong 2009; Joshi et al. 2021; Eappen et al., in prep.), chemical evolution models (Thomas et al. 2005, 2010), and SPS studies (de La Rosa et al. 2011; McDermid et al. 2015). However, it is not trivial to confirm or falsify the here estimated τ_SF–M_dyn relation (the yellow ridge line in Fig. 7) using different methods since no observation of the SFH is direct and there is a large systematic uncertainty related to all the methods.

For example, the chemical evolution models in this study assume a perfect galaxy with no gas flows, a constant SFR during the star formation epoch, and zero SFR before and after that. These approximations are likely to be valid only for massive E galaxies. The observation of dwarf galaxies with a low total stellar luminosity can be significantly dominated by their most recent star formation from a gas reservoir polluted by the local environment, invalidating our model assumptions. In addition, as mentioned in Sect. 5.2.2, the real fluctuating SFR in dwarf galaxies rather than the assumed smooth SFH would lead to a more top-heavy gwIMF, and therefore to a higher estimated SFT than in our results shown in Fig. 7.

The SFTs suggested by the SPS studies, such as that by McDermid et al. (2015), are subject to the assumed stellar population model, dust attenuation law, and ad hoc preference of a smooth SFH. The assumed gwIMF and stellar abundance ratios of different stellar populations being synthesised are prior assumptions that introduce biases to the resulting SFH. Due to the degenerate nature of inferring the SFH from the integrated galaxy light, the solution is not unique (Young et al. 2014; Aufort et al. 2020). The SPS method is not able to detect strong short-term fluctuations, and the resulting SFH is strongly affected by the prior assumed SFH shape and/or smoothness constraints (McDermid et al. 2015; Iyer & Gawiser 2017; Carnall et al. 2018, 2019; Leja et al. 2019; Iyer et al. 2019; Lower et al. 2020; Tacchella et al. 2021). With higher SFRs extending for shorter times, the total stellar mass formed with the same age and metallicity is preserved, while the gwIMF can be very different if it depends on the star formation intensity, as is suggested by the IGIMF theory (Sect. 2.1). As mentioned in Yan et al. (2019a), with non-unique SPS solutions, there is no guarantee that the preferred solution (imposed by the prior) will reproduce the observed abundance in a self-consistent chemical evolution simulation (see Bellstedt et al. 2020 as a first step to consider the SPS and the metal-evolution of the galaxy consistently).

The SFH measurements would be much more reliable and encouraging if the estimations from all different methods were be consistent with each other. More recent studies of high-redshift galaxies by Saracco et al. (2019) broadly agree with the downsizing relation suggested by McDermid et al. (2015). Independent methods such as measuring the SFH with pixel colour-magnitude diagrams (CMD, Cook et al. 2019) could also be helpful supplementary evidence supporting the decreasing SFHs of massive ETGs. The uncertainty of such a method is large, but better-resolved pixel CMD images have the potential to reveal the SFH fluctuations that are hidden in integrated-light observations. Larger future observational platforms are likely to improve our understanding of the assembly of ETGs.

With the consensus that the SFT of massive ETGs should be longer than about 1 Gyr, we are able to robustly exclude the canonical invariant IMF model in Yan et al. (2019a). We demonstrate in this work that an agreement between the estimate of the SFTs of massive ETGs using the SPS method on one hand and chemical evolution studies on the other can be achieved if the IGIMF theory is applied (Fig. 7). The required higher SNIa production efficiency for massive ETGs (Fig. 6) agrees beautifully with the independent theoretical expectation (Shara & Hurley 2002) and observation (Freundlich & Maoz 2021) for the first time. Conversely, the SFHs of dwarf galaxies, which have much shallower potentials, non-canonical gwIMFs, and dust, is harder to determine. We avoid fine-tuning our model or introducing additional model parameters. At face value, our study suggests a short SFT for the typical early-type dwarf galaxies.

7. Conclusion

In this work and for the first time, we estimate the formation timescale of E galaxies from their observed stellar mean metallicity and [Mg/Fe] values under the framework of an environment-dependent IMF, that is, the IGIMF theory. The IMF changes as a consequence of the changing metallicity and SFR of the galaxy, which in turn affects the chemical enrichment. In Eq. (15), we formulate how the total number of SNIa per unit mass of stars formed, N_Ia, should vary with the IMF shape. The results assuming an invariant overall SNIa realisation parameter, i.e., our fiducial SNIa model, are shown in Fig. 2, accounting for the observed large scatter of the observational estimates. A higher N_Ia is inferred for metal-poor galaxies that have a bottom-light IMF. On the other hand, a lower N_Ia is inferred for low-SFR galaxies that have a top-light gwIMF.

The environment-dependent gwIMF affects the SFT evaluation significantly. The resulting SFT is longer for massive E galaxies than the value estimated by the SPS studies. A natural explanation for this discrepancy is that the real total number of SNIa per stellar mass formed is higher for more massive galaxies than the number estimated by the fiducial SNIa model because they had more intense and denser star-forming activities. A plausible 2.5 times increase of the overall SNIa realisation parameter for the most massive galaxies (the error function model shown in Fig. 5), as suggested by Friedmann & Maoz (2018) and Freundlich & Maoz (2021), can resolve this discrepancy and fits the SFT values obtained from SPS. This is consistent with massive galaxies having formed more massive starburst clusters as a consequence of their higher SFRs, because such clusters are factories for SNIa progenitors and are also relevant for the rapid emergence of super-massive black holes (Shara & Hurley 2002; Kroupa et al. 2020b). Future dynamical studies may be able to further test this prediction with N-body simulations. On the other hand, adopting the invariant canonical IMF in the galaxy chemical evolution model would not make it possible to reproduce these galaxy SFTs and metal abundance constraints, because such a model would require a lower SNIa production efficiency for massive galaxies than for low-mass galaxies.

The IGIMF theory with a more bottom-light gwIMF for the low-SFR galaxy also suggests that the SFT of E galaxies should decrease for lower mass galaxies. Considering that the SFT also decreases for the most massive galaxies, as is suggested by the SPS studies, galaxies with a mass of about 10^9.5 M_⊙ may have the longest SFT (Fig. 8). This result implies that most stars in more massive E galaxies were formed in a shorter post-Big Bang gas-cloud collapse timescale, with roughly the same fraction (1/3) of gas converted to stars and a higher mass-to-light ratio, while lower-mass galaxies tend to lose their gas supply and quench star formation due to feedback over a shorter timescale. Examples of this may be the ultra-diffuse galaxy Dragonfly 44 (Haghi et al. 2019) and ultra-faint dwarf galaxy Boötes I (Yan et al. 2020).

We conclude that the IGIMF theory along with an increased production efficiency of SNIa for more massive E galaxies (described by Eqs. (15), (20), and (21)) is able to explain the observed stellar populations of E galaxies successfully, thereby also self-consistently accounting for the large dynamical mass-to-light ratios (Fig. 10) and the bottom-heavy stellar mass function in massive E galaxies. A noteworthy aspect of the theory of elliptical galaxy formation and evolution arrived at here is that many observed properties come together naturally and are in fact straightforwardly explained: the higher incidence of SNIa events in galaxies that formed over the downsizing time with high SFRs, which is the downsizing-implied systematic galaxy-wide IMF evolution accounting for the metallicity and α-element enhancement and the implied rapid formation of supermassive black holes (Kroupa et al. 2020b), all are consequences of the application of the IGIMF theory to the monolithic collapse of post-Big Bang gas clouds (Eappen et al., in prep.).

Acknowledgments

Z.Y. acknowledges financial support from the China Scholarship Council (CSC, file number 201708080069). T.J. acknowledges support through the European Space Agency fellowship programme. P.K. acknowledges support from the Grant Agency of the Czech Republic under grant number 20-21855S. The development of our chemical evolution model applied in this work benefited from the International Space Science Institute (ISSI/ISSI-BJ) in Bern and Beijing, thanks to the funding of the team ‘Chemical abundances in the ISM: the litmus test of stellar IMF variations in galaxies across cosmic time’ (Donatella Romano and Zhi-Yu Zhang).

References

Appendix A: SNIa production efficiency calculated with different IGIMF formulations and an SNIa progenitor mass range

The SNIa production efficiency calculated with a different SNIa progenitor mass range from 2 to 8 M_⊙ is shown in Fig. A.1. The efficiency is higher for the low-mass, metal-poor galaxies compared to Fig. 2. Results are similar if the considered mass ranges for the primary and the secondary star are different; for example, from 3 to 8 M_⊙ and 1.5 to 8 M_⊙, respectively (Fig. A.2).

Fig. A.1. Ten Gyr time-integrated number of SNIa per unit stellar mass formed. Same as Fig. 2, but with an SNIa progenitor mass range from 2 to 8 M⊙ instead of 3 to 8 M⊙.

Fig. A.2. Ten Gyr time-integrated number of SNIa per unit stellar mass formed. Same as Fig. 2, but with an SNIa progenitor mass range from 3 to 8 M⊙ for the primary star and 1.5 to 8 M⊙ for the secondary star instead of 3 to 8 M⊙ for both stars.

Figures A.3 to A.5 give the SNIa production efficiency calculated for the IGIMF1, IGIMF2, and IGIMF3 formulations as summarised in Jeřábková et al. (2018, Table 3). Models with different metallicities are overlapping in Fig. A.3 because the IGIMF1 model is independent of metallicity. The IGIMF2 model has an invariant IMF for low-mass stars. The IGIMF3 model has an IMF dependence on [Z] instead of Z in Eq. 4 for low-mass stars. The results shown in this section following Eq. 15 can be used to exclude theories of how the IMF varies.

Fig. A.3. Ten Gyr time-integrated number of SNIa per unit stellar mass formed. Same as Fig. 2, but adopting the IGIMF1 model of Jeřábková et al. (2018) instead of Eqs. 4 and 5.

Fig. A.4. Ten Gyr time-integrated number of SNIa per unit stellar mass formed. Same as Fig. 2, but adopting the IGIMF2 model of Jeřábková et al. (2018) instead of Eqs. 4 and 5.

Fig. A.5. Ten Gyr time-integrated number of SNIa per unit stellar mass formed. Same as Fig. 2, but adopting the IGIMF3 model of Jeřábková et al. (2018) instead of Eqs. 4 and 5.

Appendix B: Power-law function model

Fig. B.1 shows the best-fit τ_SF as a function of M_dyn for the power-law function model (see Fig. 5) by the yellow ridge line.

Fig. B.1. τSF as a function of Mdyn for the power-law function model. Same as Fig. 4, but the results shown by the colour map adopt the power-law κIa–${\overline{\psi }}_{\delta t}$ relation shown in Fig. 5.

Appendix C: Evolution of the gwIMFs

Fig. C.1 shows the gwIMF for each 10 Myr star formation epoch and the TIgwIMF for a galaxy with ${\log }_{10}({\overline{\psi }}_{\delta t}[M_{\odot }{\mathrm{yr}}^{-1}])=-0.2693$, τ_SF = 1.11 Gyr, g_conv = 0.3204 and a final mass of 10^8.502 M_⊙, corresponding to the lowest-SFR green gwIMF in Fig. 9. Fig. C.2 shows the gwIMF for each 10 Myr star formation epoch and the TIgwIMF for a galaxy with ${\log }_{10}({\overline{\psi }}_{\delta t}[M_{\odot }{\mathrm{yr}}^{-1}])=0.527$, τ_SF = 1.88 Gyr, g_conv = 0.2796, and a final mass of 10^9.4837 M_⊙, corresponding to the blue gwIMF in Fig. 9.

Fig. C.1. Same as Fig. 11, but for a galaxy with a final mass of Mdyn ≈ 108.5 M⊙. This TIgwIMF is well represented by the green model in Fig. 9.

Fig. C.2. Same as Fig. 11, but for a galaxy with a final mass of Mdyn ≈ 109.5 M⊙. This TIgwIMF is well approximated by the blue model in Fig. 9.

All Tables

All Figures

Fig. 1. Total number of SNIa explosions for a SSP after 10 Gyr of the birth of the SSP per unit stellar mass of the SSP predicted by Eqs. (13)–(16), assuming constant parameter Bbin and CIa and an IMF defined in Eq. (1), with α1 = 1.3, α2 = 2.3, as a function of α3. A smaller α3 indicates a flatter IMF/gwIMF for the massive stars, that is, a top-heavy IMF/gwIMF. Very large α3 mean that essentially no stars above 1 M⊙ are formed. The function is calibrated by the vertical and horizontal dashed lines indicating the canonical IMF/gwIMF slope, α3 = 2.3, and the value nIa(t = 10 Gyr, ξcanonical) = 0.0022/M⊙ suggested by Maoz & Mannucci (2012) for the local Universe, respectively.

In the text

Fig. 2. Ten Gyr time-integrated number of SNIa per unit stellar mass formed for the fiducial SNIa model calculated according to Eqs. (13)–(16). The fiducial SNIa model assumes a constant overall SNIa realisation parameter and the gwIMF as given by the IGIMF theory (Eq. (12)) as a function of the galaxy-wide SFR, ${\overline{\psi }}_{\delta t}$, and metallicity, [Z/X]. The black line is the relation for [Z/X] = 0. Other lines with different colours represent different values of [Z/X], as indicated by the white stripes on the colour map on the right: [Z/X] = 0.3, 0.2, 0.1, −0.1, −0.2, −0.3, −0.4, −0.5, −1, −2, −4, and −6. The black horizontal dotted line represents the canonical nIa(t = 10 Gyr, ξcanonical) value of 0.0022/M⊙ (Maoz & Mannucci 2012). The green and red horizontal dashed lines indicate observational constraints on the nIa(t = 10 Gyr) for SNIa surveys up to a certain redshift and in galaxy clusters, respectively. The shaded regions represent the uncertainty ranges of the horizontal dashed lines. References are given in Table 1.

In the text

Fig. 3. [Z/X], [Mg/Fe], and age of the E galaxies as a function of their dynamical masses. The grey circles are observed galaxies from Arrigoni et al. (2010), with the typical uncertainty shown at the lower right corner of each panel. The solid and dashed lines are, respectively, the mean and standard deviations of the observational constraints. The lines in the upper panel are given by Arrigoni et al. (2010) and smoothed by Yan et al. (2019a). The mean relation for the middle panel follows Thomas et al. (2005), while the standard deviation of the relation follows Liu et al. (2016). The points and triangles are the resulting Mdyn, [Z/X], and [Mg/Fe] values of our fiducial galaxy models adopting the IGIMF theory and a constant overall SNIa realisation parameter defined in Sect. 2.2. The triangles are models with ${\log }_{10}({\overline{\psi }}_{\delta t})=0,$ while the points are models with other SFRs. The red and the blue dots are models with a SFT τSF ≥ 1 Gyr and τSF <  1 Gyr, respectively. The darker enhanced colours represent models with a higher gconv (defined in Sect. 4).

In the text

Fig. 4. τSF as a function of Mdyn. The colour bar indicates the goodness of the fit, p(gconv, τSF) (Eq. (17)), for our fiducial SNIa model. For a given mass and SFT, τSF, any star-to-gas-mass fraction, gconv, is allowed. The most likely SFT for a given galaxy mass is indicated by the yellow ridge line. The black dotted curve (Eq. (18)) is a power-law fit of the yellow ridge line. The red dashed curve and the red dash-dotted curve (Eq. (19)) represent the SFT constrained by SPS studies by de La Rosa et al. (2011) and McDermid et al. (2015), respectively. As is introduced in Sect. 1, The SFT suggested by chemical evolution studies assuming the invariant canonical gwIMF performed by Pipino & Matteucci (2004) and Thomas et al. (2005) are shown by the blue triangles and blue dotted curve, respectively. The lower and upper yellow dotted curves correspond, respectively, to polynomial regressions of the yellow ridge lines in Figs. 5 and 7 of Yan et al. (2019a), derived using invariant canonical gwIMF. The mass ranges of the coloured curves are limited by the mass ranges of the galaxy data sets they are based on.

In the text

Fig. 5. SNIa realisation re-normalisation parameter, κIa defined in Eq. (20), as a function of the galaxy-wide SFR, ${\overline{\psi }}_{\delta t}$. Three different models are tested, and the corresponding best-fit SFT–galaxy-mass relations are shown in Figs. 4, 7, and B.1 for the fiducial model (dotted line), error function model (solid line), and power-law function model (dashed line), respectively. The error function model, formulated in Eq. (21), best reproduces the observed ETGs.

In the text

	Fig. 6. Ten Gyr time-integrated number of SNIa per unit of stellar mass formed for the error function model. Same as Fig. 2, but here the calculation assumes the variable κIa–${\overline{\psi }}_{\delta t}$ relation defined by Eq. (21) and is shown by the solid curve in Fig. 5 instead of a constant overall SNIa realisation parameter (dotted line in Fig. 5).
In the text

Fig. 7. τSF as a function of Mdyn for the error function model. Same as Fig. 4, but the results shown by the colour map adopt the variable ${\kappa }_{\mathrm{Ia}}({\overline{\psi }}_{\delta t})$ defined by Eq. (21) and shown by the solid curve in Fig. 5. The best-fit solutions (the yellow ridge line) for different Mdyn can be described by two relations, i.e., the black dotted curve (Eq. (18)) for Mdyn <  109.5 M⊙, and the red dash-dotted curve (Eq. (19)) for Mdyn >  109.8 M⊙.

In the text

Fig. 8. SFT, SFR, gas-conversion parameter, and dynamical-mass-to-stellar mass ratio of our best-fit model, from the top to the bottom panel, for the error function model. The solid black line in the top panel is the same as the yellow ridge line of Fig. 7. The black dotted and the red dash-dotted curves in the top panel follow the same relation as in Fig. 7, described by Eqs. (18) and (19), respectively. The thin horizontal lines indicate gconv values applied in our galaxy model (see Sect. 4), where gconv ≈ 1/3 is about the star formation efficiency found for embedded star clusters (Sect. 2.1).

In the text

Fig. 9. gwIMFs, defined in Eq. (12), for the best-fit galaxies of our error function model (Eq. (21)) with masses log10(Mdyn/M⊙) = 8.5, 9.0, 9.5, …, 11.5, and 11.9 from the bottom to the top line. The galaxies with log10(Mdyn/M⊙) = 8.5, 9.5, and 11.9 are highlighted in the green, blue, and red, respectively. The galactic SFRs ($\lg \psi \equiv {\log }_{10}{\overline{\psi }}_{\delta t}[M_{\odot }/\mathrm{yr}]$) of these models are given by the legend in units of M⊙/yr. The red gwIMF has an approximated α3 of about 1.2 although the top part of the IGIMF no longer preserves a power-law form.

In the text

Fig. 10. Normalised dynamical-mass-to-light ratio as a function of the dynamical mass of the galaxies (the grey circles are observational data given by Dabringhausen & Fellhauer 2016, DF16) at different galactic ages (denoted by the legend). The green, blue, and red points, i.e., galaxies with a final mass of Mdyn = 108.5, 109.5, and 1011.9 M⊙ are calculated with the gwIMF of the same colour given in Fig. 9.

In the text

Fig. 11. gwIMF, as defined in Eq. (12), of each 10 Myr star formation epoch (thin lines) and the time-integrated gwIMF for all formation epochs (TIgwIMF) for a galaxy with a final mass of Mdyn = 1011.9 M⊙. This TIgwIMF corresponds to the red model in Fig. 10. The corresponding approximating simple-population IGIMF model is shown in Fig. 9 as a red line. The gwIMF for the first and the last star formation epochs are highlighted by slightly thicker lines. The canonical IMF as given by Kroupa (2001) and the power-law IMF given by Salpeter (1955) are shown via the red dotted and blue dashed lines, respectively.

In the text

	Fig. A.1. Ten Gyr time-integrated number of SNIa per unit stellar mass formed. Same as Fig. 2, but with an SNIa progenitor mass range from 2 to 8 M⊙ instead of 3 to 8 M⊙.
In the text

	Fig. A.2. Ten Gyr time-integrated number of SNIa per unit stellar mass formed. Same as Fig. 2, but with an SNIa progenitor mass range from 3 to 8 M⊙ for the primary star and 1.5 to 8 M⊙ for the secondary star instead of 3 to 8 M⊙ for both stars.
In the text

	Fig. A.3. Ten Gyr time-integrated number of SNIa per unit stellar mass formed. Same as Fig. 2, but adopting the IGIMF1 model of Jeřábková et al. (2018) instead of Eqs. 4 and 5.
In the text

	Fig. A.4. Ten Gyr time-integrated number of SNIa per unit stellar mass formed. Same as Fig. 2, but adopting the IGIMF2 model of Jeřábková et al. (2018) instead of Eqs. 4 and 5.
In the text

	Fig. A.5. Ten Gyr time-integrated number of SNIa per unit stellar mass formed. Same as Fig. 2, but adopting the IGIMF3 model of Jeřábková et al. (2018) instead of Eqs. 4 and 5.
In the text

	Fig. B.1. τSF as a function of Mdyn for the power-law function model. Same as Fig. 4, but the results shown by the colour map adopt the power-law κIa–${\overline{\psi }}_{\delta t}$ relation shown in Fig. 5.
In the text

	Fig. C.1. Same as Fig. 11, but for a galaxy with a final mass of Mdyn ≈ 108.5 M⊙. This TIgwIMF is well represented by the green model in Fig. 9.
In the text

	Fig. C.2. Same as Fig. 11, but for a galaxy with a final mass of Mdyn ≈ 109.5 M⊙. This TIgwIMF is well approximated by the blue model in Fig. 9.
In the text