Comparing Models for IMF Variation Across Cosmological Time in Milky Way-like Galaxies

TL;DR

This study compares various models predicting the initial mass function (IMF) variation over cosmic time in Milky Way-like galaxies using high-resolution cosmological simulations, revealing that most models predict more variation than observed.

Contribution

It integrates cloud-scale IMF models with cosmological simulations to evaluate their predictions across galaxy evolution, a novel approach in IMF studies.

Findings

Most models predict greater IMF variation than observed in the Milky Way.

Feedback-regulated models do not show excessive IMF variation.

Models based on Jeans mass, turbulence, or fragmentation predict 0.6-1 dex scatter.

Abstract

One of the key observations regarding the stellar initial mass function (IMF) is its near-universality in the Milky Way (MW), which provides a powerful way to constrain different star formation models that predict the IMF. However, those models are almost universally "cloud-scale" or smaller -- they take as input or simulate single molecular clouds (GMCs), clumps, or cores, and predict the resulting IMF as a function of the cloud properties. Without a model for the progenitor properties of all clouds which formed the stars at different locations in the MW (including ancient stellar populations formed in high-redshift, likely gas-rich dwarf progenitor galaxies that looked little like the Galaxy today), the predictions cannot be fully explored, nor safely applied to "live" cosmological calculations of the IMF in different galaxies at different cosmological times. We therefore combine a suite of high-resolution cosmological simulations (from the Feedback In Realistic Environments project), which form MW-like galaxies with reasonable star formation properties and explicitly resolve massive GMCs, with various proposed cloud-scale IMF models. We apply the models independently to every star particle formed in the simulations to synthesize the predicted IMF in the present-day galaxy. We explore models where the IMF depends on Jeans mass, sonic or turbulent Bonner-Ebert mass, fragmentation with a polytropic equation-of-state, or where it is self-regulated by protostellar feedback. We show that all of these models, except the feedback-regulated ones, predict far more variation ($\sim 0.6-1$ dex $1\,\sigma$ scatter in the IMF turnover mass) in the simulations than is observed in the MW.

Tables

Table 1. Table 1: Rows: Different IMF models compared in this paper, each with the predicted scaling of the IMF turnover mass M crit subscript 𝑀 crit M_{\rm crit} with initial parent cloud properties (§  2.2 ), reference, and the predicted 1 ​ σ 1 𝜎 1\sigma dispersion in log 10 ⁡ ( M crit ) subscript 10 subscript 𝑀 crit \log_{10}(M_{\rm crit}) across the galaxy at present-day (averaging Fig.  8 over all galacto-centric radii). We measure σ M crit subscript 𝜎 subscript 𝑀 crit \sigma_{M_{\rm crit}} from five simulations: galaxies m12i , m12f , and m12m are three distinct Milky Way-mass ( ∼ 10 12 ​ M ⊙ similar-to absent superscript 10 12 subscript 𝑀 direct-product \sim 10^{12}\,M_{\odot} ) halos which produce similar disky, Milky Way-like galaxies (stellar mass ∼ 0.5 − 1 × 10 11 ​ M ⊙ similar-to absent 0.5 1 superscript 10 11 subscript 𝑀 direct-product \sim 0.5-1\times 10^{11}\,M_{\odot} ), but have different formation histories (see Fig. 2 for visualizations and Hopkins et al. 2017 for details). For each we label the mass resolution (in M ⊙ subscript 𝑀 direct-product M_{\odot} ) and minimum density n crit subscript 𝑛 crit n_{\rm crit} for creation of stellar sink particles. For m12i , we compare two alternative runs: one at lower resolution, and one including magnetic fields ( m12i+MHD ). The predicted IMF variation is remarkably robust across all these simulations.

Model	$M_{\mathrm{crit}}$	Reference	Galactic IMF variation (${\sigma }_{M_{\mathrm{crit}}}$) [dex]
			m12i	m12i	m12i+MHD	m12f	m12m
			($56000M_{☉}$)	($7000M_{☉}$)	($7000M_{☉}$)	($7000M_{☉}$)	($7000M_{☉}$)
			($1000{\mathrm{cm}}^{-3}$)	($1000{\mathrm{cm}}^{-3}$)	($1000{\mathrm{cm}}^{-3}$)	($1000{\mathrm{cm}}^{-3}$)	($1000{\mathrm{cm}}^{-3}$)
