Variation of the stellar initial mass function in semi-analytical models III: testing the cosmic ray regulated integrated galaxy-wide initial mass function

In our previous work, we derive the CR-IGIMF: a new scenario for a variable stellar initial mass function (IMF), which combines numerical results on the role played by cosmic rays in setting the thermal state of star forming gas, with the analytical approach of the integrated galaxy-wide IMF. In this work, we study the implications of this scenario for the properties of local Early-Type galaxies (ETG), as inferred from dynamical, photometric and spectroscopic studies. We implement a library of CR-IGIMF shapes in the framework of the Galaxy Evolution and Assembly (GAEA) model. Our realization includes a derivation of synthetic spectral energy distribution for each model galaxy, allowing a direct derivation of the mass fraction in the mean IMF of low-mass stars (i.e. the dwarf-to-giant ratio - $\rm f_{dg}$) and a comparison with IMF sensitive spectral features. The predictions of the GAEA model implementing the CR-IGIMF confirm our previous findings: it correctly reproduces both the observed excess of z$\sim$0 dynamical mass (mass-to-light ratios) with respect to spectroscopic (photometric) estimates assuming a universal, MW-like, IMF, and the observed increase of [$\alpha$/Fe] ratios with stellar mass in spheroidal galaxies. Moreover, this realization reproduces the increasing trends of $\rm f_{dg}$, and IMF-sensitive line-strengths with velocity dispersion, although the predicted relations are significantly shallower than the observed ones. Our results show that the CR-IGIMF is a promising scenario that reproduces at the same time dynamical, photometric and spectroscopic indications of a varying IMF in local ETGs. The shallow relations found for spectral indices suggest that either a stronger variability as a function of galaxy properties or additional dependences (e.g. as a function of star forming gas metallicity) might be required to match the strength of the observed trends.


Introduction
The stellar initial mass function (IMF) represents a fundamental ingredient for theoretical models of galaxy formation and evolution.For each given star formation episode, this statistical function defines the number of stars formed per stellar mass bin and therefore it determines the number of SNe and the fraction of baryonic mass locked in long-living stars.Direct observations of star-forming regions in the Milky Way (MW) show a remarkable consistency in the derived shape of the IMF.Different authors proposed different functional shapes for the IMF in our galaxy (Salpeter 1955;Kroupa 2001;Chabrier 2003), that mainly differ in the very uncertainty regime of brown dwarfs, i.e. at the low-mass end.Given the fact that the MW is the only galaxy where we can currently measure the IMF using direct stellar countings, the near invariance of the results prompted the idea of its universality.However, from a theoretical perspective, the physical description of gravitational collapse and fragmentation of molecular clouds (MCs) is still debated in the literature (see e.g.Krumholz 2014): several models (Klessen et al. 2005;Weidner & Kroupa 2005;Hennebelle & Chabrier 2008;Papadopoulos et al. 2011;Narayanan & Davé 2013, among others) explored the range of possible IMF shapes depending on the physical properties of the interstellar medium (ISM).
A&A proofs: manuscript no.cr_igimf.finA first indication that the IMF can indeed vary in environments different from those of the MW disc, come from Klessen et al. (2007), who show that the densest regions of the MW bulge (where stellar counts are still accessible) are better described by a different IMF (i.e. with a larger fraction of massive stars) with respect to the star forming regions in the disc.Nonetheless, the assumption of a universal IMF 1 is still used in almost all extragalactic studies.The main limitation for the study of the IMF shape over a variety of environments is our inability to resolve individual star clusters in external galaxies.Therefore, the only information currently accessible comes from the integrated light from composite stellar populations, binding our ability to detect signatures of deviations from a universal IMF to indirect evidences.
The number of evidences in favour of a varying IMF have steadily increased in recent times.Using galaxies from the GAMA 2 survey, Gunawardhana et al. (2011) show that their optical colours and H α line strengths suggest a flatter high-mass slope of the IMF at increasing SFR.Dynamical studies in a sample of early-type galaxies (ETGs) in the ATLAS 3D survey have revealed a systematic excess of dynamical mass-tolight ratios with respect to the expectations for a universal, MW-like, IMF based on observed photometry.This mass excess, α, increases with increasing galaxy velocity dispersion σ (Cappellari et al. 2012).A number of spectral features in the spectra of ETGs (such as Na I , TiO1 and TiO2), sensitive to temperature and/or surface gravity, have been proposed as tracers of the ratio between low-mass and giant stars, and therefore of the shape of the IMF they are formed from.Several groups (see Conroy & van Dokkum 2012;Ferreras et al. 2013; La Barbera et al. 2013;Spiniello et al. 2014, among the others) have adopted this approach to conclude that the IMF slope at the low-mass end becomes steeper at increasing σ and/or stellar mass M ⋆ .
The main problem in interpreting this wealth of information in a consistent picture lies in our lack of physical understanding on how the IMF shape should vary as a function of the physical conditions of star forming regions.IMFs with more low-mass or high-mass stars than the MW-like IMF are defined "bottomheavy" (or "top-light") and "top-heavy" (or "bottom-light"), respectively.Spectroscopic studies seem to prefer a "bottomheavy" solution for more massive galaxies (i.e. an IMF slope at the low-mass end steeper than the MW-like IMF).However, these studies are based only on the light of evolved stars, and therefore this approach cannot provide any constraint on the high-mass end of the IMF.In principle, the dynamical technique is sensitive to the overall shape of the IMF, however, it cannot distinguish between a "bottom-heavy" or a "top-heavy" scenario (Tortora et al. 2013).In fact, the observed α excess of the massto-light ratio can be interpreted either as a larger fraction of lowmass stars (giving rise to direct mass excess) or as a larger fraction of massive short-lived stars (giving rise to a decrease in the total light at late times).
The situation is further complicated by the detection of radial IMF gradients in ETGs (Martín-Navarro et al. 2015;La Barbera et al. 2017;van Dokkum et al. 2017;Sarzi et al. 2018;La Barbera et al. 2021).The results suggest that the lowmass end of the IMF becomes steeper towards the inner galaxy regions, while being comparable with a MW-like IMF in the 1 In the following we will consider as a reference universal IMF a Kroupa/Chabrier IMF, i.e. a multi-slope broken power-law, with a steeper high-mass end and a break at around 1 M ⊙ . 2 Galaxy And Mass Assembly outer regions.However, the same galaxies are more α enriched in the inner regions than in the outskirts, and these [α/Fe] gradients may be better explained by a larger number of Type II SNe in the center (thus suggesting an IMF becoming more "top-heavy", towards the central regions).A possible solutions to these tensions in the interpretation of the available data may be a time-varying IMF, where the two ends are allowed to vary independently over different time-scales (see e.g.Weidner et al. 2013;Ferreras et al. 2015).The fact that the two slopes can evolve in different directions at the same time makes it difficult to recast the differences with respect to the MW-like IMF into the standard definition of "bottom-heavy" or "top-heavy" IMF.In the following, we use this nomenclature only when discussing predictions related to one of the ends.
Galaxy formation models that self-consistently follow the interplay between IMF and galaxy properties represents a key tool to assess the complexities in the interpretation of the photometric/spectroscopic data.In our previous work, on the variable IMF hypothesis, we tested two different models, namely the integrated galaxy-wide IMF (IGIMF) by Weidner & Kroupa (2005) and the Cosmic Ray (CR) regulated IMF from Papadopoulos et al. (2011, PP11 hereafter).In brief, the former approach estimates the galaxy-wide IMF starting from a limited number of physically and observationally motivated axioms that can be recast as a function of the galaxy SFR.Instead, PP11 provide numerical estimates for the characteristic Jeans mass of young stars (M ⋆ J ) forming in a molecular cloud (MC) embedded in a CR energy density field of a given intensity.As energetic CRs can penetrate deeper into the core of these MCs, they can change the chemical and thermal state of the star forming gas, thus affecting the emerging IMF.In this case, the dependence can be recast as a function of the SFR surface density3 .In Fontanot et al. (2017a, F17) and Fontanot et al. (2018a, F18a) we present theoretical models implementing the IGIMF and PP11 approaches, respectively, and discuss their implications for galaxy formation and evolution.These models are able to reproduce both the [α/Fe] versus stellar mass relation and the α excess of model ETGs at z∼0.Nonetheless, the typical IMF variations in both these models mainly affect only the slope of the high-mass end, while the low-mass end of the IMF is constant by construction.Therefore, these models are not able to reproduce the variation of the low-mass end slope derived from the spectroscopic studies of nearby ETGs.
In order to overcome this limitation, Fontanot et al. (2018b, F18b hereafter) proposed a new IMF derivation (dubbed CR-IGIMF), based on a combination of the two previous approaches.The main feature of the CR-IGIMF is the predicted evolution of both the high-mass and low-mass ends of the IMF, as a function of the physical properties of model galaxies.In F18b, we presented a preliminary analysis of the impact of this scenario on the spectral properties of galaxies, based on idealized star formation histories, that we used to estimate the mean IMF shapes of galaxies at z∼0 as a function of their stellar mass.We concluded that the CR-IGIMF qualitatively reproduces the observed trends.The aim of this work is to explore in a more quantitative way the impact of this new IMF model on galaxy evolution, and in particular to explore the constraints on dynamical and spectroscopic properties of ETGs at z ∼ 0. To this end, we implement the CR-IGIMF scenario in the framework of the semi-analytic model for GAlaxy Evolution and Assembly (gaea), and employ the synthetic stellar population models of Vazdekis et al. (2012), to compute synthetic spectra of model galaxies.These implementations provide a new, self-consistently calibrated, tool to fully describe the IMF evolution along the star formation and assembly history of model galaxies.
This paper is organized as follows.In Section 2 we layout the theoretical framework of the CR-IGIMF derivation, as presented in F18b.In Section 3, we introduce the implementation of the CR-IGIMF in the gaea model.We present our main results and discuss them in the context of galaxy formation and evolution in Section 4. Finally, we outline our main conclusions in Section 5.

