Fabio Fontanot

Gabriella De Lucia

Lizhi Xie

Michaela Hirschmann

Gustavo Bruzual

Stéphane Charlot

Variations of the stellar initial mass function in semi-analytical models II: the impact of Cosmic Ray regulation

Keywords: galaxies, formation -galaxies, evolution -galaxies, abundances -galaxies, fundamental parameters -galaxies, stellar content

Recent studies proposed that cosmic rays (CR) are a key ingredient in setting the conditions for star formation, thanks to their ability to alter the thermal and chemical state of dense gas in the UV-shielded cores of molecular clouds. In this paper, we explore their role as regulators of the stellar initial mass function (IMF) variations, using the semi-analytic model for GAlaxy Evolution and Assembly (GAEA). The new model confirms our previous results obtained using the integrated galaxy-wide IMF (IGIMF) theory: both variable IMF models reproduce the observed increase of α-enhancement as a function of stellar mass and the measured z = 0 excess of dynamical mass-to-light ratios with respect to photometric estimates assuming a universal IMF. We focus here on the mismatch between the photometrically-derived (M app ⋆ ) and intrinsic (M ⋆ ) stellar masses, by analysing in detail the evolution of model galaxies with different values of M ⋆ /M app ⋆ . We find that galaxies with small deviations (i.e. formally consistent with a universal IMF hypothesis) are characterized by more extended star formation histories and live in less massive haloes with respect to the bulk of the galaxy population. While the IGIMF theory does not change significantly the mean evolution of model galaxies with respect to the reference model, a CR-regulated IMF implies shorter star formation histories and higher peaks of star formation for objects more massive than 10 10.5 M ⊙ . However, we also show that it is difficult to unveil this behaviour from observations, as the key physical quantities are typically derived assuming a universal IMF.

INTRODUCTION

Star formation represents a key, yet not fully understood, physical mechanism governing galaxy evolution. In a star formation episode, one of the fundamental parameters of the associated stellar population is the so-called Initial Mass Function (IMF), i.e. the number of stars formed per stellar mass bin. The IMF has long been seen as a nearly invariant property of star forming regions, mainly due to the remarkable consistency of its shape measured in the local neighbourhood (but apparently not in the densest regions of the Galactic Centre, see e.g. Klessen et al. 2007). Different functional forms have been proposed to describe the IMF (Salpeter 1955;Kroupa 2001;Chabrier 2003, among the most popular). These differ mainly for the abundance of stars at the low-mass end up to the brown dwarf regime. Unfortunately, direct measure-⋆ E-mail: fontanot@oats.inaf.it ments of the IMF via star counts is accessible only for local star forming regions (typically within our own Milky Way or its closest satellites), and we are thus forced to resort to the analysis of the integrated light of stellar populations in external galaxies. Indirect evidence against a universal 1 IMF has been found for different galaxy populations, such as late-type (Gunawardhana et al. 2011), earlytype (Ferreras et al. 2013) and dwarf galaxies (McWilliam et al. 2013). These claims have been recently confirmed by dynamical and spectroscopic studies suggesting relevant deviations from a universal IMF as a function of galaxy stellar mass or velocity dispersion. For a sample of early-type galaxies in the ATLAS 3D survey, Cappellari et al. (2012) contrasted integral-field maps of stellar kinematics and optical imaging against dynamical models, including the contribution of both stellar and Dark Matter (DM) components. Their results show a systematic excess of the dynamical mass-to-light ratios with respect to photometric estimates based on a universal IMF. This excess increases with galaxy velocity dispersion (σ). The interpretation of this result in terms of IMF variations is not straightforward. A possible solution is given by the assumption that massive galaxies are characterized by a "top-heavy" IMF, which implies a larger fraction of massive (short-lived) stars and a decrease in the total light expected at late times (at fixed stellar mass), with respect to a universal IMF assumption. Conversely, a "bottom-heavy" solution is also plausible: it corresponds to a larger fraction of low-mass stars and an increase of the stellar mass at fixed light, with respect to a universal IMF. Cappellari et al. (2012, see also Tortora et al. 2016) were unable to distinguish between the two scenarios. Compatible results were obtained by Conroy & van Dokkum (2012, see also Ferreras et al. 2013 andSpiniello et al. 2012) using a spectrophotometric approach. Their analysis uses spectral features in galaxy spectra sensitive to the stellar effective temperature and surface gravity, and thus able to constrain the ratio between low-mass and giant stars, which gives information on the IMF shape. In particular, they compare high-resolution spectra of a sample of compact early type galaxies against libraries of synthetic models obtained by varying the IMF slopes at the high-mass and low-mass ends. The pre-selection of compact elliptical galaxies allows them to assume that the observed σ within the effective radius is dominated by stellar kinematic (neglecting the DM contribution). They interpret their best fit solutions in favour of an increasingly bottom-heavy IMF with increasing σ and/or stellar mass.

From a theoretical point of view, it has long been argued that the IMF could and should vary as function of local properties of the interstellar medium (ISM) in star forming regions. A number of theoretical models (see e.g. Klessen et al. 2005;Weidner & Kroupa 2005;Hennebelle & Chabrier 2008;Papadopoulos 2010;Hopkins 2012) have analyzed the role of small-scale physical properties of the ISM in star forming regions in setting the IMF shape, and predict a wide range of possible IMF shapes. They differ mainly in the modelling of the physical processes that govern the fragmentation and formation of stars of different mass from the parent molecular cloud (MC). As we already mentioned, it is currently impossible to study the IMF shape directly (i.e. via star counts), but we typically constrain it using the integrated light of composite stellar populations in external galaxies. This complicates the picture, as our analysis of a composite stellar population depends on our understanding of the simple stellar populations which form at different times and locations (see e.g. Kroupa & Weidner 2003;Kroupa et al. 2013). As many physical properties of galaxies change both spatially and temporally, we expect the redshift evolution of the IMF shape to play a fundamental role in setting the physical properties of galaxies. These considerations suggests that available data can provide us only instantaneous (often galaxy-wide) information on a "mean" or "effective" IMF. Dedicated theoretical modelling is crucial to test the proposed models, following the build up of the composite stellar populations of different galaxy classes and comparing their synthetic properties with observational constraints. Using this approach, it will be possible to simulate observational samples, prepare mock catalogues and characterize selection effects. This is particularly relevant, given the current claims of possible inconsistencies between the dynamical and spectroscopic approach (see e.g. Smith 2014).