Jeans Mass	$\propto T^{3/2}{\rho }^{-1/2}$	Bate & Bonnell 2005	0.60	0.54	0.61	0.56	0.65
Turbulent/Sonic Mass	$\propto T{R}_{\mathrm{sonic}}$	Hopkins 2012	0.91	0.96	0.86	0.91	0.84
Opacity-limited, ${\rho }_{\mathrm{crit}}$	$\propto T^{3/2}$	Jappsen et al. 2005	0.63	0.53	0.58	0.54	0.61
Opacity-limited, ${\Sigma }_{\mathrm{crit}}$	$\propto T^2$	Bate 2009	0.81	0.70	0.75	0.73	0.81
Protostellar Heating	$\propto {\left(\rho T\right)}^{-1/18}$	Krumholz 2011	0.030	0.026	0.031	0.027	0.031

Equations

$$M_{\mathrm{crit},\mathrm{J}}\sim \frac{\pi c_s^3}{6G^{3/2}{\rho }^{1/2}}.$$

$$M_{\mathrm{crit},\rho }\sim \frac{\pi c_s^3}{6G^{3/2}{\rho }_{\mathrm{crit}}^{1/2}},$$

$$M_{\mathrm{crit},\mathrm{S}}=M_{\mathrm{sonic}}\sim \frac{2c_s^2{R}_{\mathrm{sonic}}}{G},$$

$${\sigma }_{\mathrm{turb}}^2(\lambda )=c_s^2\frac{\lambda }{R_{\mathrm{sonic}}}.$$

$$M_{\mathrm{crit},\mathrm{B}}\sim 0.5{\left(\frac{\rho }{1.2\times 10^{-19}g/c{\mathrm{m}}^3}\right)}^{-1/5}{\left(\frac{L_{\ast }}{150{\mathrm{L}}_{\text{sun}}}\right)}^{3/10}{\mathrm{M}}_{\text{sun}},$$

$$M_{\mathrm{crit},\mathrm{K}}\sim 0.15{\left(\frac{P/k_B}{10^{6}K/c{\mathrm{m}}^3}\right)}^{-1/18}{\mathrm{M}}_{\text{sun}}$$

Full text

1 Introduction

The (instantaneous) mass distribution of stars at their formation time, also known as the initial mass function (IMF), is one of the key predictions of any star formation model. This governs essentially all observable and theoretical aspects of star formation and stellar populations – observable luminosities and colours; effects on stellar environments via feedback in the form of stellar winds, radiation, supernovae; nucleosynthesis and galactic chemical evolution, and so on. The IMF has been well-studied within the MW, and appears to be well-fit by a simple function with a Salpeter (1955) power-law slope at high masses and lognormal-like turnover at low masses (Chabrier, 2005; Kroupa, 2002). Perhaps the most interesting feature of the IMF is its universality: it has been found that there is quite weak variation within the MW (for recent reviews, see Chabrier, 2003; Bastian et al., 2010; Offner et al., 2014; Krumholz, 2014, and references therein), albeit with a few possible outliers (e.g. Luhman et al., 2009; Kraus et al., 2017). As Offner et al. (2014) emphasize, this universality includes both very young ($\sim$ Myr-old) and very old ($\sim 10\,$Gyr-old) stellar populations; stars forming in small, nearby GMCs with masses $\sim 10^{4}-10^{6}\,M_{\odot}$ and massive complexes with masses $\sim 10^{6}-10^{7}\,M_{\odot}$; the solar neighbourhood at $\sim 10$ kpc from the galactic centre (where the gas disk surface density is $\sim 10\,M_{\odot}\,{\rm pc^{-2}}$) and the central molecular zone at sub-kpc and $\sim 100\,$pc scales (where gas surface densities are order-of-magnitude larger).

In other galaxies, the IMF usually must be assumed, and with an IMF assumption, physical properties of the stellar populations and galaxies (e.g. their stellar masses) are derived from observables (e.g. light, colours). This makes it critical to understand the IMF, in order to understand galaxy formation. Likewise it is critical for models of galaxy formation to predict or assume some IMF model, in order to make any meaningful predictions for observable quantities. The universality of the IMF in old stellar populations in the MW is widely taken as a suggestion that it may be near-universal in other galaxies, because older populations in the MW formed when the galaxy was much younger and very different, likely a typical high-redshift gas-rich, metal-poor dwarf galaxy. There are indirect constraints on the IMF both from spectral features and integrated mass-to-light constraints in nearby galaxies: these mostly also favor a universal IMF (e.g. Fumagalli et al., 2011; Koda et al., 2012; Andrews et al., 2013, 2014; Weisz et al., 2015). More recently there have been more interesting hints of variation in the centres of massive elliptical galaxies (van Dokkum & Conroy, 2010, 2011; Conroy & van Dokkum, 2012; Treu et al., 2010; Sonnenfeld et al., 2015; Cappellari et al., 2012; Posacki et al., 2015; Martín-Navarro et al., 2015b, a, c), and perhaps also in faint dwarf galaxies (Hoversten & Glazebrook, 2008; Brown et al., 2012; Geha et al., 2013). Even so, it is worth stressing that the implied variation is not radical: it implies variation of a factor $<2$ in the stellar mass-to-light ratio.