CR-IGIMF
In this paper, we take advantage of the recent derivation of the IGIMF ϕ IGIMF as proposed in Fontanot et al. (2018b), which accounts for the effect of the Cosmic Rays on the thermal and chemical properties of star-forming molecular clouds.This approach is similar to the one proposed by Weidner & Kroupa (2005, see also Dib 2023) and it is based on the ansatz of integrating the individual IMFs associated with each MC (ϕ ⋆ (m)), weighted by the MC mass function ϕ CL (M cl ): Eq. 1 critically depends both on the maximum value of the mass of a star cluster to form (M max cl ) and on the largest stellar mass forming in a given cluster (m max ).Using available evidence from the Taurus-Auriga complex (Kroupa & Bouvier 2003), we set M min cl = 5M ⊙ as the minimum mass for star clusters.Additional ansatzs are then used to define the key quantities entering Eq. 1 as a function of the assumed SFR.(2) Eq. 2 describes the scaling of M max cl with SFR, and it has been derived by Weidner et al. (2004) Eq. 3 is used to predict the maximum stellar mass (m max ⋆ ) to form in a cluster of mass M cl .Pflamm-Altenburg et al. (2007) derive it assuming that each stellar cluster contains exactly one m max star.
(4) Eq. 4 represents the functional shape for the star cluster mass function.
Eq. 5 provides the allowed values for the β parameter in Eq. 4. Lada & Lada (2003) suggest β = 2 based on observational results in the local Universe.The flattening of β at high-SFRs is suggested by Gunawardhana et al. (2011) in order to reproduce the stellar population of galaxies in the GAMA survey.
α 3 = 2.35 ρ cl < 9.5 × 104 M ⊙ /pc 3 1.86 − 0.43 log( ρ cl 10 4 M ⊙ /pc 3 ) ρ cl ≥ 9.5 × 10 4 M ⊙ /pc 3 (6) Eq. 6 quantifies the dependence of the high-mass end (α 3 ) on the cluster core density (ρ cl ) as proposed in Marks et al. (2012).This deviation from a Salpeter high-mass end slope has been reported in the literature (see Kroupa et al. 2013 for a review of the subject).It is important to keep in mind that in Eq. 6 we choose to focus on the dependence of α 3 on ρ cl and neglect any possible dependence of α 3 on metallicity.Our choice is functional to provide a framework for the IGIMF derivation based on SFR and SFR density.Jeřábková et al. (2018) explore an alternative derivation of IGIMF considering metallicity and SFR as the main drivers of its variation, finding similar conclusions as in F18b (see discussion below).Finally, the system of equations is closed by considering the link between the cluster core density and cluster mass as given by Eq. 7 (Marks & Kroupa 2012): In the original Weidner & Kroupa (2005) framework, the IGIMF associated with a given SFR event is numerically derived by assuming a universal ϕ ⋆ (m) in each individual star-forming MC.Following our previous works, we assume that the IMF shape in each MC is well described by the Kroupa (2001) multislope functional form: In Eq. 8 we already introduce the 3-slope formalism used in F18b, by considering a variable inner break at m br ; we then fixed the remaining parameters at m low = 0.1 M ⊙ , m 1 = 1.0 M ⊙ , α 1 = 1.3 and α 2 = 2.35.The variable break m br represents the novelty introduced by F18b, as it relaxes the Weidner & Kroupa (2005) ansatz of a universal IMF in individual MCs 4 .Instead, we assume that stars in a given MCs form from an IMF well described by Eq. 8, where the position of the break at m br depends on the characteristic Jeans mass of young stars (M ⋆ J ), that in turn depends on the MC density and on the environment where it lives, and in particular on the CR density field.In detail, for a given CR energy density (U CR ) and core density given by Eq. 7, M ⋆ J is derived using the numerical solutions to the chemical and thermal equations regulating the CR-dominated ISM (i.e.Fig. 4 in P11): Table 1.Analytical fits to the CR-IGIMF.The quantities f ⋆ dg , (1-f ret ) ⋆ and f ⋆ SN2 represent, respectively, the value of the dwarf-to-giant ratio, of the mass fraction locked in low-mass stars and of the Type II SNe fraction, relative to the same quantities computed for a MW-like IMF with a Kroupa (2001) multi-slope functional form.