In a previous paper, Fontanot et al. (2017, F17a hereafter, see also Fontanot 2014) we tested the impact of the so-called Integrated Galaxy-wide stellar Initial Mass Function (IGIMF) theory, first proposed by Kroupa & Weidner (2003, see also Weidner & Kroupa 2005), on predictions from our GAlaxy Evolution and Assembly (GAEA) model. This model is particularly well designed for this kind of studies as it features (i) a detailed chemical enrichment model (De Lucia et al. 2014), following the evolution of several species and abundance ratios and (ii) an improved treatment for stellar feedback, gas ejection and re-accretion (Hirschmann et al. 2016) that reproduces the redshift evolution of the galaxy stellar mass function over a wide range of stellar masses and redshifts (Fontanot et al. 2017). The IGIMF theory directly relates the mean shape of the IMF to the instantaneous star formation rate (SFR) of each galaxy, so that it is well suited for theoretical models not resolving the details of star forming regions in galaxies. The IGIMF variation as a function of SFR corresponds to a change of the highstellar mass slope, with stronger SFR episodes corresponding to top-heavier IMFs.

In this work, we will study the implications of an independent approach, first proposed by Papadopoulos (2010) and later developed in Papadopoulos et al. (2011, P11 hereafter). The role of Cosmic Rays (CR), associated with supernovae (SNe) explosions and stellar winds around young stars, in the process of star formation has been long debated in the literature Krumholz (see e.g. 2014), given their higher efficiency in penetrating the inner dense, UV-shielded regions of MCs, with respect to UV-photons. As an intriguing possibility, PP11 postulate that the local CRs energy density could indeed regulate the minimum temperature, ionization state and chemical composition of these dense MC cores and, ultimately, their characteristic Jeans mass. This is a quite different assumption with respect to the IGIMF theory, that postulates that individual molecular clouds are characterized by a canonical (Kroupa-like) IMF. Nonetheless, it is worth stressing that the effect of CRs should be implicitly accounted for in the IGIMF theory, as it is based on an empirical calibration of the stellar IMF and its dependence on local properties (Marks et al. 2012;Kroupa et al. 2013). The predicted evolution of the IMF shape is however sensibly different from the IGIMF theory: in the PP11 formalism, the main evolution is seen in the characteristic mass (i.e. the position of the knee of the IMF). These two models represent two independent ways of varying the IMF shape: a comparison of their predictions will thus provide insight on their different implications in terms of star formation history and mass assembly.

This paper is organized as follows. In Section 2 we will outline the basics of the CR theory as presented in PP11 and [START_REF] Papadopoulos | Space Science Proceedings, The Initial Conditions of Star Formation: Cosmic Rays as the Fundamental Regulators[END_REF]. We will then describe its semianalytic implementation in Section 3. We will present our results and compare them with the IGIMF expectations in Section 4. Finally, we will summarise and discuss our results in Section 5.

COSMIC RAYS AS REGULATORS OF STAR FORMATION

In order to estimate the IMF variations as a function of the CR density, we follow the approach described in PP11 (see also [START_REF] Papadopoulos | Space Science Proceedings, The Initial Conditions of Star Formation: Cosmic Rays as the Fundamental Regulators[END_REF] cluding the CR heating term (ΓCR) in a thermal balance equation (Goldsmith 2001):

ΓCR = Λ line (T k ) + Λ gd (T k ) (1)
where Λ line and Λ gd represent the contribution to gas cooling from molecular lines and gas-dust interactions 2 , and ΓCR depends on the MC gas core density (nH2) and on the CR ionization rate (γCR).

ΓCR ∝ nH2γCR (2)
Following PP11, we assume that γCR is proportional to the CR energy density UCR. Using a simple analytic model based on a fixed chemical composition, PP11 showed that the minimum temperature in dense MC gas cores is much higher in galaxies with SFR densities larger than our own Milky Way. A more general solution for the chemical and thermal equations regulating the CR-dominated ISM requires the modelling of the evolution of ISM chemistry, which is also strongly dependent on γCR; this solution can be determined numerically (see e.g. Thi et al. 2009). PP11 discussed the evolution of the gas temperature as a function of CR energy densities T k (ΓCR) and MC gas core densities. The characteristic Jeans mass of young stars (M ⋆ J ) can be written as:

M ⋆ J = kB T k (ΓCR) GµmH2 3/2 ρ -1/2 c = 0.9 T k (ΓCR) 10K 3/2 nH2 10 4 cm -3 -1/2 M⊙. (3)
In PP11 there is no attempt to estimate a global IMF shape (i.e. to be associated with the whole sta forming galaxy): their results hold for individual star forming MCs, given their gas and CR densities. In order to translate these results into a prescriptions to be implemented in GAEA, we thus need some additional assumptions.

First of all we assume that the overall IMF shape is well described by a broken power law (see e.g. Kroupa 2001), characterized by two slopes:

ϕ⋆(m) = ( m m low ) -α 1 m low m < m break ( m 0 m low ) -α 1 ( m m 0 ) -α 2 m break m < mmax (4)
where m low = 0.1M⊙, mmax = 100M⊙, α1 = 1.3, α2 = 2.35. Such a general shape is motivated by local observations of individual clouds, and it agrees well with theoretical calculations based on fragmentation of giant molecular clouds (see e.g. Hennebelle & Chabrier 2008, and reference herein). In our models, we will assume that the high-and low-mass ends of the IMF have fixed slopes, and that only the characteristic mass of young stars, i.e. the knee m break of the IMF is affected by the CR density. At a given CR density (i.e. T k (ΓCR)), we assume that m break corresponds to M ⋆ J at n(H2) = 10 5 cm -3 , a representative density for typical molecular clouds (Papadopoulos, private communication). In detail, we adopt the numerical solutions by PP11 (their Fig. 4, see also Tab. 1), where M ⋆ J is computed as a function of n(H2) for 6 different values of CR energy density (normalized to the MW value, UCR/UMW). We approximate the CR energy density in our models galaxies by the SFR surface density (ΣSFR). Following PP11, we define 6 different IMF shapes (Fig. 1) assuming:

UCR UMW = ΣSFR ΣMW (5)
2 Explicit expressions for each cooling term can be found in Appendix A of PP11

Table 1. Analytic description for the IMF variations where ΣMW = 10 -3 M⊙ yr -1 kpc -2 is the estimated SFR surface density for the MW disc. It is important to keep in mind that this description of the IMF variation has not been tested yet against resolved stellar populations, such as globular clusters and/or young open clusters. In fact, these observations suggest an invariant break mass. The main line of evidence in favor of cosmic ray regulated IMF variations come from the properties of starburst galaxies [START_REF] Papadopoulos | Space Science Proceedings, The Initial Conditions of Star Formation: Cosmic Rays as the Fundamental Regulators[END_REF]Romano et al. 2017). CRs may also shift the average stellar mass in a population rather than only the break mass and may thus partially or wholly be responsible for the changes in the stellar IMF deduced by Marks et al. (2012).