As a result, there is a long history of both theoretical and empirical models for galaxy formation which have attempted to predict the IMFs that should arise in different galaxy populations, as a function of either galaxy-scale or $\sim$kpc-scale properties (what can be resolved in most previous calculations) within the galaxies (see e.g. Baugh et al., 2005; van Dokkum, 2008; Hopkins, 2013; Davé et al., 2012; Narayanan & Davé, 2012; Recchi & Kroupa, 2015; Lacey et al., 2016; Blancato et al., 2016). At the same time, the specific hints of galaxy-to-galaxy variation discussed above have prompted a new wave of theoretical models which argue the IMF could vary under certain conditions at the GMC or sub-GMC scale, in a way that may connect to the systematic variation inferred in different galaxies (e.g. Weidner et al. 2013; Bekki 2013; Chabrier et al. 2014; Ferreras et al. 2015).

However, these models in every case rely on very strong simplifying assumptions – the IMF is predicted as a function of the cloud properties out of which the stars form, such as its temperature, density, turbulent velocity dispersion, virial parameter, etc. (from which properties like the Jeans mass, or the turbulent Bonner-Ebert mass, or the IGIMF turnover mass, are determined). It is impossible at present to know these empirically because all the clouds that formed the old stellar populations in a galaxy (or even for most clouds within a galaxy at present day), so instead some strong additional assumptions are usually applied. For example, authors assume isothermal gas with $T=10\,$K (or some other temperature) at all densities, a universal linewidth-size relation across all galaxies, redshifts, and regions within galaxies, or a Jeans or sonic mass within clouds that somehow scales proportionally to that measured from the gas at the $\gtrsim$kpc scales resolved in the cosmological calculations. But if these properties vary across cosmic time, or cloud-to-cloud, then any such model will produce variation in the predicted IMF which can be compared to the observational limits within the MW.

In this paper we therefore investigate the predicted variation in the IMF peak imprinted by these physics in a number of IMF models. We combine high-resolution simulations of MW-like galaxies (where the cloud-scale properties can be at least partially resolved) with the relevant small-scale models for IMF variation as a function of cloud properties.

2 Model and Methodology

2.1 Simulation

We utilize a set of numerical simulations of MW-like galaxies (see Table 1) presented in Hopkins et al. (2017), from the Feedback in Realistic Environments (FIRE) project (Hopkins et al. 2014).111http://fire.northwestern.edu The simulations are fully cosmological “zoom-in” simulations (where resolution is concentrated on one galaxy in a large cosmological box, run from redshift $z>100$ to today) and are run using GIZMO (Hopkins 2015)222http://www.tapir.caltech.edu/~phopkins/Site/GIZMO.html, with the mesh-free Godunov “MFM” method for the hydrodynamics (Hopkins, 2015). Self-gravity is included with fully-adaptive force and hydrodynamic resolution; the simulation mass resolution is fixed at $7000$ or $56000\,M_{\odot}$ (Table 1). The simulations include detailed metallicity-dependent cooling physics from $T=10-10^{10}\,$K, including photo-ionization/recombination, thermal bremsstrahlung, Compton, photoelectric, metal line (following Wiersma et al. 2009), molecular, fine structure (following Ferland et al. 2013), dust collisional and cosmic ray processes, including both a meta-galactic UV background and each star in the simulation as a local source.

Individual stars are not resolved in the simulations; but star formation is approximated from resolved scales via a sink-particle method. Gas which is locally self-gravitating, self-shielding, Jeans unstable, and exceeds a minimum density $n>n_{\rm crit}=1000\,{\rm cm^{-3}}$ (Table 1) is transformed into “star cluster sink particles” on its dynamical time. Each such particle represents an IMF-averaged single stellar population of the same age and metallicity, with mass equal to the mass resolution.

Once formed, the simulations include feedback from these star particles via OB & AGB mass-loss, SNe Ia & II, and multi-wavelength photo-heating and radiation pressure; with inputs taken directly from stellar evolution models (Leitherer et al., 1999), assuming (in-code) a universal IMF (Kroupa, 2002).