SFR
The numerical solutions we consider are characterized by an increase of M ⋆ J with increasing U CR , due to the higher CR heating rate associated with the increased energy density.At fixed U CR , PP11 show also that M ⋆ J decreases at increasing ρ cl .This trend may appear counter intuitive, as high-mass MCs are expected to be hotter and the Jeans mass increases with temperature.However, the numerical results suggest that, when the contribution of CR energy density is taken into account, the dependence of the Jeans mass on density is dominant with respect to temperature.

Library
Consistently with our previous work, in the following we employ a CR-IGIMF binned library (see Tab. 1).The use of a predefined, limited, set of IMFs allows us to retain the computational efficiency of the semi-analytic approach without losing precision.We consider 7 equally spaced bins in log(SFR) in the range [-3,3] and 6 U CR values relative to the corresponding MW CR density (U MW , estimated from the MW SFR density).While the 6 U CR /U MW are limited by the original computations in P11, we employ 7 bins in log(SFR) after having tested the numerical convergence of broad band magnitudes predicted by the model (see Sec. 3.2).We numerically integrate Eq. 1 in each of the 42 total bins combinations, and we fit these numerical results with a multi-component power law, consistent with an extension of the Kroupa IMF formalism.We start from Eq. 8 and we then add slopes and break points until we get the required convergence of the fit (i.e. using an additional slope and break point results in a non-significant increase of the reduced χ 2 probability).Our procedure indicates that a multi-component power law, with 5 slopes (α I , α II , α III , α IV , α V , going from the faint end to the bright end) and 4 break points (m I , m II , m III , m IV , going from low mass to high mass) is required to fully describe the IMF shapes predicted by the CR-IGIMF.We collect the best fit values, describing the multi-component power-laws, in Table 1: some of the most relevant CR-IGIMF shapes have been shown in Fig. 1 of F18b.
The possibility of an intrinsic IMF variation of individual MCs has a fundamental impact on the shape of the emerging galaxy-wide IMF predicted by the IGIMF framework.
-For the typical conditions of a MW-like star forming region (U CR ∼ U MW and SFR∼1M ⊙ yr −1 ), the CR-IGIMF predicts a MW-like shape for the IMF, in agreement with local observations.
-At fixed SFR, the low-mass end of the CR-IGIMF steepens as U CR decreases; however, it becomes indistinguishable from a Kroupa-like α 1 slope for U CR larger than ten times the typical MW value.Indeed, the characteristic Jeans mass of young stars decreases as U CR decreases (Fig. 4 in PP11), thus leading to smaller m br (Eq.9).-At fixed U CR , the high-mass end of the CR-IGIMF follows the typical IGIMF expectation, i.e. it becomes shallower at increasing SFR.This is driven by Eq. 6, through the density of the cluster (following Eq. 7).-At U CR /U MW < 10, the low-mass end of the CR-IGIMF steepens as SFR increases.This is due to the combined effect of the flattening of β at high-SFR (Eq.5), with the decreasing of the characteristic Jeans mass of young stars in denser environments (Fig. 4 in PP11).
-If U CR /U MW > 100, the low-mass end slope is again indistinguishable from an α 1 slope.Indeed, in this U CR range, m br is always larger than 1 M ⊙ at all ρ cl scales.
The combination of the last two items implies that, at least for U CR smaller than 10 times the MW environment, a change in the SFR affects both the high-mass and low-mass end slope at the same time.In particular, an increase of the SFR above the typical MW-like values flattens the high-mass end slope and steepens the low-mass end slope of the CR-IGIMF at the same time, leading to a change in the concavity of the IMF, for SFR> 100 M ⊙ yr −1 .
Assuming a variable IMF affects different observable properties of galaxies, as well as their evolution.In order to better quantify this effect we report in Table 1 the values of some key simple stellar population (SSP) properties5 such as the dwarf-togiant ratio (f dg ), the mass locked into long-lived stars (defined as the complementary of the returned fraction f ret ) and the Type II SNe fraction (f SNII ) per unit solar mass formed.In Table 1, we report the ratio of these quantities to the corresponding values for a MW-like IMF with a Kroupa (2001) functional form (i.e.Eq. 8).The fraction of Type II SNe, f SNII , shows the largest spread both in the entire library and at fixed U CR , ranging from more than twice to roughly one third of the expected value for a Kroupa IMF.This spread is driven by the predicted evolution of the high mass end slope.The mass locked in long-lived stars and stellar remnants decreases almost monotonically at increasing SFR (due to the α 3 evolution) and U CR , reflecting mostly the increased number of high-mass stars predicted by the CR-IGIMF.Finally, f dg is systematically lower than the Kroupa IMF expectation for U CR > ∼ 100U MW (due to m br being always larger than 1 M ⊙ at all ρ cl scales), while showing a wider range of values for lower densities.It is worth stressing that, by construction, the maximum and minimum values for the low-mass end slope cannot exceed the assumed α 1 and α 2 , and the latter case implies a f dg value 1.18 larger that corresponding to the Kroupa IMF (see also Sec. 4.3).These values have been set by local observations of individual clouds and are consistent with theoretical expectations (see e.g Hennebelle & Chabrier 2008).It is important to keep in mind that it is unclear if the IMFs of individual MCs should respect these constraints in all physical environments, especially in case strong CR external fields are present (like the typical condition in starburst galaxies, Zhang et al. 2018).
As we already mentioned, Jeřábková et al. (2018, see also Yan et al. 2021 for a more recent extension of the formalism) propose an alternative approach for an expansion of the IGIMF framework, which considers a secondary dependence on metallicity (alongside with SFR).Their results are qualitatively consistent with the scenario considered in this work, the main differences lying at the low-mass end.In particular, in the CR-IGIMF scenario the low-mass end slopes steepen at decreasing SFR densities, and it equals the Kroupa IMF value (i.e.α 1 ) at Σ SFR > ∼ 100 Σ MW .On the other hand, in Jeřábková et al. (2018) the galaxywide IMF (their model IGIMF3) can either be bottom-light or bottom-heavy if for sub-Solar or super-Solar metallicities, respectively.Both approaches associate a shallower high-mass end slope (corresponding to a larger number high-mass stars with respect to the canonical α 2 slope) to S FR > 1M ⊙ yr −1 events.Finally, it is worth stressing that none of these IGIMF derivations predicts, by construction, low-mass end slopes steeper than α 2 .