U CR /U MW m low α 1 m break α 2 mmax [M ⊙ ] [M ⊙ ] [M ⊙ ] 0.2 0.

SEMI-ANALYTIC MODEL

In order to assess the influence of the proposed variable IMF model on the chemical and physical properties of galaxy populations, we implement it in the GAlaxy Evolution and Assembly (GAEA) semianalytic model (SAM). This method follows the evolution of galaxy properties along cosmic epochs by modelling the network of relevant physical mechanisms acting on the baryonic component of dark matter haloes, whose hierarchical evolution (traced using Nbody cosmological simulations) represents the gravitational backbone of the model. These baryonic processes (which include the cooling and heating of baryonic gas, star formation, accretion of gas onto Super-Massive Black Holes and the related feedback processes) are described using approximated analytic and/or numerical prescriptions that are observationally and/or theoretically motivated: this results in a very flexible tool to predict galaxy properties for galaxy samples in cosmological volumes.

Our GAEA model provides significant improvements with respect to previous versions of the code (De Lucia & Blaizot 2007). In this section we will quickly sketch the relevant changes implemented for the present analysis, and we refer the interested reader to the original papers for more details. GAEA implements a detailed treatment of chemical enrichment, described in De Lucia et al. (2014), which accounts explicitly for the different lifetimes of stars of different mass and their peculiar enrichment patterns. The differential lifetimes and yields for individual chemical species from Asymptotic Giant Branch (AGB) stars, Type II SNe and Type Ia SNe are explicitly tracked by the code.

A second relevant change in the code corresponds to the updated modelling of stellar feedback discussed in Hirschmann et al. (2016). In that paper, we show that some form of preventive or ejective feedback is needed to reproduce the observed space density evolution of galaxies below the knee of the stellar mass function. Our reference feedback scheme (H16F in the following) corresponds to a prescription combining (i) gas reheating and energy injection schemes mutuated from high-resolution hydrodynamical simulations (Muratov et al. 2015, i.e.), with (ii) a gas ejection 2011) and (iii) a timescale of gas re-incorporation that depends on halo mass (Henriques et al. 2013). In Hirschmann et al. (2016), we show that this model is also able to reproduce the gas fractions and mass metallicity relations at z < 3, while Fontanot et al. (2017) demonstrated that the agreement of the H16F run with the evolution of the stellar mass function and cosmic star formation rate extends to z ∼ 7 -10. It is worth stressing that the H16F feedback prescriptions do not include any explicit dependence on the number of SNe. However, in a variable IMF scenario, the fraction of SNe per unit stellar mass formed (fSN) is not constant. In order to account for this effect we compute fSN for the 6 discrete IMF shapes considered in this paper (see below). For each SF event, we then rescale the feedback efficiencies with the ratio between the fSN of the chosen IMF and the MW-like IMF (i.e. UCR/UMW = 1).

In the following, we will refer to predictions of this model as the PP11 run3 .