There are two reasons for using cosmological simulations instead of present-day observational data or a more localized cloud simulation. (1) We wish to test and validate the approach of using cloud-scale IMF models dynamically in next-generation cosmological simulations. Because stellar feedback and observable properties depend on the IMF, truly self-consistent predictions should include some sub-grid IMF model. These cosmological simulations were run assuming a universal IMF, but others (see references in § 1) adopt a dynamical IMF model based on resolution-scale properties. But it has not been asked whether the models they considered violate observational constraints in the MW. (2) Stars at a given present-day position in a galaxy can form at wildly different times/places (some even in other dwarf galaxies). This is especially true for the stars in old MW clusters, which appear to have formed at high redshifts, probably in distinct dwarf galaxies with entirely distinct radiation fields, turbulent velocity dispersions, gas masses, etc. It is impossible to know the distribution of progenitor cloud properties at star formation for old stars (needed for a given IMF model to make predictions) from observations (or localized simulations) alone.

The galaxies studied here, shown in Figs. 1-2 – have been studied extensively in previous work: they are similar to the MW in their stellar mass, present-day gas fractions and SFRs (Hopkins et al., 2014), and metallicity (Ma et al., 2016b). Our “fiducial case” m12i is also similar to the MW in its stellar kinematics, thin+thick disk morphology, metallicity gradient and metal abundance ratio gradients (in both vertical and radial directions) stellar age distribution (Ma et al., 2016a, 2017), $R$-process element distribution (van de Voort et al., 2015) and galactic stellar halo and satellite dwarf population (Wetzel et al., 2016). The other two examples represent a slightly later-type (m12f) and earlier-type (m12m) galaxy, at the same stellar mass and SFR. This is particularly useful because of course no simulation will exactly match the formation history of the MW, so it is important to understand whether our predictions are sensitive to this.

These and other FIRE simulations have also been shown to reproduce the observed Kennicutt-Schmidt relation (Sparre et al., 2015; Orr et al., 2017), properties of galactic outflows (Muratov et al., 2015) and (in higher resolution, non-cosmological simulations) the observed mass function (and CO luminosities), size-mass, and linewidth-size distributions of GMCs (Hopkins et al., 2012; Hopkins et al., 2013). One might reasonably worry that this cannot be captured at the lower resolution necessary in cosmological simulations. Therefore in Fig. 4 we plot the mass function and linewidth-size relation of GMCs identified at present-day in the actual simulations studied here. They appear to agree at least plausibly with observed properties (Dobbs et al., 2014; Heyer & Dame, 2015). Note that our mass resolution introduces a cutoff at the low mass end of the GMC mass function because these clouds can not be resolved by the simulation. However, all simulations included in this paper do resolve the most massive GMCs ($>10^{6}\,\mathrm{M}_{\rm\sun}$), in which most of the mass is concentrated (owing to the shape of the GMC mass function), allowing us to recover galactic properties even at lower resolutions. This is clearly apparent in the linewidth-size relation which shows a good agreement with Bolatto et al. (2008). All of this is not to say that the simulations are perfect analogues to the MW; however they are at least a reasonable starting point (see Fig. 4 for details).

Our choice of $n_{\rm crit}=1000\,{\rm cm^{-3}}$ (at mass resolution $7000-56000\,M_{\odot}$) can be justified by assuming that GMCs follow the mass-size relation of Bolatto et al. (2008) ($M_{\rm cloud}\sim\pi\,(85\,M_{\odot}\,{\rm pc^{-2}})\,R_{\rm cloud}^{2}$), which implies the density threshold for star formation is slightly higher than the mean density of the most-dense resolved clouds. More specifically, we chose the density threshold to correspond to the typical density where the Jeans/Toomre fragmentation scale falls below our mass resolution. In either case the GMC mass function and SFR is dominated by the most massive (hence well-resolved) clouds. This is evident in Fig. 4, where we show the GMC mass function and linewidth-size relation predicted at present-day in the galaxy both (a) agree reasonably well with observations (within a factor $\sim 2$ at all cloud sizes/masses resolved), and (b) are insensitive to resolution (except, of course, that at higher resolution they extend to smaller GMCs). This gives us some confidence that our predictions are not strongly resolution-dependent. In Table 1 we show that varying resolution and physics (in an otherwise identical run including magnetic fields, m12i+MHD) do not significantly alter our predictions.