Semi-analytic model
In this paper, we discuss the implications of the CR-IGIMF in a cosmological framework, by implementing it in the state-of-theart gaea semi-analytic model (SAM).SAMs represent a modelling technique aimed at describing the redshift evolution of galaxy populations across cosmic times.The physical processes responsible for the energy and mass exchanges among the different baryonic components are treated by means of a system of differential equations.Specifically, each process is described using analytical prescriptions that could be either empirically, numerically or theoretically motivated.Several of these prescrip- Fig. 1.Calibration set used for the CR-IGIMF run.Grey points represent the galaxy luminosity functions in different wavebands and at different redshifts including the SDSS g, r and i-band, K-and V-band (the same compilation of observational estimates used in F17 and F18a, see these papers for detailed references).In all panels, the blue line refers to the reference model predictions.
tions require the calibration of a fixed number of free parameters against a selected set of observational constraints.Observational data outside the calibration set then provide genuine predictions.When coupled with a statistical realization of the evolution of Large Scale Structure as traced by Dark Matter Haloes (i.e. a merger tree), SAMs represent a flexible tool to study galaxy evolution in a cosmological volume: they require only a fraction of the computational time usually associated with hydrosimulations, while providing comparable representations of the statistical properties of large galaxy samples.Moreover, the limited computational request allows SAMs to sample their associated multi-dimensional parameter space very efficiently.The main disadvantage with respect to hydro-simulations lies in the intrinsic difficulty to describe the distribution of baryons inside haloes, and therefore to trace the spatial properties of galactic structures, like discs and bulges.The coupling of the CR-IGIMF with a full fledged theoretical model of galaxy evolution is a key improvement with respect to our preliminary analysis in F18b, where we resort to simplified assumptions and idealized star formation histories to have an estimate of the IMF-sensitive spectral features.Indeed, the availability of a self-consistently calibrated model allows us to follow in detail the assembly of ETGs, by assigning a different IMF shape from our library to each SSP that describes the galaxy progenitors.This aspect is crucial in order to correctly establish the contribution of each IMF to the final model ETG spectra.
In order to estimate the IMF shape to be associated with a given star formation event, we use the integrated properties of each model galaxy at the corresponding cosmic time.While the definition of a galaxy-wide SFR is straightforward, the second parameter required by the CR-IGIMF, namely the CR energy density (U CR ) requires some extra discussion.In this work, we chose the same approach as in Fontanot et  focus on a model implementing the PP11 results), and we approximate the U CR by the SFR surface density (Σ SFR ), computed using the disc optical radius (i.e.3.2 times the exponential scale length).This implies that we are assuming U CR to be homogeneous over the entire galactic disc and proportional to its SFR surface density.The PP11 simulations are defined for U CR levels relative to the MW density field; in detail we assume: where we use the reference value of Σ MW = 10 −3 M ⊙ yr −1 kpc −2 for the MW disc.It is important to keep in mind that the two parameters regulating the shape of the CR-IGIMF are not independent in our approach.In particular, gaea model galaxies populate well defined regions of the SFR-Σ SFR space (as shown e.g. in Fig. 3 of Fontanot et al. 2018a).

gaea
In order to test the effect of the CR-IGIMF hypothesis on the physical and observable properties of galaxies, we implement it in the gaea model.2014); (ii) an updated prescription for ejective stellar feedback6 (Hirschmann et al. 2016, HDLF16); (iii) an updated scheme to predict disc size evolution (Xie et al. 2017).This choice allows a direct comparison with our previous results.Additional modules of the gaea model not considered in this paper include a treatment for the partition of cold gas in neutral and molecular hydrogen (Xie et al. 2017), an updated modelling of environmental effects (Xie et al. 2020) and of gas accretion onto Super-Massive Black Holes (Fontanot et al. 2020).We plan to study the influence of these recent updates in the variable IMF framework in future work.
The version of this model using a universal IMF has been shown to be able to reproduce a number of fundamental properties of galaxy populations, including the evolution of the space density of M ⋆ < 10 10.5 M ⊙ galaxies (HDLF16), the evolution of the galaxy stellar mass function (GSMF) and cosmic star formation rate at z < ∼ 7 (Fontanot et al. 2017b), the evolution of the mass metallicity relations (De Lucia et al. 2020;Fontanot et al. 2021).The code modifications required to deal with a variable IMF scenario mainly refer to the baryonic mass fraction locked into low-mass stars (affecting the total stellar mass and luminosity of model galaxies) and to the relative abundance of Type Ia and Type II SNe (responsible for the amount of metals and energy restored into the ISM).We refer the reader to F17 for a discussion of the technical details.