We compare these predictions with already published results from F17a (namely the "High-αSF model in F17a -IGIMF run in this paper). In F17a we introduced the code modifications needed to deal with a variable IMF. Those include the generalisation of lookup tables used by the code to track the enrichment of individual chemical elements and the energy given by each Single Stellar Population at different cosmic epochs. As in F17a, we construct photometric tables corresponding to the binned IMFs in Tab. 1, using an updated version of the Bruzual & Charlot (2003) model (Bruzual & Charlot, in preparation). This includes a prescription for the evolution of thermally pulsing AGB stars (Marigo et al. 2008). The treatment of dust extinction is the same as in De Lucia & Blaizot (2007).

In addition to these changes, the variable IMF runs we develop in this paper require further improvements in the code. Star formation in GAEA takes place in galaxy discs, with only a minor contribution from collisional starbursts in mergers (Fontanot et al. 2013). In the PP11 run, we take advantage of the updated modelling of disc sizes for GAEA described in Xie et al. (2017, see also Guo et al. 2011). While gas and stellar discs grow in mass, the code tracks their angular momentum evolution following the mass and energy exchanges between the different galaxy components (halo, disc and bulge). Xie et al. (2017) showed that the updated model predicts gas and stellar disks larger than the original H16F run, especially for more massive galaxies, in better agreement with observational measurements. Overall, the updated disc size model does not affect dramatically the predictions of the H16F run, and in particular statistical properties like the stellar mass functions and the mass-metallicity relations.

In this paper, we will use GAEA predictions based on the merger trees extracted from the Millennium Simulation (Springel et al. 2005), i.e. a ΛCDM concordance cosmology, with parameters4 derived from WMAP1 (i.e. ΩΛ = 0.75, Ωm = 0.25, Ω b = 0.045, n = 1, σ8 = 0.9, H0 = 73 km/s/Mpc). We will consider three different models, i.e. the original H16F model (as in Hirschmann et al. 2016), a model implementing the IGIMF theory, and the new model based on the PP11 approach. All these runs implement the same feedback scheme (as in H16F). However, the changes in galaxy evolution induced by a variable IMF approach are such that a recalibration of the key GAEA parameters is needed.

As for the IGIMF runs, we choose to recalibrate our model by requiring it to reproduce the evolution of multiwavelength luminosity functions. In fact, we cannot use the evolution of the galaxy stellar mass function (GSMF) and/or cosmic SFR as calibration set, since the observational estimates for these physical quantities are derived under the assumption of a universal IMF. We give more details on the recalibration process in Appendix A. The values of the relevant parameters are compared to the H16F and IGIMF parameters in Tab. 2.

In next sections, we will compare GAEA predictions with observational constraints involving physical properties of galaxies (like M⋆ or SFRs). These are usually derived from broadband photometry (spectral energy distribution fitting or colour scalings) or spectroscopy under the assumption of a universal/invariant IMF. In order to perform a proper comparison with our model predictions, we follow the same approach as in F17a and define an apparent -Chabrier ( 2003) IMF equivalent -stellar mass (M app ⋆ ) using a mass-to-light vs colour relation calibrated on synthetic magnitudes (see F17a and Appendix B for more details). In brief, M app ⋆ represents the stellar mass that an observer would derive from the synthetic photometry under the assumption of a universal IMF, while M⋆ is the intrinsic stellar mass predicted by the model. Both M⋆ and M app ⋆ account for the mass of stellar remnants.

RESULTS

In this section we will compare predictions based on the PP11 approach, with our previous runs implementing either the IGIMF framework or the canonical IMF (H16F). We aim to explore common trends between the variable IMF scenarios and highlight discrepant predictions that can be used to disentangle between the different scenarios. The predicted evolution of key physical and chemical properties for model galaxies in the PP11 run confirm the main conclusions by F17a in the framework of the IGIMF theory. Namely, the inclusion of a variable IMF provides a viable explanation for some recent and puzzling observational results and long standing problems. These models are able to recover, in particular, the observed relation between the [α/Fe] enrichment and the stellar mass of elliptical galaxies, and the correct trend for the z = 0 mass-metallicity relation. De Lucia et al. (2017) discuss the difficulties in reproducing both relations at the same time in theoretical models of galaxy formation (both SAMs and hydrodynamical simulations) using a universal IMF. Moreover, both IMF variation models are able to qualitatively reproduce the observational evidence in favour of a non-universal IMF, coming from dynamical analysis (see e.g. Cappellari et al. 2012) or spectral synthesis models (see e.g. Conroy & van Dokkum 2012 or La Barbera et al. 2017). The prediction of the IGIMF and PP11 models are quite similar for these observables, and hardly distinguishable within the errorbars. The interested reader will find the plots corresponding to these predictions in Appendix B, together with a more detailed discussion. In the following we will focus on the main differences we find between the two models and the reference H16F run.

In Fig. 2, we compare the GSMFs at different redshifts in the H16F, IGIMF and PP11 runs. In each panel, coloured solid lines refer to the GSMFs as a function of M app ⋆ . Using photometricallyequivalent stellar masses, the GSMF evolution is quite similar between H16F, IGIMF and PP11 runs (black dashed lines show the evolution in the H16F run). This is a very interesting conclusion, that highlights the intrinsic difficulty of using the stellar mass function as a discriminant of galaxy evolution, in a varying IMF context. However, the situation is dramatically different when considering the predicted GSMFs as a function of the intrinsic stellar mass M⋆ (dot-dashed coloured lines). As shown in F17a, in the IGIMF framework the differences are small for most redshifts, with the main deviations seen at low redshifts and high masses. In particular the GSMFs based on M app ⋆ are in good agreement with the z = 0 GSMF, thus easing the discrepancy between the intrinsic GSMF and observational data. The picture is completely different for the PP11 approach: the intrinsic GSMFs correspond to space densities systematically smaller than those obtained with M app star . The systematic difference holds for the highest redshift probed (i.e. z ∼ 7), and reduces at z < 1. This mismatch is a clear indication that the mass assembly is radically different in the PP11 run with respect to H16F (and IGIMF) runs. In the following sections we will focus on this aspect.

Comparison between PP11 and IGIMF theory.

The differences in the GSMF evolution can be better understood if we consider the distribution of GAEA model galaxies in the SFR -UCR/Umw space. Fig. 3 shows this distribution for galaxies from the PP11 run and highlights the differences between the two IMF scenarios. Contour plots mark the number density levels corresponding to 1, 5, 20 and 50 percent. Blue dots in Fig. 3 show the position of massive galaxies (i.e. M⋆ > 10 11 M⊙). In each panel, the red vertical and green horizontal lines correspond to the SFR and UCR/Umw values typical of a MW-like galaxy, i.e. to a canonical IMF in the PP11 and in the IGIMF theory, respectively. Above the green lines and to the right of red lines model galaxies are characterized by "Top-Heavy" IMFs in the corresponding variable IMF scenario, and viceversa. The figure show that most galaxies in the IGIMF run form stars with an IMF that is comparable and/or bottom-heavier than the MW one, while massive galaxies tend to have a top-heavy IMF, especially at high-redshifts. Deviations from a MW-like IMF are more marked in the PP11 model (i.e. there is a stronger evolution along the y-axis). In this case, the fraction of model galaxies forming stars with a MW-like (or bottom-heavy IMF) is decreasing with increasing redshift. A more quantitative comparison between the two models is complicated by the different evolution predicted for the IMF shape. Indeed, a "Top-Heavy" IMF in the PP11 runs is not the same as in the IGIMF framework.

To better compare the two scenarios, we use the fraction of SN per unit of stellar mass formed ( fSN). It should be noticed, however, that given the different shape of the IMF, similar values of fSN correspond to different ratios between SNIa and SNII. We compute fSN for the 6 IMF shapes considered in this work and for the 21 ones in F17a. The dashed magenta line in Fig. 3 connects the points corresponding to similar values of fSN in the two runs. If a model galaxy lies above this line, it is forming stars with an IMF that is top-heavier in the PP11 run than in IGIMF theory: this is the case for most sources at high-redshift. At z < ∼ 1 massive model galaxies are more evenly distributed around the magenta line, i.e. the two approaches predict comparable IMF shapes. Fig. 3 clearly shows that the PP11 implementation corresponds to larger deviations from the universal IMF hypothesis than in the IGIMF theory.

Mass assembly history.

Given the previous discussion, it is not surprising that the largest discrepancies between the two variable IMF models are seen in the mass assembly history of model galaxies. In order to highlight this effect, in Fig. 4 we show the mean star formation histories, cumulative mass assembly and evolution of the [O/Fe] ratio for four different z = 0 mass bins (corresponding to M⋆ ∼ 10 12 , 10 11.5 , 10 10.5 and 10 9.25 M⊙). Fig. 4 can be compared with a similar figure in F17a, where we showed that the star formation and mass assembly histories in the H16F and IGIMF runs are quite similar. The situation is completely different for the PP11 model (right column). The evolution of galaxies more massive than M⋆ ∼ 10 10.5 M⊙ is somehow accelerated, with a higher peak of star formation at high redshift and a shorter star formation time-scales with respect to their analogous in the IGIMF model (left column), thus enhancing the more rapid formation of more massive galaxies in hierarchical models of galaxy evolution (De Lucia et al. 2006). Model galaxies in the smallest mass bin show a reversal of this trend. While in the IGIMF (and H16F) runs they are characterized by flat star formation histories, in the PP11 run they clearly show a rising SFR at late times. This behaviour is driven by the implementation of the PP11 formalism, and only mildly affected by the recalibration of the model. In order to test this, we run a model realization (PPNT) using the PP11 scenario with the same GAEA parameters as in H16F: the corresponding results are shown in Fig. 4 (middle column) and clearly demonstrate that the change in mass assembly history for massive galaxies is driven by the IMF variation.

Mass excess distributions.

As we mentioned earlier, one of the key successes of runs implementing IMF variations (either assuming the IGIMF or the PP11 theory) is the prediction of a discrepancy between the intrinsic M⋆ and the apparent M app ⋆ , growing at increasing M⋆. If we interpret the former as a dynamical mass estimate and the latter as a photometric mass estimate, we can associate the M⋆/M app ⋆ ratio to the so-called "mass excess" found in dynamical and spectroscopic studies. At high-M⋆, the dispersion of the predicted relation is wide enough that a non-negligible fraction of model galaxies has M⋆/M app ⋆ ∼ 1, i.e. they do not show a mass excess. This result is interesting by itself, as it potentially explains contrasting results obtained for massive lensed galaxies (Smith et al. 2015;Leier et al. 2016), that are found to have mass-to-light ratios consistent with a universal IMF. In order to understand the origin of this galaxy population in our models, we consider the mass excess distribution for model galaxies in a slice corresponding to the mass range 11.4 < log(M⋆/M⊙) < 11.6 (Fig. 5). The position of the peak of the distribution is clearly different with respect to the H16F run (featuring a universal IMF) in both variable IMF scenarios. There are also differences between the PP11 and IGIMF runs, with the latc 0000 RAS, MNRAS 000, 000-000 ter showing a broader distribution (and more model galaxies compatible with the universal IMF scenario). We use the information in Fig. 5 to define two subsamples of model galaxies in the PP11 and IGIMF runs. The first one (subsample H -red shaded area) includes galaxies around the peak of the distribution, while the second one (subsample N -blue hatched region) corresponds to the tail of the distribution, i.e. galaxies with small mass excesses and whose global properties are compatible with the hypothesis of a universal IMF. These two subsamples are defined using different M⋆/M app ⋆ ratios, due to the different shape of the parent distribution in the two runs.

We then repeat the same analysis as in Fig. 4 for the galaxies in the two subsamples. In Fig. 6 we compare the resulting mean normalized star formation histories and cumulative mass assembly with those of the parent M⋆ ∼ 10 11.5 M⊙ population (black dotted lines). Subsample N galaxies (blue dashed lines) are characterized by a redshift evolution which is quite different with respect to both the total sample and subsample H (red solid line). In particular they . Mass excess distributions for M⋆ ∼ 10 11.5 M ⊙ galaxies in the PP11, IGIMF and H16F runs. The red shaded area (subsample H) includes galaxies around the peak of the distribution, while the blue hatched region (subsample N) selects galaxies with small mass excesses (i.e. compatible with the hypothesis of a universal IMF). . Parent halo mass distributions for M⋆ ∼ 10 11.5 M ⊙ galaxies in the PP11, IGIMF and H16F runs. Red shaded and blue hatched regions refer to subsample H and N, as in Fig. 5. Black, red and blue arrows mark the mean halo mass for the whole sample and subsample H and N, respectively tend to form the bulk of their population at later times, and over a longer star formation timescale with respect to the total population. This different evolutionary history relates to the different environment these galaxies live in. If we consider the distribution of their parent DM halo masses (Fig. 7), we clearly find a systematic difference between subsample N and H: galaxies with small (or no) mass excess tend to live in less massive haloes with respect to those close to the peak of the mass excess distribution. The black, blue and red arrows in Fig. 7 mark the mean halo mass associated with the global population, subsample H and subsample N, respectively. This environmental signature is particularly evident in the IGIMF run, where the two distributions are well separated, while in the PP11 run there is a relevant overlap between the two populations.

DISCUSSION & CONCLUSIONS

In this paper, we present a version of the GAEA semi-analytic model, including a model for a variable IMF based on PP11, who characterized the impact of the CR energy density on the properties of star forming regions. We assume that the main effect on the IMF shape is an evolution of the characteristic mass (i.e. the knee of the function), that increases with increasing SFR surface density (we use this quantity as a proxy for CR energy density). We contrast predictions from this model with those from our previous work, based on the IGIMF framework (Weidner & Kroupa 2005). This alternative scenario predicts a top-heavier IMF with increasing SFR. This represents a further development for GAEA, which is now able to handle two different prescriptions for IMF variations, driven by different physical properties of model galaxies (although SFRs and SFR surface densities are correlated quantities, as we show and discuss in Sec. 4.1). In this paper we discuss similarities c 0000 RAS, MNRAS 000, 000-000 and differences between the two IMF variation scenarios, and we discuss the implications in terms of galaxy evolution and assembly.

The PP11 and IGIMF formalisms predict different shapes for the IMF. The latter model is based on a deterministic derivation of the IMF, based on the assumption of a universal (Kroupa-like) IMF in individual clouds. The instantaneous galaxy-wide SFR is the main driver of the variation, and mainly affects the slope(s) of the IMF at the high-mass end. In the PP11 approach, the star formation conditions in dense UV-shielded MC cores are determined by the CR energy density. This sets the minimum temperature of the gaseous phase and modifies the Jeans mass estimate. The CR density can be connected to the star formation surface density, which we use as the physical parameter governing the IMF variation in this scenario. It is important to realise that in the PP11 model the IMF variation is a local process, taking place inside MCs, and possibly varying across different locations within galaxies. This is a clear difference with respect to the IGIMF theory, that offers an empirical derivation of the galaxy integrated IMF, based on and compatible with the observed properties of individual star clusters (Yan et al. 2017).

In this work, we model the IMF variation driven by CR regulation as a change in the knee of the IMF, keeping the high-mass and low-mass slopes fixed. Moreover, we compute the IMF variation in the entire star forming disc (i.e. we are using the total SFR surface density as the leading parameter). Further studies, with hydrodinamical simulations able to resolve the internal structure of galaxies, will allow us to understand the effect of these variable IMF scenarios on the radial properties of galactic discs.

The variable IMF models considered here provide consistent predictions for the evolution of key physical quantities, especially when these quantities are computed using an apparent -Chabrier IMF equivalent -stellar mass from synthetic photometry. In particular, the PP11 run is able to reproduce the [α/Fe]-stellar mass relation (see e.g. Johansson et al. 2012) and the mass-to-light (stellar mass) excess inferred for local elliptical galaxy samples. These results are in agreement with those found in our previous work, in the framework of the IGIMF model. Moreover, we show that these observational constraints are not able to disentangle between the two scenarios. Nonetheless, our results show the robustness of the predictions of galaxy formation models against different variable IMF frameworks.

In particular, in F17a, we showed that in the GAEA run implementing the IGIMF theory the mass assembly and star formation histories of model galaxies are hardly changed with respect to the H16F run adopting a universal IMF. The situation is different for the PP11 model. In this run, the intrinsic mass assembly of model galaxies is dramatically altered. Above 10 10.5 M⊙, mean star formation histories become more peaked, and the peak moves at higher redshift with respect to the H16F predictions. Albeit these differences are obvious and significant when considering model results, they are difficult to unveil from photometric data. Indeed, the stellar mass functions derived from photometrically derived M app ⋆ are indistinguishable between H16F, IGIMF and PP11 runs.

As in the case of the IGIMF model, PP11 predictions are in qualitative agreement with the results presented in Conroy et al. (2013), i.e. with an excess of the intrinsic stellar mass with respect to the photometric expectation. This mass excess is increasingly larger at increasing stellar mass. Our runs contain a significant number of model galaxies characterized by synthetic photometry and mass-to-light ratios compatible with a universal IMF scenario. This finding can explain results based on the lensing analysis of a limited samples of early-type galaxies (Smith et al. 2015;Leier et al. 2016). We study in detail the properties of this subsample, and we show that these model galaxies assemble on a longer timescale with respect to the bulk of the population, and preferentially live in less massive DM halos.

It is worth stressing, that our analysis is based on the mismatch between synthetic photometry and intrinsic physical properties of model galaxies. In particular, the mass excess in our runs is not due to an intrinsic "bottom-heavier" IMF in massive galaxies. This contrasts with the interpretation of Conroy & van Dokkum (2012) and La Barbera et al. (2013), based on the analysis of spectral features sensitive to the giants to low-mass stars ratio. In future work, we plan to deepen this aspect by generating synthetic spectra based on the star formation histories extracted from our variable IMF runs. We will then analyze these synthetic spectra using the same approach used in the observed samples, in order to better characterize the level of disagreement between theoretical predictions and reconstructed IMF shapes.

Another interesting approach is based on the study of the relative abundances of CNO isotopes in the ISM. Romano et al. (2017) showed that the cosmic evolution of isotopes of these chemical species is extremely sensitive to the fraction of massive stars (AGB stars and novae) and thus to the IMF shape. They argued that in order to explain the isotope ratios (mainly 12 C/ 13 C and 16 O/ 18 O) observed in starbursts, an IMF skewed towards high-stellar masses is required. Line ratios of isotopologues (i.e. molecules that differ in their isotopic composition) have practical advantages to study IMF variations: most of them are accessible in the submmillimetre regime and are mostly insensitive to dust obscuration. Current facilities like ALMA already extend the accessible redshift range up to z ∼ 3. All these attempts represent the necessary next steps in theoretical studies of IMF variations, and will be at the centre of our future work.

g-sdss @z=0 r-sdss @z=0 i-sdss @z=0 K-band @z=0 K-band @z~1.0 K-band @z~3.0 V-band @z=1.3 V-band @z=2.4 V-band @z=3.0

APPENDIX B: PHYSICAL PROPERTIES OF MODEL GALAXIES

In this Appendix we present the physical properties of model in the PP11 run, and contrast them against the H16F and IGIMF results. In all plots the green, red and black solid lines refer to the the PP11, IGIMF and H16F runs respectively. Unless differently stated, we always show predictions as a function of the apparent stellar mass (M app ⋆ ), computed from synthetic magnitudes using a mass-to-light vs colour relation as routinely done in observational samples (see e.g. Zibetti et al. 2009). In detail, we adopt a relation proposed by Zibetti et al. (in prep.):

log Υi = υ(g -i) + δ + ǫ (B1)
where Υi represents the stellar mass-to-light ratio in the i-band.

The best-fit coefficients υ = 0.9 and δ = 0.7 are derived using a Monte Carlo library of 500,000 synthetic stellar population spectra, based on the revised version of the Bruzual & Charlot (2003) simple stellar population tracks, and using the age-dependent dust attenuation prescription from Charlot & Fall (2000). Based on F17a . Bottom panel: evolution of the cold gas fraction in star-forming galaxies (i.e. SFR/M app ⋆ > 10 -2 Gyrs -1 ); left upper panel: z = 0 total stellar metallicity; right upper panel: z = 0 cold gas metallicity metallicity. Lines and colours as in Fig. A1. In all panels, grey points and areas represents observational constraints as in F17a (see references herein). results, we also include the additional factor ǫ = 0.13 in order to account for a shift in M app ⋆ with respect to M⋆, due to spatial resolution effects (see F17a for a complete discussion on the origin of this shift). This choice ensure that M app ⋆ and M⋆ are statistically equivalent in the H16F run (modulo some scatter). Eq. B1 provides a Chabrier (2003) IMF equivalent stellar mass; in the following, we neglect the small difference in normalization between the Kroupa (2001, which better describes the general form of the PP11 IMFs) and the Chabrier (2003) IMF.

Overall, we find similar trends in the physical properties of model galaxies predicted in the IGIMF and PP11 runs. The consistencies between the two approaches are particularly interesting, since they predict similar behaviour of the IMF in strongly star forming sources (either via high-SFR or high-ΣSFR), i.e. a topheavier IMF with respect to the local neighbourhood, but with a different shape evolution. In the IGIMF theory, the main evolution of the IMF is predicted at its high-mass end, while the low-mass end and the low-mass characteristic mass are unaffected. In PP11, both slopes are fixed, and the effect of a different UCR is mainly seen as a change in the characteristic mass.

However, some interesting differences are seen, which could be in principle used to discriminate between the two approaches. In Fig. B1, we consider the relation between [α/Fe] enrichment and stellar mass for B/T > 0.7 galaxies (which are considered a good proxy for elliptical galaxies). As in F17a, we will refer to [α/Fe] ratios for the observational data and to the [O/Fe] for theoretical predictions. This choice is related to the fact that oxygen represents the most abundant (so a good tracer) among α-elements, even if most observational estimate are actually derived from Magnesium lines. The predicted slope (green solid line -the hatched area represents the 1-σ scatter) of the [α/Fe]-stellar mass relation in PP11 is closer to the fits to the data by Thomas et al. (2010) and Johansson et al. (2012), while IGIMF theory predicts a somehow steeper relation, but still compatible with the data.

In Fig. B2, we present a collection of the additional relevant physical properties in galaxy evolution studies, namely the evolution of the cold gas fractions of star-forming galaxies (SFR/M app ⋆ > 10 -2 Gyr -1 -bottom panels) and z = 0 massmetallicity relations (total metallicity -upper left panel -and cold gas metallicity -upper right panel). Cold gas fractions (and their redshift evolution) appear to be the direct observable most sensitive to IMF variation. Nonetheless, the present uncertainties on observed samples prevent a clear separation between models: PP11 better reproduce the observational data, but the IGIMF run is still within the 1-σ scatter. It is also worth pointing out the peculiar behaviour of PP11 predictions on the cold gas mass-metallicity relation (Fig. B2 -upper right panel). The model predict a clear turnover at ∼ 10 11 M⊙, which is slightly outside the range covered by the data. It is however unclear if this feature could be used as a test for a variable IMF model, since we test that this is particularly sensible to the strength of AGN feedback and its interplay with stellar feedback schemes.

In Fig. B3 we show the ratio between photometrically derived and intrinsic quantities as a function of intrinsic quantities for B/T > 0.7. Those plots are meant to be compared with the corresponding observational data suggesting deviations from the assumption of a universal (Chabrier-like) IMF. The left panel considers the ratio of proper stellar mass-to-light ratios in the i-band (M⋆/Li) and the photometric equivalent Υi from Eq. B1: this plot has to be compared to the dynamical analysis from Cappellari et al. 2012 in the ATLAS 3D sample of early-type galaxies. The right panel shows the evolution of the M⋆/M app ⋆ ratio as a function of the intrinsic stellar mass and roughly corresponds to the spectroscopic analysis as in Conroy et al. (2013). Interesting differences between the IGIMF and PP11 runs are also evident in these plots and are particularly relevant for the dynamical estimate (Fig. B3, right panel), with the PP11 runs showing a clear steeping of the relation with respect to the IGIMF run. The steeping is particularly c 0000 RAS, MNRAS 000, 000-000 relevant at low-mass-to-light ratios, where the PP11 model predicts a deficit in the photometric with respect to the dynamical estimate.

Figure 1 .

 1 Figure1. Variable IMF scenarios. Left panel: evolution of the IMF shape according to PP11. Labels close to lines report the corresponding cosmic ray densities (normalized to the Milky Way value, CR = U CR /U MW ). Right panel: IMF variation as predicted by the IGIMF theory. Labels close to the lines report the reference SFRs (see alsoFig. 1 in F17a). In both panels, each IMF is normalized to 1 M ⊙ in the stellar mass interval 0.1-100 M ⊙ .

Figure 2 .

 2 Figure 2. Redshift evolution of the galaxy stellar mass function. In each panel, the black long-dashed line refers to the predictions of the H16F model; the solid lines to the GSMF as a function of M app ⋆ and the dotdashed lines to the GSMF as a function of the intrinsic galaxy stellar mass M⋆. Grey points show a collection of observational measurements as in Fontanot et al. (2009, see detailed references herein).

Figure 3 .

 3 Figure3. Distribution of model galaxies in the SFR vs U CR /U MW space in the PP11 run, at different redshift. Contours correspond to number densities for the global population, while blue dots mark the position of Mstar > 10 11 M ⊙ galaxies. Red vertical and green horizontal lines mark the values corresponding to a MW-like IMF in IGIMF and PP11 formalism, respectively. Magenta dashed line connects points corresponding to IMFs with similar f SN in the two formalisms (see text for more details).

Figure 4 .

 4 Figure4. Mean normalized star formation histories (upper panels), cumulative mass assembly (middle panels) and evolution of the [O/Fe] ratio (lower panels) for galaxies in the different logarithmic stellar mass intervals (see caption). The predictions relative to the IGIMF and PP11 models are shown in the righthand and left-hand columns, respectively. The middle column shows the prediction for a model (PPNT) implementing the PP11 model, but using the same parameters as in H16F.

  Figure5. Mass excess distributions for M⋆ ∼ 10 11.5 M ⊙ galaxies in the PP11, IGIMF and H16F runs. The red shaded area (subsample H) includes galaxies around the peak of the distribution, while the blue hatched region (subsample N) selects galaxies with small mass excesses (i.e. compatible with the hypothesis of a universal IMF).

Figure 6 .

 6 Figure6. Mean normalized star formation histories (top panels) and cumulative mass assembly (bottom panels) for M⋆ ∼ 10 11.5 M ⊙ galaxies. Black dotted, blue dashed and red solid lines refer to the total sample, to subsample N and subsample H, respectively.

Figure A1 .

 A1 Figure A1. Model calibration set. Predicted luminosity functions in different wavebands and at different redshifts. Solid Black, red and green lines refer to the H16F, IGIMF and PP11 runs respectively. Grey points show the compilation of observational estimates in the SDSS g, r and i-band, K-and V-band from F17a (see detailed references herein).

Figure B1 .

 B1 Figure B1. Predicted [O/Fe] ratios in B/T > 0.7 model galaxies compared to the observed [α/Fe] ratios in local elliptical galaxies. Data collection (grey points, contours, dot-dashed and long-short dashed lines) as in F17a (see detailed references herein).

Figure B2 .

 B2 Figure B2. Evolution of key galaxy properties as a function of photometrically estimated M app ⋆

Figure B3 .

 B3 Figure B3. Left panel: ratio of proper M⋆/L i and the photometric equivalent Υ i , as a function of M⋆/L i . Right panel: M⋆/M app ⋆ versus M⋆. In each panel, only B/T > 0.7 galaxies have been considered and colours and linetypes corresponds as in Fig. A1. Grey dashed contours mark galaxy number densities (in percentage of the maximum density) in the PP11 run.

Table 2 .

 2 Parameter values adopted for the runs considered in this study (parameters are defined in Appendix B)

		Parameter	H16F	IGIMF PP11
		α SF	0.03	0.19	0.08
	.	ǫ reheat	0.3	0.575	0.30
		ǫ eject	0.1	0.12	0.04
		γ reinc	1.0	1.0	0.50
		κ radio /10 -5	1.0	1.78	1.18

c 0000 RAS, MNRAS 000, 000-000

We explicitly tested that results from our reference PP11 run are qualitatively similar to those of a run implementing the same variable IMF model, but without any f SN scaling.

As shown in previous papers (e.g.Wang et al. 2008;Guo et al. 2013), we do not expect the mismatch of these cosmological parameters with respect to the most recent measurements (Planck Collaboration XVI 2014) to change our main conclusions. c 0000 RAS, MNRAS 000, 000-000

ACKNOWLEDGEMENTS

We thank P.P. Papadopoulos for enlightening discussions on the details of the CR-regulated model for the IMF. FF and GDL acknowledge financial support from the grants PRIN INAF 2014 "Glittering kaleidoscopes in the sky: the multifaceted nature and role of Galaxy Clusters" (1.05.01.94.02). GDL acknowledges financial support from the MERAC foundation. MH acknowledges financial support from the European Research Council via an Advanced Grant under grant agreement no. 321323 (NEOGAL). GB acknowledges support for this work from the National Autonomous University of México (UNAM), through grant PAPIIT IG100115.

APPENDIX A: MODEL CALIBRATION

The original H16F model is tuned against the reconstructed evolution of physical quantities such as M⋆ and SF R. However, such constraints are usually derived under the hypothesis of a universal IMF and therefore cannot be used in a variable IMF scenario. As already discussed in F17a, our model, based on the PP11 approach, requires a recalibration of key model parameters, given the different impact on galaxy evolution of the assumed IMFs (in terms of baryons locked in long-living stars, fraction of SNe and ratio between SNIa and SNII). In detail, the relevant parameters are the SFR efficiency (αSF), AGN feedback (κ radio ), stellar feedback reheating (ǫ reheat ) and ejection rate (ǫeject), and the reincorporation rate (γreinc). We refer the interested reader to Hirschmann et al. (2016) for a more detailed discussion of the role these parameters play in the GAEA context. The assumed calibration set is the same defined in F17a (see there for a full reference list), i.e. the evolution of the Kand V -band LF at z < ∼ 3, and the z = 0 LFs in the Sloan Digital Sky Survey (SDSS) g, r and i-bands. Fig. A1 compares the resulting luminosity functions (green lines) with those obtained in the context of the IGIMF theory (red lines -High-αSF model in F17a) and using the reference H16F model (black lines). c 0000 RAS, MNRAS 000, 000-000