Because the simulations resolve down to cloud scales, but no further, we treat each star-forming gas element as an independent “parent cloud”, which sets the initial conditions for its own detailed IMF model (in accordance with the IMF models we investigate). Specifically, whenever a sink particle is spawned, we record all properties of the parent gas element from which it formed, and use these in post-processing to predict the IMF. Fig. 5 shows the properties of gas elements at one instant, $z=0$, weighted by star formation rate. Integrating over all times and all galaxies which form stars that ultimately reside in the final galaxy, Fig 6 shows the density and temperature distribution of these “star forming particles” (gas at the moment the simulation assigned its mass to a sink particle) at the time of their formation, from our high-resolution m12i run. Not surprisingly most sinks form around the simulation density threshold from this particular run ($\sim 1000\,{\rm cm^{-3}}$). This choice has no effect on the scatter in both $n_{\rm cloud}$ and $T_{\rm cloud}$, which are the relevant parts to our study. Scatter in these quantities translate to a scatter in the local IMF according to the IMF models we are studying. Note that Narayanan & Hopkins (2013) show the inferred temperature range from mock CO-ladder observations will tend to be significantly smaller than the range plotted here. We wish to emphasize that what is plotted in Fig. 6 is not the density/temperature of the core or proto-stellar gas that which directly collapses and forms stars; that is not resolved in these simulations. Instead these are the properties of the progenitor molecular clouds, measured at the smallest resolved scales, which will (and must, physically) fragment into denser sub-clumps that can directly form stars. Also, the width in this distribution is expected to be higher than in present day clouds because of the longer lifetime of stars which preserves the progenitor cloud properties in their IMF for cosmological timescales.

2.2 From Parent Cloud to IMF Properties

From this point we infer the IMF turnover mass from the initial conditions of these parent clouds. This exercise has been done in detail by Guszejnov et al. (2016) where the semi analytical framework of Guszejnov & Hopkins (2016) was utilized to create a mapping between GMC properties and the IMF. Fig. 7 shows how the IMF peak scales with initial temperature in an equation of state (EOS) and a protostellar feedback based IMF model. Such scaling relations can be analytically derived for other IMF models (e.g. Jeans mass) as well – we focus here on how each model predicts the turnover or “critical” mass $M_{\rm crit}$ scale, because this is the most identifiable feature of the IMF (it sets the mass-to-light ratio, and varies significantly between models). In contrast the bright-end slope varies negligibly between models333Note that observations do indicate variations in the IMF slopes in extra-galactic populations (e.g. Cappellari et al., 2012; Spiniello et al., 2012; Shu et al., 2015) but these measurements only sample the relatively low mass region of the IMF ($<\mathrm{M}_{\rm\sun}$). The IMF in the MW, however, is well sampled at higher masses and appears to be consistent with a near-universal power-law tail (Offner et al., 2014), which most IMF models are able to roughly reproduce (see references in Table 1 for specifics in each case), so it is not useful as a diagnostic.

In this paper we investigate the sensitivity to initial conditions for the following classes of IMF models (summarized in Table 1):

• •

Jeans mass models: The Jeans instability is the primary mechanism for the collapse of gas clouds into stars, so these models assume that IMF properties are set by local mean Jeans mass of the parent molecular cloud complex (e.g. Bate & Bonnell 2005). Therefore, the critical mass is

[TABLE]

Note that the models may still assume sub-fragmentation to smaller scales, but the key assumption (for our purposes) is simply that the turnover mass somehow scales proportional to the parent cloud Jeans mass.

• •

Opacity limit equation of state (EOS) models: As the molecular gas becomes denser it reaches the point where it becomes opaque to its own cooling radiation, leading to a transition from isothermal to adiabatic behavior, terminating fragmentation at the Jeans mass at this density. This can occur at a critical volume density $\rho_{\rm crit}$ (e.g. Low & Lynden-Bell 1976; Whitworth et al. 1998; Larson 2005; Glover & Mac Low 2007; Jappsen et al. 2005; Masunaga & Inutsuka 2000). Motivated by radiation transfer simulations like Bate 2009 we also investigated the case where the transition occurs at a critical surface density $\Sigma_{\rm crit}$. The critical masses in these cases are:

[TABLE]

where $\rho_{\rm crit}$ and $\Sigma_{\rm crit}$ are the critical densities for the isothermal-adiabatic transition.

• •

Turbulent/sonic mass models: A number of analytical theories derive the CMF and IMF from the properties of the turbulent medium, in which they form (e.g. Padoan & Nordlund 2002; Hennebelle & Chabrier 2008; Hopkins 2012; Hennebelle & Chabrier 2013). In these models, both the CMF and IMF peaks are set by the “sonic mass” $M_{\rm sonic}$, namely the turbulent Jeans or Bonner-Ebert mass at the sonic scale ($R_{\rm sonic}$) below-which the turbulence becomes sub-sonic and therefore fails to generate large fluctuations (which seed fragmentation). The critical mass is

[TABLE]

where $R_{\rm sonic}$ is defined through the linewidth-size relation

[TABLE]