Runs & Calibration
In this work, we run gaea on merger trees extracted from the Millennium Simulation (Springel et al. 2005), that represents a realization of a cosmological 500 Mpc 3 volume assuming a ΛCDM concordance model (with WMAP1 parameters Ω Λ = 0.75, Ω m = 0.25, Ω b = 0.045, n = 1, σ 8 = 0.9, H 0 = 73 km/s/Mpc).We do not expect the differences between these values and more recent determinations (i.e.Planck Collaboration XVI 2014) to affect our conclusions (as shown by Wang et al. 2008;Guo et al. 2013).
A variable IMF has a relevant impact on the amount of energy released by SNe, the chemical patterns and on the amount of baryonic mass locked in low-mass, long-living, stars.As in our previous work, we calibrate our models primarily considering multiwavelength luminosity functions (LFs, Fig. 1).The choice of using LFs to calibrate the model derives essentially from the fact that in this scenario the GSMF is no longer a robust indicator of galaxy evolution, since most of the M ⋆ estimates available in the literature have been derived under the assumption of a universal, MW-like, IMF.As in our previous work, we use photometric tables corresponding to the binned IMFs, using an updated version of the Bruzual & Charlot (2003) models (see e.g.Plat et al. 2019), that employs stellar evolutionary tracks including the Marigo et al. (2008) prescriptions for the evolution of the thermally pulsating AGB stars and stellar libraries as specified in Appendix A of Sánchez et al. (2022).Dust extinction is modelled as in De Lucia & Blaizot (2007).We consider the same calibration set we used in F17 and based on the local SDSS g, r and i-band LFs, and on the redshift evolution of the K-and V-band LFs.In Fig. 1 we compare observational estimates with predictions of our reference CR-IGIMF model.ǫ eject , the reincorporation rate γ reinc and the Radio-mode AGN feedback efficiency (κ radio ).
As discussed in the previous paragraph, the assumption of a variable IMF requires an additional layer in the analysis, as it introduces a mismatch between the intrinsic stellar mass M ⋆ predicted by the model, and the observational datasets, which are usually derived assuming a universal IMF.Therefore, in order to provide a quantitative comparison between the predictions of our model and existing observational estimates we define an apparent stellar mass (M app ⋆ ), using the same approach as in F17 and F18a, that we summarize in the following.Since our synthetic photometry is derived self consistently from SSP models that consider the correct IMF of each star formation event, we can use this information to robustly estimate the stellar mass an observer would derive using magnitudes in different bands and assuming a MW-like IMF.In particular, we employ a theoretically derived7 relation between the mass-to-light ratio in the iband and the g − i colour: with υ = 0.9 and δ = 0.7 (Zibetti, private communication, based on an analysis similar to Zibetti et al. 2009).An additional factor ǫ = 0.13 has been introduced in F17 to account for spatial resolution effects (see F17a for a complete discussion on the origin and justification of this shift).

Results
We first consider model predictions against a set of observational results that have already been tested in F17 and F18a, in order to check the consistency of the new IMF variability scheme with our previous findings (Sec.4.1 and 4.2).We then move to explore additional spectroscopic constraints to validate our proposed variability scenario (Sec.4.3 and 4.4).

Stellar mass assembly
The resulting GSMFs at different redshifts, computed using the intrinsic M ⋆ and the apparent M app ⋆ are shown in Fig. 2: as shown and discussed in F17, the photometrically derived GSMF in a variable IMF scenario exhibits systematic deviations from the intrinsic one, in particular at low redshift (z<1), where the growth of massive structures is apparently slowed down considering M app ⋆ .At higher redshift, the main differences between the intrinsic GSMF and the one defined using M app ⋆ are found around the knee and at the low-mass end.
Intrinsic and apparent stellar masses for individual systems can be compared to assess the relevance of deviations from the MW-like IMF.In order to reproduce the same diagnostic as in observational samples, we consider bulge-dominated model galaxies, i.e. a bulge-to-total stellar mass ratio (B/T) larger than 0.7.We try to mimic the dynamical analysis of Cappellari et al. (2012) by considering the ratio of the true stellar mass-to-light ratio in the i-band (M ⋆ /L i ) and the apparent one (Eq.11), as a function of the proper M ⋆ /L i (right panel of Fig. 3).We also consider the M ⋆ /M app ⋆ ratio as a function of M ⋆ (left panel of Fig. 3), which should correspond to Fig. 4 in Conroy et al. (2013).The predictions for the CR-IGIMF run (blue solid lines) are in qualitative agreement with observational results showing an increase of the so-called α-excess with stellar mass and massto-light ratios.For reference, we also consider predictions of the same quantities from our standard run (HDLF16 -black dashed lines) adopting a universal, MW-like, IMF, which predicts flatter relations.These results confirm the conclusions of our previous work.As in F17, the positive α-excess at M ⋆ > 10 10 M ⊙ (M ⋆ /L i >0) is driven by variations of the IMF high-mass end (being shallower than Salpeter in high-SFR events connected with the formation of massive ETGs).In particular, we do not interpret the α-excess in massive ETGs as a result of a typical IMF bottom-heavier than the MW-like, but instead as due to the mismatch between proper and synthetic mass-to-light ratios.Nonetheless, it is worth noting that assuming that also the low-mass end can vary (as in our CR-IGIMF scenario) does not increase the maximum predicted values for the α-excess (and/or the normalization of the relation).We will deepen this point in Sec.4.3.The main difference with the standard IGIMF case is found at the low-mass-to-light ratio end of the right panel.In the IGIMF scenario, the constancy of the low-mass-end results in a flattening of α-excess (see i.e.Fig. 8 in F17), while the CR-IGIMF realization predicts an (almost) constant slope of the relation extending to negative α-excess values.

Metallicities
We then consider the predicted evolution of the mass-metallicity relations (MZRs, both for the stellar component and for the gaseous component -left and right panels in Fig. 4, respectively).We compare predictions from the CR-IGIMF run (blue solid line -the shaded area represents the 15-85th percentiles of the distributions) with the collection of observational estimates used in Fontanot et al. (2021, see that paper for a complete reference list).We also include predictions from our standard HDLF16 run (black dashed line).As discussed in Fontanot et al. (2021) we assume a 0.1 dex shift downwards of gas metallicities in order to get the right normalization with respect to z ∼ 0 data and better highlight the fact that our standard model is able to reproduce the redshift evolution of the gas MZR.It is worth noting that this assumed shift is still within the uncertainty in the normalization of the MZR (see e.g.Kewley & Ellison 2008).
The results shown in Fig. 4 clearly highlight that the CR-IGIMF run predicts an overall normalization of the MZR relations in good agreement with the observational estimates and with the universal IMF scenario: the main tensions are seen for gas metallicities of low-mass galaxies at z ∼ 0 and there is a clear overprediction of the stellar metallicities at z∼3.5 (which is already present in the HDLF16 run).Nonetheless, the CR-IGIMF run systematically predicts steeper MZRs at z > 1, both in the stellar and gaseous phases.The increase in metallicity is particularly relevant at M ⋆ > 10 10 M ⊙ , with model galaxies reaching (super-)solar metallicities even at considerable redshifts (z∼2-3).These results for more massive galaxies are in tension with observational estimates reported in Fig. 4. Nonetheless, some recent results suggest that, at least a fraction of the most massive galaxies at z ∼ 3.4 have already super-solar metallicities (Saracco et al. 2020), in qualitative agreement with the predictions of the CR-IGIMF scenario.Moreover, Saracco et al. (2023) recently find no metallicity evolution for massive ETGs over the redshift range 0.0 < z < 1.4 (see also Vazdekis et al. 1997, and in particular their Fig.21).
We also check that the CR-IGIMF run reproduces reasonably the [α/Fe]-stellar mass relation for Early Type Galaxies (that we select in model predictions again applying a B/T>0.7 criterion).In detail, we compare observational estimates for [α/Fe], calibrated using magnesium lines, against [O/Fe] theoretical predictions (since oxygen represents the most abundant α element).Fig. 5 shows that predictions for the CR-IGIMF run (blue solid line) are perfectly in line with our previous findings, showing that the IGIMF scenario naturally predicts an increase of [α/Fe] with stellar mass, at variance with the model assuming a universal, MW-like, IMF.In F17 (see their Fig.7), we show that the positive slope is a combination of (a) the increase of the αenhancement for massive galaxies, due to the shallower highmass end slope (with respect to α 2 ) associated with the high-SF events characterizing their early evolution and (b) the predominance of IMFs bottom-heavier that the MW-like in low-mass galaxies (leading to a decrease of the [α/Fe]).
Also in this case, as for the α-excess, there is no relevant change in the results obtained within the CR-IGIMF framework with respect to the standard IGIMF scenario.We thus confirm all our previous conclusions about the impact of a variable IMF on the assembly of galaxy stellar masses and metallicities: our interpretation still holds in the CR-IGIMF scenario.In the following we will now consider additional spectroscopic constraints.