In our calculations $\sigma_{\rm turb}^{2}$ is estimated from the simulations when a star particle forms by measuring the velocity dispersion (after subtracting the mean shear) between neighboring particles in a sphere of radius $\lambda$ (taken to be that which encloses the nearest $\sim 32$ gas neighbours).

• •

Protostellar feedback models: Although there are a number of ways newly-formed stars can regulate star formation, most studies have concluded that at the scale of the IMF peak (early protostellar collapse of $\sim 0.1\,M_{\odot}$ clouds) the most important self-regulation mechanism is radiative feedback from protostellar accretion (Bate, 2009; Krumholz, 2011; Guszejnov et al., 2016). This sets a unique mass and spatial scale within which the protostellar heating has raised the temperature to make the core Jeans-stable, terminating fragmentation. The resulting critical masses are

[TABLE]

where $L_{*}$ is the average luminosity of accreting protostars and $P$ is the pressure. These different formulas come from Bate (2009) and Krumholz (2011), respectively; the differences are due to the detailed uncertainties treating radiation. However for our purposes they give nearly identical results, so we will focus on the model from Krumholz (2011).

3 Results and Discussion

Fig. 8 shows that there is significant variation the properties of the progenitor GMC complexes which formed stars that ultimately end up at a specific galacto-centric radius. We stress that this is not the variation of properties in present-day star-forming clouds, but includes all variations in time as well: if the galaxy progenitor was gas-rich (gas fraction $\sim 1/2$) at $z\sim 1-2$ for example, then the midplane gravitational pressure ($\sim G\,\Sigma_{\rm gas}^{2}$) would have been a factor $\sim 100$ larger than in the galaxy today. Fig. 9 shows that this, in turn, produces large IMF variations in all models here except those accounting for protostellar heating. Such variations ($>0.5\,$dex in $M_{\rm turnover}$) are strongly ruled-out by observations (Bastian et al., 2010). Note that these results are robust to variations in simulation parameters (see Table 1).

The variations in the IMF predicted by some of the simple models here (e.g. the Jeans-mass models) have often been substantially underestimated in previous work in the literature. In analytic models of the IMF (see references in Sec. 1) or galaxy-scale models which fail to resolve individual “parent clouds”, but post-process the entire galaxy (with $>$ kpc-scale resolution) to determine an IMF (Narayanan & Davé, 2012; Hopkins, 2013; Blancato et al., 2016), it is commonly assumed that all star-forming clouds are uniformly at the same isothermal temperature (e.g. $T=10\,$K at all densities), virial parameter, and lie exactly on the same linewidth-size relation. For example, if the gas at all densities and all cosmic times had exactly the same temperature, then the variation in the IMF for the opacity-limited EOS models would vanish (all clouds and cores lie on exactly one adiabat). This assumption is not correct, however, as even in the present-day MW (e.g. fixed redshift and galaxy properties) both GMC and clump temperatures (e.g. Bergin & Tafalla, 2007; Mills & Morris, 2013; Sánchez-Monge et al., 2014; Ott et al., 2014; Nishimura et al., 2015)444Note that Nishimura et al. (2015) only focuses on Orion A and B so these results are not necessarily representative of the entire MW. Also, as our stars form primarily around $n_{\rm crit}=1000\,{\rm cm}^{-3}$, the average temperature of the progenitor clouds is higher than the observed GMCs because they have not reached the cooler, higher-density fully-molecular phases (see Figs. 5-6). and virial parameters (e.g. Kauffmann et al., 2013; Svoboda et al., 2016) vary substantially. As expected variations are more pronounced in other nearby dwarf or star-forming galaxies (Gorski et al., 2017; Tang et al., 2017) or redshift $z\gtrsim 1-2$ galaxies and starburst systems (see e.g. Ott et al., 2011; González-Alfonso et al., 2012; Narayanan & Krumholz, 2014; Mangum et al., 2013; Miyamoto et al., 2015; Zschaechner et al., 2016), which are better analogues to the progenitors where many of the stars in the present-day Galaxy formed.

It is certainly possible that we (and these observations) have over-estimated the range of temperatures of GMCs in different environments. But the strong temperature sensitivity of the EOS models (e.g. $\propto T^{2}$) means that the temperature of all progenitor clouds, at all redshifts and in all progenitor galaxies, which formed stars in today’s MW, would have to lie within a scatter of just $\sim 20\%$ in temperature (smaller than that observed in just solar-neighborhood clouds) in order to avoid exceeding the allowed IMF variation in the MW.