Dwarf-to-Giant ratio
In order to see if the CR-IGIMF run can reproduce the observational finding of a bottom-heavy IMF in more, relative to less, massive galaxies at z ∼ 0, for each model galaxy, we compute the IMF at z ∼ 0, and the corresponding dwarf-to-giant ratio, f dg .The f dg is defined as the ratio of the mass of stars below 0.6 M ⊙ , to that of stars below 1.0 M ⊙ , in the IMF.The reason for this normalization is that stars more massive than ∼1 M ⊙ are already dead for old populations, and do not contribute to the integrated galaxy light (Martín-Navarro et al. 2019).The f dg has been found to increase with galaxy velocity dispersion, reflecting the variation of IMF slope La Barbera et al. (2013).
Fig. 6 plots the relation between f dg and velocity dispersion, σ, for model galaxies in the CR-IGIMF run.As in Fig. 3 of F18b, σ is estimated using the relation between velocity dispersion and stellar mass from Zahid et al. (2016).The green stars in Fig. 6.The ratio of dwarf-to-giant stars in the IMF, f dg , is plotted against the velocity dispersion of model galaxies, σ, for the CR-IGIMF run (black dots).The top and bottom panels are obtained by considering all stars in the model, and those formed in-situ, respectively.Green stars show the binned median trends, while the blue line corresponds to the IMF -σ relation from La Barbera et al. (2013).The orange circles with error bars show the latter relation corrected to an aperture of 1 R e , while the red line corresponds to an infinite aperture (see text).The expected f dg values for a Kroupa and Salpeter IMF are marked by dashed horizontal lines, while the dotted line marks the f dg for a MW-like IMF in the model.
the Figure show the median in bins of σ.We note that in the CR-IGIMF framework, the MW-like IMF (i.e.U CR /U MW = 1 and SFR∼ 1 M ⊙ yr −1 ) has f dg = 0.718, slightly above the value for a Kroupa IMF (f dg = 0.713, see the Figure ).We find that f dg tend to increase (by 0.01-0.02)with σ, approaching the value for a MW IMF at σ ∼ 80 km s −1 .The trend is consistent, although shallower, than that of F18b.Moreover, it is slightly more clear for in-situ stars8 , where f dg gets to a value as low as ∼ 0.72 at σ ∼ 80 km s −1 (in contrast to f dg ∼ 0.73, when all stars are considered).We stress that in the present work, at variance with the preliminary analysis presented in F18b, we are analyzing a fully calibrated, self-consistent, CR-IGIMF gaea realization.
In order to compare the predicted f dg -σ trend with observations, we consider the relation obtained by La Barbera et al. (2013) (see their section 7, for the "2SSP+X/Fe" fitting case), for the central regions of ETGs as observed by the SDSS (blue curve in Fig. 6).This goes from a Kroupa-like distribution at σ < ∼ 170 km s −1 , to an f dg consistent, or even steeper, than that corresponding to a Salpeter IMF (f dg ∼ 0.84), at σ ∼ 300 km s −1 .In order to perform a fair comparison with gaea, we have to consider that the model predicts global quantities for each galaxy, while observations refer to the SDSS aperture.To perform an aperture correction, we assume an IMF-radial profile as in La Barbera et al. (2017), where the IMF slope is bottom-heavy in the center, and decreases to a Kroupa-like distribution at about the effective radius (R e ).This trend is also consistent with the IMF radial profiles derived by van Dokkum et al. (2017).The IMF profile is normalized in each σ bin to match the IMF slope within the SDSS fiber aperture, and then used to recompute the IMF slope within a given aperture, and the corresponding f dg value.The orange points and red lines in Fig. 6 show the f dg -σ relations corrected to 1 R e , and an infinite aperture9 , respectively.As expected, the relation flattens as one moves to a larger aperture.For an infinite aperture, the range of predicted values for f dg (green stars) is similar to that observed.However, the model trend is still shallower than the observed one: at high sigma (σ ∼ 300 km s −1 ), the model tends to underpredict f dg , while at 100 < σ < 220 km s −1 , most model galaxies have f dg above the observed values.Therefore, the present CR-IGIMF implementation appears to be unable to quantitatively reproduce the IMF-σ relation observed at z ∼ 0.