Moreover the linewidth-size relation is observed to vary systematically, both within the MW and galaxy-to-galaxy, with high-redshift galaxies (the progenitors of the MW) differing by more than an order of magnitude (see e.g. Swinbank et al., 2011, 2015; Canameras et al., 2017, and references therein). Even if temperature variations are neglected entirely, in the “Turbulent/Sonic Mass” models the turnover mass is proportional to the square of the deviation ($(\sigma_{\rm cloud}/\langle\sigma[R]\rangle)^{2}$) of each cloud from the linewidth-size relation (Hopkins, 2012; Hennebelle & Chabrier, 2008), but these deviations are observed to be $\sim 0.3-0.5\,$dex even within the MW at present-time (Bolatto et al., 2008) implying $>0.6$ dex scatter. Likewise the density dependence in “Jeans Mass” models predicts $>0.3$ dex scatter even if all temperature variations, time variations, and progenitor-galaxy variations are neglected (e.g. if we use only the scatter in cloud densities observed in the solar neighborhood of the MW at the present instant).

Recall, the cloud properties we use to predict the IMF are measured at a density scale of $\sim 1000\,{\rm cm^{-3}}$ and mass scale $\sim 7000-56000\,M_{\odot}$. Obviously the clouds must continue to evolve and fragment to form actual stars – this is what our cloud-scale IMF models attempt to model. One might wonder, however, whether during this process some of the scatter might be reduced (if, for example, the clouds all converged to the same temperature eventually, owing to some additional physics). In the opacity-limited models, the equation of state (EOS) is specified (generally the cloud cools with $T\propto\rho^{-0.3}$ to some density, becomes approximately isothermal, then becomes adiabatic above the opacity-limit density), so this is already built into the model explicitly. In the “Jeans Mass” or “Turbulent/Sonic Mass” models, we have implicitly assumed an isothermal EOS within each cloud so their temperature was assumed to be constant throughout their evolution and set by the initial conditions. One might, therefore, consider a more complicated version of these models (different from the simple scalings used thus far). Let us assume star formation occurs above some critical density $\rho_{\rm crit}$ and the gas follows a polytropic EOS with index $\gamma$. The critical mass (Table 1) will then depend on $T_{\rm crit}(\rho=\rho_{\rm crit})=T_{\rm cloud}\,(\rho_{\rm crit}/\rho_{\rm cloud})^{\gamma-1}$, as $M_{\rm crit}\propto T_{\rm crit}^{\alpha}$ where $\alpha=3/2,\,1$ for the Jeans and Turbulent/Sonic models, respectively. Some simple algebra then gives us logarithmic variance in $M_{\rm crit}$, $S_{\log{M_{\rm crit}}}=\alpha\,(S_{\log{T_{\rm cloud}}}+\,[\gamma-1]\,S_{\log{\rho_{\rm cloud}}})$. Putting in the actual values (Fig. 8) this gives a dispersion $\sigma_{\log{M_{\rm crit}}}\approx 0.6,\,0.4$ dex for the Jeans and Turbulent/Sonic models (for any $\gamma\sim 0.5-1.5$). This reduces the predicted IMF variation, but still leaves it far larger than observed.

Thus we have shown that some additional physics on cloud or sub-cloud scales must be accounted for to reconcile the predictions of the “no-feedback” IMF models with the (weak) IMF variations observed in the MW. The “protostellar heating” models represent one physically-motivated class of models that do exactly this. Of course there may be others, but, broadly-speaking, they would need to either (a) strongly reduce the level of dependence of the predicted IMF on cloud properties (as the protostellar heating models do), or (b) strongly reduce the variation in GMC-scale properties predicted across cosmic time in the progenitor galaxies that form the MW. The latter is not impossible but seems to contradict the direct observations cited above, showing large variations in cloud properties in distant galaxies.

3.1 Caveats

Of course detailed, complex simulations (like the cosmological FIRE runs we are using) employ a large number of approximations to make problems numerically tractable. Although these simulations have been extensively vetted numerically (for details see Hopkins et al., 2017) some caveats worth noting include:

• •

Our analysis uses a somewhat arbitrary $n_{\rm crit}=1000\,\rm cm^{-3}$ minimum density threshold for star formation, based on numerical considerations. Using a much higher threshold would require much greater mass resolution (or else it would introduce severe numerical artifacts), which is not computationally feasible (these are the highest-resolution cosmological simulations of MW-mass galaxies ever run, at present). However, within the range we can probe, our results do not appear to depend sensitively on the density threshold or other numerical criteria for star formation.555In Orr et al. (2017), we show the results of a number of simulations where we re-run our m12i galaxy from $z\approx 0.1-0.0$, as in Fig. 4, but vary the numerical SF criteria. This includes changing the minimum SF density (from $\sim 10-1000\,{\rm cm^{-3}}$), removing requirements that the gas be molecular and/or self-gravitating, and changing the efficiency per free-fall time with which gas that meets this criteria will turn into stars (from $\sim 1-100\%$). We have re-run our analysis, restricted to just those stars formed over the period the simulations were re-run, and find these changes do not significantly influence the predicted IMF variations.