IMF-sensitive indices at z∼0
In order to further compare predictions of the CR-IGIMF run with observations, we compute synthetic spectra and spectral indices for gaea model galaxies.First, using the spectral synthesis code of Vazdekis et al. (2012), we compute a large library of SSP models, with varying age, metallicity, and IMF.For the present work, we generate only models in the optical spectral range, relying on the MILES stellar library (3540 < λ < 7410 Å; Sánchez-Blázquez et al. 2006), at a spectral resolution of 2.5 Å (FWHM; Falcón-Barroso et al. 2011).We note that the spectra use a different spectrophotometric synthesis code with respect to the one used to calibrate the model (i.e. the updated version of BC03, see above).Nonetheless, since we use broad-band photometry for calibration, we do not expect significant differences between the two codes, that could affect our conclusions.Generating synthetic spectra with the Vazdekis code allows us to perform a more consistent comparison to results obtained from observed spectra using the same code (see e.g.La Barbera et al. 2013).On the other hand, using BC03 for calibration, we can have a direct comparison with our previous results (F17 and F18a).We construct SSPs based on Teramo isochrones, with ages in the range from 0.03 to 14 Gyr, and metallicities from [M/H] = −2.27 to +0.26 (see Vazdekis et al. 2015 for details).For each age and metallicity, SSPs are computed with the 42 different IMFs defined in Table 1, corresponding to seven (six) values of log(SFR) (U CR /U MW ).For each model galaxy, we store its full star formation history, consisting of the list of mass elements that form in all galaxy's progenitors at different time steps of the simulation.In order to disentangle the effect of a variable IMF from that of the SFH in gaea, we fix, for simplicity, the age and metallicity of all elements to their mass-weighted values over the galaxy SFH.This should also allow us to perform a more fair comparison to observations, where usually the IMF is inferred by assuming a single SSP.For each element, we compute, by linear interpolation, the corresponding synthetic spectrum, given its age, metallicity, and IMF.The spectrum of each element is normalized to its stellar mass, and all elements are summed up to produce the final synthetic spectrum of the galaxy.run, the former with an f dg of ∼ 0.75 (i.e.heavier than MW), and the latter with a MW-like IMF (f dg ∼ 0.71).The bottom panel of the Figure shows the ratio of the two spectra to those obtained by assuming a MW IMF.For the low-mass galaxy, as expected, the ratio is close to one, while for the high-mass system we see features typical of a more bottom-heavy IMF (see La Barbera et al. 2013), such as stronger NaD (λ ∼ 5900 Å) and TiO2 SDSS (λ ∼ 6200 Å) absorptions, and weaker CaH+K lines (λ ∼ 3900 Å).However, the effect is small, and variations are below ∼ 1 % for most of the spectrum.This is consistent with what found in F18b based on a preliminary implementation of the CR-IGIMF framework.
In order to look in more detail at absorption features, we compute spectral indices from the synthetic spectra of all model galaxies.We focus here on two prominent IMF-sensitive features in the optical spectral range, i.e.TiO2 SDSS and NaD, whose definitions, i.e. central passband and pseudocontinua, are taken from La Barbera et al. (2013) and Trager et al. (1998), respectively.Fig. 8 plots the difference of TiO2 SDSS (top; δ(TiO2 SDSS )) and NaD (bottom; δ(NaD)) with respect to the case of a MW IMF, as a function of the galaxy velocity dispersion.Note that, following the original definition of the Lick indices, NaD is computed in units of Å, while TiO2 SDSS is given in magnitudes.Green stars represent median-binned trends for model galaxies, while red lines show the observed trends for the SDSS stacked spectra of La Barbera et al. (2013).The observed values of δ(TiO2 SDSS ) and δ(NaD) correspond to the case of an infinite aperture (see Fig. 6), and were obtained as follows.For a given stacked spectrum, we compute the TiO2 SDSS and NaD indices for an SSP model with the age, metallicity, and IMF slope measured for the given spectrum (see La Barbera et al. (2013) for details).The IMF slope is corrected to an infinite aperture (see Sec. 4.3), while for age and metallicity we consider the values measured within the SDSS fiber aperture.Then we compute the difference of these line-strengths with respect to the case where the IMF is set to a MW-like distribution, while age and metallicity are not varied.These differences give the red curves in Fig. 8.In the model, both NaD and TiO2 SDSS tend to increase with σ, as expected from the behaviour of f dg .However, the predicted variation is too small compared to observations.For instance, at the highest σ (∼ 300 km s −1 ), the observed δ(NaD) is ∼ 0.12 Å , while for model galaxies, it is, on average, only ∼ 0.02 Å.In order to understand this discrepancy, we should keep in mind that the strength of IMF-sensitive features in the optical is mostly driven by f dg , but, to second order, it is also affected by the detailed shape of the IMF.For instance, using different IMF parametrizations (unimodal and low-mass tapered IMFs), LB13 showed that optical indices imply a similar trend of f dg with galaxy velocity dispersion, with small differences, at fixed sigma, between different parametrizations (see fig. 19 of LB13).Indeed, we found that the difference between the observed trends in Fig. 8 (red curves) and model predictions (green stars) is twofold.At σ ∼ 300 km s −1 , the average f dg of the CR-IGIMF run is ∼ 0.745, i.e. lower than the observed value of ∼ 0.76 (for the case of an infinite aperture, see Fig. 6).However, assuming f dg = 0.745 and the same IMF parametrizations10 as in La Barbera et al. (2013), and following the same procedure adopted to derive the red curve in Fig. 8, we would infer a value of δ(NaD)∼ 0.06 Å, which is smaller than the estimate obtained from SDSS spectra, but still significantly larger than the value of ∼ 0.02 Å in the model (green stars in the Figure ).Therefore, we conclude that the discrepancy seen in Fig. 8 is due to the fact that both the f dg -σ relation is too shallow in the models, compared to the data, and that the IMF shape in the CR-IGIMF framework is not able to produce a significant variation of the spectral features, as found in the data.We point out that the CR-A&A proofs: manuscript no.cr_igimf.finIGIMF framework does not assume any direct scaling of IMF with σ (or galaxy mass).As already discussed in F18b, the dependence on σ results from the competing effect of the increase of SFR with galaxy mass (which leads to a larger f dg at fixed U CR , see e.g.lines 4 to 7 or lines 11 to 14 in Tab. 1) and the U CR increase with galaxy mass (which favours lower f dg values, compare e.g.lines 12, 19 and 26 in Tab. 1).The fact that the trends of δ(TiO2 SDSS ) and δ(NaD) with σ in Fig. 8 are too weak compared to the data, implies that some further ingredient is still missing in the current implementation of the varying IMF framework.One possible solution would require an explicit scaling with σ, e.g., assuming a variation of the low-mass end of the IMF as a function of some physical parameter, such as metallicity (e.g.Jeřábková et al. 2018;Yan et al., in preparation) and/or surface density (see Tanvir & Krumholz 2023), consistent with observational results (e.g.Martín-Navarro et al. 2015, La Barbera et al. 2019).Since more massive galaxies are also more metal-rich and denser than their low-mass counterparts, we expect this would produce, by construction, a stronger increase of δ(TiO2 SDSS ) and δ(NaD) with σ, as well as a stronger trend of f dg with galaxy mass.Nonetheless, such implementation goes beyond the scope of the present work: we defer a more throughout investigation of the effect of these additional dependencies to future work.

Conclusions
In this work, we present predictions for an updated version of the state-of-the-art semi-analytic model gaea, which implements our recently proposed variable IMF framework (CR-IGIMF).The CR-IGIMF scenario combines the Integrated galaxy-wide IMF approach (see e.g.Kroupa et al. 2013) with the results of high-resolution numerical simulations of star formation in giant molecular clouds (Papadopoulos et al. 2011): the main improvement with respect to the standard IGIMF approach, lies in the fact that we allow individual MCs to be characterized by different IMFs, under the assumption that their characteristic Jeans mass depends on the Cosmic Ray background they are embedded in.As a result of our calculations, we show in Fontanot et al. (2018b) that this approach is able to predict CR-IGIMF shapes that deviate from a universal, MW-like, IMF at the high-and low-mass ends at the same time.We deem this property very promising to explain in a self-consistent way several observational findings of a non-universal IMF, that are usually difficult to reconcile under a single top-heavy or a bottomheavy IMF scenario.The coupling of the CR-IGIMF with gaea allows us to exploit the features of our galaxy formation model (and in particular its improved modelling of chemical enrichment -De Lucia et al. 2014) to consider a variety of predictions for the physical, dynamical, chemical and photometric properties of model galaxies.For the purposes of this study we also couple gaea with the spectral synthesis code by Vazdekis et al. (2010) in order to derive synthetic spectral energy distributions, that self-consistently account for IMF variation, thus allowing a direct comparison of synthetic spectral features with observed ones.
We can summarize our main conclusions as follows.
-The CR-IGIMF scenario confirms our previous findings (F17 and F18a): the expected variability at the IMF high-mass end allows us to match the so-called α-excess, i.e. the reported excess, with respect to expectations based on a universal IMF, of stellar mass and mass-to-light ratio derived from dynamical and spectroscopic studies.Moreover, our models ∼ 1, both the predicted stellar and cold-gas MZR are steeper than observational estimates.In particular, our model predicts substantial metallicities at the high-mass end (i.e.M ⋆ > 10 10 M ⊙ ) already at intermediate redshift, in tension with the observed relations.Nonetheless, these predictions may explain recent findings of the super-solar metallicities for a sample of massive ETGs at z < ∼ 3 (Saracco et al. 2023;D'Ago et al. 2023;Bevacqua et al. 2023).
-The properties of synthetic spectral energy distributions for model ETG galaxies are in qualitative agreement with results derived from observed spectra.In particular, our synthetic spectra show a similar increasing trend of the dwarfto-giant ratio (f dg , i.e. the mass fraction in the IMF of stars with masses below 0.6 M ⊙ ) as a function of galaxy velocity dispersion.Moreover, also IMF-sensitive synthetic spectral features, like TiO2 SDSS and NaD show a slightly increasing trend with velocity dispersion.However, all these trends, despite showing the correct dependence with velocity dispersion, are considerably weaker than observed, possibly demanding further developments of the CR-IGIMF framework.
Overall, our CR-IGIMF gaea realization is able to reproduce at the same time a variety of observational findings, based on dynamics, and photometric/spectroscopic studies suggesting a time variability of the IMF of ETGs galaxies.This represents a fundamental step forward in our understanding of the overall impact of a variable IMF on galaxy properties, and in constraining variability scenarios.The main problem of our current model lies in the fact that it is not able to reproduce the correct strength in the evolution of spectral features such as TiO2 SDSS and NaD with velocity dispersion, although the predicted relations trend in the correct direction.It is important to stress that this is not a limitation of the IGIMF framework (and its extension such as the CR-IGIMF), but it is more tightly connected with the underlying assumptions.For example, our CR-IGIMF library is limited to low-mass end slopes that range between a Kroupa and a Salpeter IMF by construction, i.e. we assume that the IMFs associated with individual clouds cannot be steeper than a Salpeter IMF or shallower than a Kroupa IMF below 1 M ⊙ .Our results shows that this might not be sufficient to reproduce the observed f dg trends as a function of velocity dispersion, although they typically range between these two extremes (see e.g.Fig. 6).In order get steeper slopes both for f dg and spectral indeces, we likely have to act at the level of IMF slopes in individual clouds.For example, more extreme IMF, super-Salpeter α 1 ∼ −3 slopes have been suggested in La Barbera et al. (2016) to explain radial IMF trends in the innermost regions of local ETGs.The origin of such extreme IMF slopes depends on the physical conditions of star forming regions (see e.g.Chabrier et al. 2014), that may deviate largely from local counterparts in the MW environment, in particular for what concerns molecular gas density and metallicity.We will explore the impact of a wider range of low-mass end slopes in future work.

Fig. 2 .
Fig.2.Galaxy stellar mass function evolution, as derived from different stellar mass estimates.The blue solid lines correspond to the estimate using the intrinsic stellar masses from the gaea run.Red dashed lines show the predicted evolution using M app ⋆ , i.e. the stellar mass derived from synthetic photometry under the hypothesis of a universal, MWlike, IMF (see text for more details.Grey points represent the observational datapoints as in the compilation fromFontanot et al. (2009).
Fig.3.α-excess as a function of stellar mass (Left-hand panel) and stellar mass-to-light ratio (Right-hand panel).In each panel, blue solid and black dashed lines correspond to the predictions of the CR-IGIMF and the universal, MW-like, IMF runs, respectively; light blue contours mark galaxy number densities levels (normalized to the maximum density) corresponding to 1, 5, 10 and 25 per cent, in the CR-IGIMF run.

Fig. 4 .
Fig. 4. Redshift evolution of the cold gas (left panel -a downwards 0.1 dex shift has been applied to the predictions to match the overall normalization -see text for more details) and stellar (right panel) MZR relations.Lines and shaded areas refer to the CR-IGIMF and HDLF16 models as in Fig.5.Datapoints correspond to the compilation used inFontanot et al. (2021).

Fig. 7 .
Fig.7.Examples of synthetic MILES spectra (top panel) for two model galaxies in the CR-IGIMF run of gaea.The black and grey curves correspond to a high-(M ⋆ ∼ 11.3) and low-mass (M ⋆ ∼ 9) galaxy, respectively.The bottom panel shows the ratio of each spectrum to that obtained for a MW IMF.Both spectra are at the nominal resolution of MILES (∼2.5 Å FWHM).

Fig. 8 .
Fig.8.Difference of spectral indices, between the case of variable-and a MW-IMF, for model galaxies in the CR-IGIMF run (black dots).The top and bottom panel are for the TiO2 SDSS and NaD indices, respectively.Green stars show the binned trend, with error bars marking the 5-th and 95-th percentiles of the distribution in each bin.The red curves correspond to the observed trends from LaBarbera et al. (2013), corrected to an infinite aperture (see text for more details).
Table 2 contains the values of the relevant parameters used in the recalibration (see HDLF16 for a detailed discussion): the SFR efficiency α SF , the stellar feedback reheating efficiency ǫ reheat and ejection rate A&A proofs: manuscript no.cr_igimf.fin