• •

In the simulations, gas elements are replaced by star particles instantaneously once all star formation criteria and timescales are satisfied, so star formation happens in discrete steps. In the large GMCs where most stars form ($\sim 10^{6}-10^{7}\,M_{\odot}$), this means that the first generation of stars formed can continue to alter the GMC properties while subsequent star formation occurs. However star formation in the smallest GMCs will be artificially “abrupt” (although GMCs with masses this low contribute negligibly to the variation in the IMF).

• •

Feedback processes from low-mass stars, e.g. proto-stellar outflows, are not explicitly included in the simulations. We only consider the effects of massive stars, which dominate on GMC scales provided there are sufficient stars to sample the IMF.

• •

The turbulent velocity dispersion in the code is calculated from using a kernel interpolation over the the relative velocities between the nearest $\sim 32$ resolution elements, after subtracting the coherent shear and contraction/expansion terms. This means that for very small GMCs with masses $\lesssim 10$ times the resolution, internal motions are not well-resolved. In this limit the general tendency is to under-estimate the turbulent velocity dispersions (see e.g. the detailed turbulence studies in Hopkins 2015). But again, these do not contribute significantly in our predictions.

• •

The simulations do not explicitly follow non-equilibrium chemistry (e.g. molecular hydrogen formation/destruction), instead relying on pre-tabulated equilibrium cooling rates as a function of density, temperature, metallicity, and the strength of the local radiation field in several bands. It was shown by Hopkins et al. (2012) that these approximations have little to no effect on galactic star formation properties but they could conceivably alter the scatter in small-scale cloud properties.

4 Conclusions

In this paper we explore the application of broad classes of IMF models to high-resolution fully-cosmological galaxy formation simulations. Stars at a some present-day location might have formed at very different times and places, in an environment radically different from today’s MW: only by using a cosmological simulation instead of local simulations or observations can we predict the properties of their progenitor star-forming clouds at these times and places, and therefore use these models to predict e.g. the variations in the IMFs predicted for old stellar populations in the present-day galaxy. This also provides an important consistency and validity check for future attempts to incorporate these IMF models into such simulations dynamically, as stellar feedback plays a critical role in the simulations and it, obviously, depends on the IMF.

In summary, we find that only models accounting for protostellar heating produce sufficiently weak IMF variations, in these simulations, to be compatible with observations. This discrepancy is not obvious in many previous studies (either analytic or idealized single-cloud simulations) as they artificially assume all clouds (at all locations and cosmic times) have the same temperature and obey the same linewidth-size relation (without scatter or systematic variation), whereas observations find significant variations in molecular gas temperatures and velocity dispersions (both within the MW and in nearby and high-redshift star-forming galaxies, which may more closely resemble the MW progenitors where these stars formed).

The protostellar heating models, on the other hand, actually predict IMF variations significantly below the observational upper limits (see Fig. 9). Additional sources of variance are therefore easily accommodated in these models, such as those that should come from a combination of (a) stochastic statistical sampling effects (see Bastian et al. 2010; these may be especially important in small clouds such as Taurus which are not resolved by our simulations, see Kraus et al. 2017), (b) measurement uncertainties, or (c) additional physics not accounted for by the model (e.g. bursty accretion or other physics may modify the radiative efficiency and heating effects of protostars, introducing some IMF variation).

In future work, we will examine whether the protostellar heating models considered in this study should produce observably-large IMF variation under more extreme conditions. Preliminary comparison of single-cloud conditions in Guszejnov et al. (2016) suggests these models can produce as much as factor $\sim 2$ shifts in the turnover mass under extreme starburst conditions analogous to Arp220, but this needs to be explored in more detail. We will also explore in more detail IMF shape variations, the predicted IMF in different sub-regions of the galaxy (e.g. the galactic nucleus), and the IMF in specific populations (e.g. metal-poor globular clusters versus present-day stellar populations).

Acknowledgments

Support for PFH, DG, and XM was provided by an Alfred P. Sloan Research Fellowship, NASA ATP Grant NNX14AH35G, and NSF Collaborative Research Grant #1411920 and CAREER grant #1455342. Numerical calculations were run on the Caltech compute cluster “Zwicky” (NSF MRI award #PHY-0960291) and allocation TG-AST130039 granted by the Extreme Science and Engineering Discovery Environment (XSEDE) supported by the NSF.