# The metallicity dependence of the stellar initial mass function

Tabassum S. Tanvir,1 Mark R. Krumholz${}^{\href
https://orcid.org/0000-0003-3893-854X\,1,2}$
1Research School of Astronomy and Astrophysics, Australian National
University, Canberra, ACT 2611, Australia
2Australian Research Council Centre of Excellence in All Sky Astrophysics
(ASTRO3D), Canberra, ACT 2611, Australia<EMAIL_ADDRESS>mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Dust is important for star formation because it is the crucial component that
couples gas to stellar radiation fields, allowing radiation feedback to
influence gas fragmentation and thus the stellar initial mass function (IMF).
Variations in dust abundance therefore provide a potential avenue by which
variation in galaxy metallicity might affect the IMF. In this paper we present
a series of radiation-magnetohydrodynamic simulations in which we vary the
metallicity and thus the dust abundance from 1% of Solar to 3$\times$ Solar,
spanning the range from the lowest metallicity dwarfs to the most metal-rich
early-type galaxies found in the local Universe. We design the simulations to
keep all dimensionless parameters constant so that the interaction between
feedback and star-forming environments of varying surface density and
metallicity is the only factor capable of breaking the symmetry between the
simulations and modifying the IMF, allowing us to cleanly isolate and
understand the effects of each environmental parameter. We find that at a
fixed surface density more metal-rich clouds tend to form a slightly more
bottom-heavy IMF than metal-poor ones, primarily because in metal-poor gas
radiation feedback is able to propagate further, heating somewhat larger
volumes of gas. However, shifts in IMF with metallicity at a fixed surface
density are much smaller than shifts with surface density at fixed
metallicity; metallicity-induced IMF variations are too small to explain the
variations in mass-to-light ratio reported in galaxies of different mass and
metallicity. We, therefore, conclude that metallicity variations are much less
important than variations in surface density in driving changes in the IMF and
that the latter rather than the former are most likely responsible for the IMF
variations found in early-type galaxies.

###### keywords:

magnetic fields — radiative transfer — stars: formation — stars: luminosity
function, mass function — stars: protostars — turbulence

††pubyear: 2015††pagerange: The metallicity dependence of the stellar initial
mass function–The metallicity dependence of the stellar initial mass function

##  1 Introduction

The stellar initial mass function (IMF) is the mass distribution of stars at
the point of their birth. The IMF is one of the most important distributions
in astrophysics because it at least partly determines everything from chemical
evolution to the strength of feedback during galaxy formation. The IMF has
been found to be nearly universal in the Milky Way and its closest neighbours,
which are the only locations where measurements of the IMF by direct star
counting are possible (Offner et al., 2014, and references therein). A number
of hypotheses have been proposed to explain this lack of variation, mostly
focusing on the universality of turbulence (Padoan et al., 1997; Padoan &
Nordlund, 2002; Hennebelle & Chabrier, 2008, 2009; Hopkins, 2012, 2013; Nam et
al., 2021) and the stabilising effects of stellar radiation feedback (e.g.,
Bate, 2009, 2012; Krumholz, 2011; Krumholz et al., 2012; Guszejnov et al.,
2016; Cunningham et al., 2018) or the isothermal-adiabatic transition (e.g.,
Lee & Hennebelle, 2018; Hennebelle et al., 2019). One important outcome of
these studies is that, while turbulence alone can explain why the high mass
portion of the IMF always has the same slope, some additional physical process
is likely required to explain the existence of an IMF peak at a particular
mass (e.g., Krumholz, 2014; Krumholz et al., 2016; Guszejnov et al., 2016,
2020). This finding is particularly significant because tentative evidence has
started to emerge that in the extreme star forming environments found in
early-type galaxies (ETGs), the IMF peak shifts slightly towards lower masses
than the IMF peak found in the local Universe. There are multiple lines of
evidence for such a shift derived from different measurement techniques –
spectroscopic (e.g., van Dokkum & Conroy, 2010; Spiniello et al., 2012; La
Barbera et al., 2013; Conroy et al., 2017), dynamical (e.g., Cappellari et
al., 2012; Newman et al., 2017; Oldham & Auger, 2018), and gravitational
lensing (Treu et al., 2010; Spiniello et al., 2015) – all pointing towards at
least qualitatively consistent conclusions; see Smith (2020) for a
comprehensive review. There is also much more tentative evidence for the
possibility that the IMF might be more bottom-light in ultra-faint dwarf
galaxies (e.g., Geha et al. 2013; Gennaro et al. 2018; however, see El-Badry
et al. 2017 for a contrary perspective).

These observations raise the theoretical question of which physical processes
might be responsible for inducing the shift in the IMF peak in ETGs, and at
least potentially in dwarfs as well. In Tanvir et al. (2022, hereafter Paper
I), we performed a series of carefully controlled radiation
magnetohydrodynamic (RMHD) simulations to study the role of the environment
and its interaction with stellar feedback mechanisms in setting the IMF peak.
We control these experiments by keeping all the dimensionless parameters
constant (e.g., virial parameter, Mach number, Alfvén Mach number) while
altering only one other parameter so that it is easier to deduce the role of
that parameter along with the feedback mechanism. In Paper I we explored the
role of surface density and showed that with increasing surface density the
IMF peak shifts towards lower masses than the ones found in the Milky Way, but
by less than would be expected purely from a shift in the mean Jeans mass with
surface density because of the enhanced effectiveness of stellar radiation
feedback in denser environments. The resulting shift in the IMF peak plausibly
matches what is required to explain the mass-to-light ratios measured in ETGs.

However, in Paper I, we did not alter another potentially important parameter
that differs between ETGs and other galaxies: the metallicity, and therefore
dust properties. Dust matters since it is the component of the ISM that
couples gas to stellar radiation. Some previous authors have conducted
numerical simulations to understand how changing metallicity, and therefore
changing dust abundance, might alter the IMF. For example, Bate (2019) carries
out radiation hydrodynamical simulations to study the metallicity-dependence
of the properties of stellar populations formed at metallicities $Z=0.01-3$
$\rm Z_{\odot}$, and finds that the stellar mass distribution is largely
insensitive to metallicity. This result is in qualitative agreement with a
number of similar, earlier numerical studies (Myers et al., 2011; Bate, 2014;
Bate & Keto, 2015). More recently, however, Guszejnov et al. (2022) found that
the characteristic mass of a forming stellar population does vary with both
metallicity and the strength of the interstellar radiation field (ISRF). They
find that with a typical Galactic ISRF, lower metallicities produce a higher
characteristic mass. Chon et al. (2021) also studied the IMF in metal-poor
environments and found that the number of low-mass stars increases with
metallicity and that the mass function is top-heavy in low-metallicity
environments. Sharda & Krumholz (2022) use analytic models to survey a wide
range of parameter space, and find that the characteristic stellar mass is
comparatively high at metallicities low enough that metal line cooling
dominates, but begins to decrease with metallicity once the metallicity is
high enough ($Z\gtrsim 0.01\mathrm{Z}_{\odot}$) for dust-gas coupling to
dominate gas thermodynamics Bate (2023) conducts simulations taking into
account how the combined effects of increasing cosmic microwave background
intensity at high redshift and variation in metallicity impacts the IMF. He
finds that for the CMB intensity at redshift $z=5$, increasing metallicity
increases the characteristic mass of stars.

While all these authors have studied the effects of varying metallicity, none
of the previous simulations (as opposed to analytic models) have done so for
the very high surface-density environments that likely characterised star
formation in ETGs. Moreover, the approach in most of these papers was to
consider a wide range of environmental changes – for example not just
metallicity, but also CMB intensity, ISRF intensity, total star-forming cloud
mass, etc. While this has the advantage of being maximally realistic since in
fact all of these parameters do vary with redshift and galactic environment,
it has the disadvantage of making the results extremely difficult to
interpret. Since these simulations vary many parameters at once, it is not
clear which of the many possible knobs that are being turned is responsible
for any particular aspect of the simulation outcomes, nor is it easy to
identify by what mechanism any particular knob operates.

In this paper we take a different approach, following on from that in Paper I:
we present a series of RMHD simulations where we explore the role of
metallicity variation in setting the IMF by keeping the dimensionless
parameters that describe our star-forming clouds constant except for
parameters that we vary systematically and one at a time so that we can create
a clean experiment that isolates the effects of environment and its
interaction with stellar feedback. A critical aspect of our simulation
approach is that our experiments are engineered so that, except for the one
change that we make, our simulations are all simply re-scaled versions of one
another, such that if we were to turn off stellar feedback and complex
thermodynamics, and simply adopt an isothermal equation of state, the results
of all simulations would be identical. This allows us to solve the problem of
interpretability since we can then unambiguously connect effects and causes.

We describe the numerical method and initial conditions we use to achieve
these effects in Section 2. In Section 3 we discuss the results from our
simulations and their implications for the physics behind the IMF. We
summarise our conclusions in Section 4.

##  2 Numerical Methods and Initial Condition

###  2.1 Numerical Methods

The numerical methods we employ in this study are identical to those used in
Paper I, and we refer interested readers to that paper for more detailed
information. Here we provide a brief summary. We carry out our simulations
using the orion2 adaptive mesh refinement code (Li et al., 2021). The code
uses the approach of Li et al. (2012) to solve the equations of ideal
magnetohydrodynamics with self-gravity (Truelove et al., 1998; Klein et al.,
1999) and radiation transfer (Krumholz et al., 2007) in the two-temperature,
mixed-frame, grey, flux-limited diffusion approximation.

The code includes sink particles (Krumholz et al., 2004) to replace regions
where protostars are forming and that are collapsing beyond our ability to
resolve. Each sink particle runs a one-zone protostellar evolution model as
described in Offner et al. (2009), which provides the instantaneous properties
of the star that it represents, such as radius, luminosity, and polytropic
index; these depend on the accretion history determined self-consistently from
the simulations. The luminosity of each sink particle is then used as a source
term in the radiative transfer equations. We also take into account the
feedback caused by protostellar outflows through momentum sources around each
sink particle. The outflow model we use is described in Cunningham et al.
(2011). In this model, whenever mass is accreted onto a sink particle, a
fraction $f_{w}$ of it is ejected back into the simulation in the form of an
outflow. This outflow material is launched with a speed of $v_{w}$. In our
simulation, we adopt the same wind model parameters used in Hansen et al.
(2012) and Cunningham et al. (2018), with $f_{w}=0.3$ and $v_{w}=\min(v_{\rm
kep},60\;\mathrm{km/s})$, where $v_{\rm kep}$ is the Keplerian speed at the
protostar’s surface (also determined self-consistently from the accretion
history via the protostellar evolution model). These parameters are based on
observations of the momentum budget of protostellar outflows (Richer et al.,
2000).

A critical component of the simulations is the opacity of the dusty gas, which
is responsible for coupling the gas flow to the radiation field generated by
the stars and by the thermal radiation from the dusty gas itself. As in Paper
I, we take the Rosseland and Planck mean opacity of the dusty gas as a
function of density and temperature from the tabulated results of Semenov et
al. (2003). In order to study the effects of varying the metallicity, we scale
these tabulated opacities by the metallicity relative to Solar, i.e., for runs
with metallicity $Z=0.1Z_{\odot}$ we take the opacities to be 10% of Semenov
et al.’s tabulated values. We should therefore understand the metallicity we
quote when describing our simulations as the dust metallicity; the gas-phase
metallicity is unimportant as long as the density is high enough for dust and
gas to be thermally coupled by collisions since in this case dust heating and
cooling completely dominates the gas thermodynamics.

There is a limitation to the radiation transfer method: we assume the dust and
gas temperatures are equal. This is generally a good assumption for gas and
dust temperature at densities above $\rm 10^{4}-10^{5}cm^{-3}$ (e.g.,
Goldsmith, 2001). However, in low-density regions, dust-gas collisions may
occur too infrequently to allow efficient gas-dust coupling, allowing gas to
be either hotter or cooler than the dust. This is however unlikely to be a
factor in determining the IMF since in regions where gas is collapsing and
fragmenting, the density is high enough for gas and dust to be well coupled.
Studies that have treated gas and dust temperatures separately (Bate & Keto,
2015) find that its effect on fragmentation is minimal (Bate, 2019).

###  2.2 Initial and boundary conditions

We construct our initial conditions following Paper I, and we refer readers to
that paper for full details. Our initial conditions are designed to produce a
carefully controlled experiment, whereby we hold all simulation parameters
fixed except for one, which we vary systematically in order to isolate the
effects of that parameter. In this study, we consider two different series of
runs. The “medium” density case (M runs hereafter) consists of a periodic
box111While all MHD quantities use periodic boundaries, the radiation field
uses Marshak boundaries, with a inward radiation flux corresponding to that of
an isotropic blackbody radiation field with a radiation temperature of 10 K;
see Paper I for full details. containing a mass $M_{\rm box}=1000$ M⊙ with
surface density $\Sigma=1$ g cm-2 (corresponding to a mean density
$\rho_{0}=7.0\times 10^{-19}$ g cm-3, box length $L=0.46$ pc), an initial
temperature $T_{0}=10$ K, an initial gas velocity dispersion $\sigma=2.4$ km
s-1, and an initially-uniform magnetic field of strength $B_{0}=0.73$ mG; for
this combination of parameters, the box free-fall time $t_{\mathrm{ff}}=160$
kyr, the crossing time $t_{\mathrm{cross}}=0.4$ Myr, the virial ratio
$\alpha_{\mathrm{vir}}=1$, the Mach number $\mathcal{M}=12.6$, and the plasma
$\beta=0.012$. The “low” density series (L runs hereafter) are a rescaled
version of the M runs with a surface density that is lower by a factor of
$f=1/2$, and a box length, volume density, and magnetic field strength that
are multiplied by factors of $f^{-1}$, $f^{2}$, and $f$ relative to the M
series, respectively; the gas velocity dispersion and temperature are
unchanged. These transformations have the property that they leave
$\alpha_{\mathrm{vir}}$, $\mathcal{M}$, and $\beta$ unchanged between the L
and M series, so the only dimensionless number that varies between the L and M
series is the optical depth of the gas.

In addition to varying surface density and thus the optical depth at fixed
metallicity, as in Paper I, in this study we also vary the metallicity
independently of the surface density. The metallicity range we explore is from
1% of Solar to 3$\times$ Solar; this covers the range from the lowest
metallicity dwarf galaxies in the local Universe (Madden et al., 2013) to the
most metal-rich early-types (Gu et al., 2022). We carry out runs at
$Z/Z_{\odot}=1\%$, $10\%$, $1$, and $3$ for both the L and M cases; we refer
to these runs as L1%, L10%, L1$\times$, and L3$\times$, and similarly for the
M series, and we summarise their full properties in Table 1. Note that the
L$1\times$ and M$1\times$ runs are identical to the L1 and M1 runs in Paper I.

Metallicity variation will change how the gas interacts with stellar radiation
feedback since metals (in the form of dust grains) are what couple the
radiation to the gas. They will also change how the gas interacts with
protostellar outflow feedback, since outflows shock the gas, and the rate at
which the gas is then able to cool back down to its equilibrium temperature
depends on the box optical depth. Since the dimensionless parameters are the
only ones that describe the system in the absence of stellar feedback
(radiation and protostellar outflow), any differences we find between the
simulations can only be a result of the interaction of stellar feedback with
the differing surface density and metallicity, and thus the differing optical
depth. By comparing to the Solar metallicity runs reported in Paper I, we can
further separate the effects of metallicity and surface density.

Name |  $M_{\rm box}$ (M⊙) |  $L$ (pc) |  $\rho_{0}$ (g cm-3) |  $B_{0}$ (mG) |  $t_{\rm ff}$ (kyr) |  $\Delta x$ (AU) |  $\rm t_{\rm cross}$ (Myr) |  $\Sigma$ (g cm-2) | $Z/Z_{\odot}$ | $\tau_{100K}$
---|---|---|---|---|---|---|---|---|---|---
L1% | 2000 | 0.92 | $\rm 1.74\times 10^{-19}$ | 0.36 | 160 | 46 | 0.4 | 0.5 | 0.01 | $5.38\times 10^{-3}$
L10% | 2000 | 0.92 | $\rm 1.74\times 10^{-19}$ | 0.36 | 160 | 46 | 0.4 | 0.5 | 0.1 | 0.054
L$1\times$ | 2000 | 0.92 | $\rm 1.74\times 10^{-19}$ | 0.36 | 160 | 46 | 0.4 | 0.5 | 1 | 0.54
L$3\times$ | 2000 | 0.92 | $\rm 1.74\times 10^{-19}$ | 0.36 | 160 | 46 | 0.4 | 0.5 | 3 | 1.6
M1% | 1000 | 0.46 | $\rm 6.96\times 10^{-19}$ | 0.73 | 80 | 23 | 0.2 | 1 | 0.01 | 0.0108
M10% | 1000 | 0.46 | $\rm 6.96\times 10^{-19}$ | 0.73 | 80 | 23 | 0.2 | 1 | 0.1 | 0.108
M$1\times$ | 1000 | 0.46 | $\rm 6.96\times 10^{-19}$ | 0.73 | 80 | 23 | 0.2 | 1 | 1 | 1.08
M$3\times$ | 1000 | 0.46 | $\rm 6.96\times 10^{-19}$ | 0.73 | 80 | 23 | 0.2 | 1 | 3 | 3.2
Table 1: Simulation parameters, from left to right: run name, mass in the
computational box, size of the computational box, mean density in the
computational box, initial magnetic field strength, mean-density free-fall
time, cell size at the finest AMR level, turbulent crossing time, surface
density, metallicity, and optical depth computed using the Planck-mean opacity
evaluated a temperature $T=100$ K. All simulations use the same initial
velocity dispersion $\sigma_{0}=2.4$ km s-1 and temperature $T_{0}=10$ K.

###  2.3 Resolution, refinement, and sink particles

In order to ensure that we can compare our current runs to those carried out
in Paper I, we use identical resolution, refinement, and sink particle
creation criteria, and we refer to that paper for a detailed description,
which we only summarise briefly here. The AMR hierarchy in these simulations
is set on a $512^{3}$ base grid which we denote $\mathcal{L}=0$. The
simulation takes place in two stages; during the first, “driving” phase, with
lasts for two crossing times, we disable self-gravity and radiation transport,
and drive the turbulence to a steady velocity dispersion, providing time for
the turbulence to achieve a statistically steady state. During the driving
phase we disable AMR, so no higher level exists. After the driving phase we
start the “collapse” phase, during which we disable driving and re-enable
self-gravity and radiation. Once we turn on gravity the grid is allowed to
adaptively refine to a maximum level $\rm\mathcal{L}_{max}=2$, and we refine
in any cell where the Jeans number

$J=\sqrt{\frac{G\rho\Delta x^{2}}{\pi c_{s}^{2}}}$ (1)

rises above $J=1/8$; here $\rho$ is the gas density, $\Delta x$ is a cell
size, and $c_{s}$ is the gas isothermal sound speed. We report the size of the
cells on the finest level in Table 1.

Sink particle formation is triggered in any zone on the finest AMR level where
the gas is dense enough to reach a local Jeans number $J>1/4$. Once a sink
particle is formed it interacts with the gas via gravity, accretion, and
stellar feedback only.

##  3 Results

We present the results of our simulations here. First, we give an overview of
the simulations in Section 3.1. Next, we discuss the mass distribution of the
sink particles and identify how different metallicity at a fixed surface
density impacts the IMF in Section 3.2. In Section 3.3 and Section 3.4 we
interpret these results in terms of the effects of radiation and protostellar
outflow feedback, and the interaction of these two feedback mechanisms with
gas of varying metallicity.

###  3.1 Overview of simulations

In Figure 1 and Figure 2 we show the column density and the density-weighted
temperature of runs L and M for our four different metallicities: 1$\%$,
10$\%$, $1\times$, and $3\times$ Solar metallicity. As in Paper I we show
these plots at matching SFE (star formation efficiency) instead of at matching
times since star formation occurs at different times in these runs. It is
clear from Figure 1 that turbulence has created dense filamentary structures
and that star formation activity is confined within these structures.
Morphologically the runs are very similar to one another, which is not
surprising since by construction in the absence of a feedback mechanism these
runs would be identical. There is one small difference visible between the two
sets of runs with different column densities: run L produces a filamentary
structure that is more straight and narrow than the filamentary structures
produced in run M. By contrast, at different metallicities but fixed surface
density the simulations are almost identical to one another. This indicates
that, at least with regard to morphology, metallicity has a relatively minor
impact.

Figure 1: Column densities of simulations L (top) and M (bottom) series at 5%
SFE at metallicities 1%, 10%, 1$\times$, 3$\times$ solar metallicity. The
colour scale goes from $\log(\Sigma/\Sigma_{0})=-1$ to $1$, where
$\Sigma_{0}=\rho_{0}L$ and $\rho_{0}$ is the mean density in the simulation
domain. Circles show star particles, and are colour-coded by mass $m_{\star}$
from $\log(m_{\star}/M_{\rm box})=-5$ to $-3$, where $M_{\rm box}$ is the
total mass of the simulation box.

The temperature structure of the gas shown in Figure 2 behaves quite
differently. Here differences with metallicity are much more apparent than in
Figure 1. For both the L and M series, lower metallicity runs are warmer
compared to the higher metallicity runs. This indicates that the lower
metallicity runs stellar radiation is able to escape further from the dense
regions around individual stellar sources, leading to more widely distributed
heating and a warmer mean temperature over most of the volume. This is in
contrast to the effects of surface density, in which lower surface density
runs are cooler because radiation is trapped less efficiently in the full
simulation box.

Figure 2: Same as Figure 1, but showing density-weighted projected temperature
rather than column density.

We show the time evolution of the star formation efficiency (SFE) in Figure 3,
and the total number of stars present as a function of the SFE in Figure 4; we
define the SFE as the ratio of total stellar mass in the simulation to the
total initial mass present in the volume. We find that, in all runs, once the
star formation activity begins it takes approximately 0.5 free fall times to
reach 5$\%$ SFE. There is very little difference between the runs with varying
metallicities at a fixed surface density. This indicates that whatever effects
the feedback mechanisms have on the star formation rate in the simulations are
independent of metallicity and surface density. By contrast, we see larger
differences in Figure 4, which shows the number of stars as a function of SFE.
At a fixed surface density, the two more metal-poor runs in both the L and M
series, corresponding to 1% and 10% of Solar metallicity, produce the same
number of stars at a given SFE. On the other hand, the runs with Solar and
$3\times$ Solar metallicity produce more stars at fixed SFE, indicating that
the stars formed in these runs have systematically lower masses; the effect is
not large, but is clear, particularly in the L runs. We examine these trends
further in the next section.

Figure 3: Star formation efficiency as a function of time since the formation
of the first star $t_{\star}$, measured in units of the free-fall time $t_{\rm
ff}$. The left panel shows the L series of runs at different metallicities,
and the right panel shows the corresponding results for the M series. Figure
4: Number of stars formed as a function of star formation efficiency. As in
Figure 3, the left panel shows the L run series and the right panel shows the
M series.

###  3.2 Stellar mass distribution

Figure 5: The top panels shows the evolution of the median of the sink
particle mass distribution for run L (left) and M (right) at different
metallicities in absolute mass expressed in M⊙. The lower panel shows the
evolution relative to the box mass $M_{\rm box}$. Here the medians are
calculated with respect to mass rather than number.

To explore the metallicity dependence of the stellar mass function further, we
present two sets of plots. The first, Figure 5, shows the evolution of the
median of the sink particle mass distributions for all the runs. As in Paper I
we measure the median with respect to the stellar mass rather than the number,
i.e., the median stellar mass $m_{*}$ is defined by the condition that half
the total stellar mass is found in stars with masses $<m_{*}$. Looking at
Figure 5 we see that the median masses have reached nearly steady values at
5$\%$ SFE. When measured on both absolute (top panel) and relative scales
(bottom panel) the L runs at varying metallicity show greater variation in
median mass than the M runs; however, the direction of variation is the same
in both sets of runs, which is that as the metallicity increases the median
mass decreases. However, we further note that the differences between the
different metallicities are smaller than the differences between runs L and M,
particularly when expressed in terms of absolute rather than relative mass.

Figure 6: Cumulative distribution function (CDF) of the sink particle masses
of the simulations, the top panel is the run L at different metallicities, and
the bottom panel is run M at different metallicities. The left side is the CDF
with respect to absolute stellar mass and the right side is the CDF for mass
measured relative to the box mass $M_{\rm box}$. In both the L and M series,
lower metallicity runs produce a slightly top-heavy IMF than the higher
metallicity runs

To investigate these differences further, in Figure 6 we show the cumulative
mass functions of the runs at 5$\%$ SFE on both absolute and relative mass
scales; as in Figure 5, we measure the CDF with respect to mass rather than
number, since this is much more numerically stable. Consistent with the trend
as a function of time shown in Figure 5, in the L runs as the metallicity
increases there is a slight shift in the CDF shape and mass range it covers,
with metal-poor runs showing slightly heavier mass distributions than metal-
rich runs at the high mass end. The trend is less visible, and may in fact
reverse, at the lowest stellar masses, so that the metal-rich distribution is
narrowest. The variations in run M are similar to or perhaps slightly smaller
than those in run L, and are in qualitatively the same direction. Consistent
with Paper I, we also see that, on an absolute scale, the L runs have a
slightly heavier IMF than the M runs.

###  3.3 The effect of radiation feedback at varying metallicity

Recall that our simulation series is constructed so that if the gas were
isothermal and not subject to protostellar feedback, the different metallicity
and surface density runs would all simply be re-scaled versions of the same
simulation, and we would therefore expect identical results. Thus the
differences in IMF we have observed between the runs must be a result of
feedback, and variations in IMF at fixed surface density must be due to the
interaction of feedback with the gas of varying metallicity. In this section,
we examine how radiation feedback alters gas temperature distribution in runs
of varying metallicity and how this, in turn, influences fragmentation, and in
the next section, we perform a similar exercise for protostellar outflow
feedback.

Figure 7 shows the gas mass distribution with respect to both density and
temperature in runs L and M at different metallicity, all at 5% SFE. From the
plot, it is clear that there are differences between the metal-poor and the
metal-rich runs. The most significant differences are in the high-density
parts of the distribution since this is the gas that is closest to star
formation and thus whose fragmentation will most directly affect the IMF. In
the metal-poor runs, we see an upturn in the minimum temperature of the dense
gas, while in the more metal-rich runs the temperature remains flat as high
density. This difference occurs because at low metallicity the gas is less
opaque and therefore is not able to shield itself from heating by the
radiation produced by nearby stars, and therefore becomes warmer as it
collapses. The most extreme case here is the M run at 1% of Solar metallicity,
where there is essentially no gas at densities $>10^{3}\rho_{0}$ and
temperatures $\lesssim 15$ K. On the other hand, the temperature distribution
at high density extends to noticeably higher temperatures in the metal-rich
runs than in the metal-poor ones, most likely for the same reason: because
when the gas is more opaque, radiation is more easily “bottled up”, leading to
hotter warm regions around protostars that have already started radiating.
While all of these temperature differences might seem minor, recall that the
Jeans mass varies as $T^{3/2}$, so a factor of $1.5$ corresponds to a factor
of $1.8$ in mass – comparable to or larger than the IMF shifts we have
observed between the runs, and of roughly the size required to explain the
observations of ETGs.

In order to understand how the opacity affects fragmentation more directly, we
show in Figure 7 the 1D mass-weighted cumulative distribution function of mass
with respect to $M_{\rm J}/M_{\rm box}$, where $M_{\mathrm{J}}$ is the Jeans
mass computed from the local gas density and temperature. For a detailed
description of how to construct this figure we refer readers to Paper I), but
to summarise, note that the horizontal axis in Figure 8 corresponds to the
dashed lines of constant $M_{\mathrm{J}}/M_{\mathrm{box}}$ in Figure 7. The
vertical axis in Figure 8 then shows what fraction of the mass in the
simulation lies to the right of the corresponding dashed line in Figure 7,
i.e., it shows that fraction of the mass in a given simulation has a Jeans
mass smaller than the indicated value, and thus has the potential to fragment
to produce stars with mass smaller than that value. From this plot, we see a
clear trend between the metal-poor and metal-rich runs. The metal-rich runs
tend to have a more mass at lower $M_{\rm J}/M_{\rm box}$ than the metal-poor
runs, consistent with our observation that the metal-rich runs to produce a
more bottom-heavy IMF. The location at which the dashed 1-1 line in the figure
intersects the CDF gives a rough estimate of the minimum object mass that is
possible for gravity to produce: to the left of this intersection point, the
box contains less than a single Jeans mass of material for which $M_{J}$ is
that small, and is therefore not capable of producing collapsed objects of
such small mass. The intersection point is on average shifted to slightly
lower mass at higher metallicity, indicating that radiation feedback is more
efficient in metal-poor runs at suppressing the formation of lower-mass
objects.

Figure 7: The joint distribution of normalised density $\rho/\rho_{0}$ and
temperature $T$ for runs L and M at different metallicity. The top row shows
the state of the L runs and the bottom row shows the M runs; metallicity
varies from 1% to 3$\times$ Solar from left to right, as indicated above the
columns. All plots show the state of the simulations when they reach 5% SFE.
The colour bar shows the mass in each density-temperature bin. Dashed lines
indicate loci of constant box-normalised Jeans mass $M_{\rm J}/M_{\rm box}$;
lines are spaced logarithmically at intervals of 0.5 dex, with the left-most
line corresponding to $\log M_{\rm J}/M_{\rm box}=-2$. Figure 8: Cumulative
distribution functions of mass with respect to $M_{\rm J}/M_{\rm box}$ for the
L run series (left) and the M run series (right). The dashed straight lines in
each panel show the one-to-one relation, and indicate a minimum condition for
fragmentation: for any mass $M_{\rm J}/M_{\rm box}$ for which the CDF falls
below the dashed line, there is less than a single Jeans mass of material with
$M_{\mathrm{J}}$ that small in the box, and thus it is impossible to create an
object of that mass via gravitational collapse.

###  3.4 The effect of outflow feedback

Figure 9: Accretion rate-weighted mean surface escape speed of stars formed
the simulations as a function of SFE. The top panel shows the L series of run
at varying metallicity, and the bottom panel shows the M series.

Figure 10: Outflow kinetic energy, momentum, and volume relative to the box
kinetic energy, momentum, and volume, all as a function of SFE. The top set of
panels shows the L series of runs and the bottom shows the M series.

As discussed in Paper I, radiation and outflow feedback play complementary
roles in shaping the IMF: the former sets a lower bound on the masses of stars
that can form, shaping the low-mass end of the IMF, while the latter inhibits
the accumulation of mass by the most massive stars, shaping the high-mass end
of the IMF. This difference explains why it is possible for the low-mass and
high-mass parts of the CDF to shift in opposite directions, as we observe in
Figure 6. In this case, if radiation feedback were the only actor in our
simulations, then we would expect a uniform shift of the IMF towards the
higher masses in the metal-poor runs. The actual pattern of IMF shift with
metallicity is more complex, indicating that protostellar outflow feedback is
significant as well, and that its dependence on metallicity is not identical
to that of radiation feedback.

To understand how outflow feedback affects our runs we first look at Figure 9,
which shows the evolution of accretion rate-weighted mean surface escape speed
in the runs as a function of SFE. This quantity matters because the outflow
prescriptions we use in these simulations link the outflow velocity to the
stellar surface escape speed, as is observed to be the case (Richer et al.,
2000, and references therein) and as is expected theoretically (Konigl &
Pudritz, 2000, and references therein), since outflows are launched from close
to the stellar surface. See Cunningham et al. (2011) for a detailed
description of our outflow prescription. From the plot, it is clear that
$v_{\mathrm{esc}}$ varies with surface density (run L has a higher surface
escape speed than run M) but for a fixed surface density there is no
systematic difference in the surface escape speed between the different
metallicity runs; this is not surprising, given that the primary determinant
of escape speed is the ratio of the time stars have had to contract toward the
main sequence to the stellar Kelvin-Helmholtz timescale; the former depends on
surface density but not on metallicity (since the simulation star formation
rates are almost independent of metallicity), and the latter to good
approximation depends on neither. Thus we find that at any given metallicity
the L series of runs has higher mean escape speed than the M series, since the
L runs take longer to reach a given SFE, but that there is no difference with
metallicity.

Thus the differences we find in how outflows influence the IMF at different
metallicity cannot be a result of differences in the outflow strength.
Instead, in order to understand how protostellar outflows are breaking the
symmetry in the simulations, we next examine Figure 10, which shows the total
kinetic energy, scalar momentum, and volume occupied by outflow material
normalised to the box kinetic energy ($E_{\mathrm{box}}=M_{\rm
box}\sigma_{v}^{2}/2$, where $\sigma_{v}$ is the initial 3D velocity
dispersion), scalar momentum by $p_{\mathrm{box}}=M_{\rm box}\sigma_{v}$ and
volume $V_{\mathrm{box}}=L^{3}$; we identify outflow material by a passive
scalar that we add to the gas launched into outflows, and refer readers to
Paper I for details on the procedure. Looking at this plot we see that in the
metal-rich runs the outflowing material on average has more kinetic energy and
scalar momentum, and occupies more volume, than for the more metal-poor runs
at the same surface density. In Paper I we found that greater outflow momentum
results in less efficient fragmentation of the gas, because the outflows punch
through their environments and escape more efficiently, therefore creating an
environment more favourable to the formation of higher-mass stars. However,
this does not appear to be the case here: the metal-rich runs have more
outflow momentum and kinetic energy, but have slightly fewer massive stars, as
can be seen by examining the CDFs shown in Figure 6. One possible explanation
is that the outflow momentum and kinetic energy vary between the L and M run
series for different reasons than they vary with metallicity at fixed surface
density. The outflow energy and momentum are larger in the L runs than the M
ones because the L runs have more powerful outflows due to the higher surface
escape speeds achieved by stars that have longer to contract toward the main
sequence. By contrast, the differences between the runs at different
metallicities are almost certainly driven by the different cooling rates of
the outflow gas and the dense cloud material with which it mixes. In the
metal-poor runs, cooling times are longer due to the lower metal content of
the gas, and as a result, regions of shocked cloud gas may build up that
pressure-confine the outflow more effectively, in turn making outflows less
efficient at breaking up the gas clumps. This allows slightly more massive
stars to form at the high-mass end of the IMF, as we observe in Figure 6.

###  3.5 Implications for the mass to light ratio in early type galaxies

One of the primary motivations of this study is to investigate the role of
metallicity in the variations of the IMF that have been observed in ETGs. The
mass-to-light ratio is the most direct line of evidence we can collect from
our simulations to study these variations. In this section we explore the
mass-to-light ratios of our simulations with varying metallicity. We use the
slug stellar population synthesis code (da Silva et al., 2012; Krumholz et
al., 2015) to generate isochrones at stellar population ages of 5 to 10 Gyr
and at the metallicities used in the simulations (see Paper I for a detailed
description of the procedure). We use the isochrones to calculate the mass-to-
light ratio of the stellar populations formed in our simulations at ages of 5
to 10 Gyr in the SDSS $r$ band, which is commonly used for measurements of
$M/L$ in ETGs.

The top panel of Figure 11 shows the mass-to-light evolution of the runs from
age 5 to 10 Gyr, while the lower panel shows the IMF mismatch parameter
$\alpha$, defined as the ratio of the actual $M/L$ ratio in the simulations to
the $M/L$ ratio expected for a population with a Solar metallicity and a
Chabrier (2005) IMF at the same age; this latter quantity is the index most
commonly used in observations to study IMF variations in ETGs, though we
emphasise that the absolute value here is less significant than the
differences between the runs due to the systematic uncertainties. For both the
absolute $M/L$ or the IMF mismatch parameter, we see that, at any given
metallicity, there is a difference of $\sim 0.3-0.5$ dex in mass to light
ratio between the L and M runs; this finding is consistent with that in Paper
I for the Solar metallicity case, and extends the results to non-Solar
metallicities. At a fixed surface density, by contrast, there is a smaller
difference in the mass-to-light ratio between the runs, consistent with our
finding that surface density is a more influential factor than metallicity
when it comes to determining the IMF.

Moreover, Figure 11 in fact somewhat exaggerates the effect of metallicity,
because some of the difference in mass-to-light ratio between the differing
metallicity runs occurs because the metallicity itself affects stellar
evolution and atmospheres; thus we would not expect to identical $M/L$ values
for two populations of different metallicity even if they had exactly the same
IMF. To remove this effect, in Figure 12 we show the same quantities as in
Figure 11, but where we have calculated the mass-to-light ratio for all runs
using the isochrones for Solar metallicity; while this is clearly artificial,
it enables us to isolate the effects of metallicity on the IMF from the
effects on stellar evolution and atmospheres. In Figure 12 the runs at
different metallicity but the same surface density cluster even more tightly,
and are further from the corresponding runs at different surface density, than
in Figure 11. This reinforces our conclusion that the effects of metallicity
on the IMF and the observational diagnostics used to assess it are fairly
minor compared to the effects of surface density.

Our qualitative conclusion that metallicity effects cannot be responsible for
the IMF variations seen in ETGs is the same as that reached in a recent paper
by Bate (2023), who explored the role of varying metallicity on the stellar
properties at redshift $z=5$. Bate found that at this redshift, with its
higher CMB temperature, increasing metallicity leads to an increase in the
characteristic mass of stars, exactly the opposite of what would be required
to reproduce observed IMF variations in ETGs. However, our conclusions differ
somewhat in detail; Bate varied both the metallicity and the CMB temperature,
whereas we, in keeping with our philosophy of controlled experiments, have
varied only the former, a decision that is likely to be significant at
$z\gtrsim 3$, when the CMB temperature begins to exceed our adopted background
IR radiation field temperature of 10 K. Thus there is no contradiction between
the findings of Bate (2023) and of this work, since the experiments we have
carried out are different. On the other hand, our finding that there is a very
weak IMF variation with metallicity at fixed surface density is qualitatively
consistent with the findings of Bate (2019), though they explore only a single
surface density case, one similar to our L series.

Figure 11: The top panel shows the mass to light ratio $(M/L)_{r}$ in the SDSS
$r$ band as a function of population age computed for the stellar population
in formed in each of the simulations, using isochrones computed with the same
metallicity as used in the simulations. The black lines represent the $L$
series of runs and the red lines represent the $M$ series runs, respectively,
with different line styles corresponding to different metallicity as indicated
in the legend. The bottom panel shows the IMF mismatch parameter
$\alpha=(M/L)/(M/L)_{Z_{\odot},\rm Chabrier}$, where
$(M/L)_{Z_{\odot},\mathrm{Chabrier}}$ is the mass to light ratio expected for
a Solar-metallicity stellar population with a Chabrier (2005) IMF, again as a
function of stellar population age.

Figure 12: Same as Figure 11, but now computed using the same set of
isochrones (corresponding to those for $Z=Z_{\odot}$ for all simulations),
regardless of the simulation metallicity.

.

##  4 Conclusions

In this paper, we present a set of radiation-magnetohydrodynamic simulations
of star formation including radiation and protostellar outflow feedback from
young stars. We carry out a systematic exploration of how the mass function of
the stars formed in the simulations varies as a function of surface density
and metallicity, exploring surface densities from those typical of Milky Way-
like galaxies to those typical of starbursts or early-type galaxies, and
metallicities that range from those typical of the most metal-poor local
dwarfs, $Z=\mathrm{Z}_{\odot}/100$, to those typical of the centres of the
most metal-rich early-types, $Z=3\mathrm{Z}_{\odot}$. The setup of the
simulations allows us to separate out the effects of metallicity at fixed
surface density and of surface density at fixed metallicity, and to identify
the specific mechanisms by which these parameters interact with stellar
feedback. This extends the results of Tanvir et al. (2022), who used the same
methodology to explore the effects of surface density alone, without
metallicity variation. As in that paper, they key to approach is to set up a
series of simulations where we keep all dimensionless parameters except the
cloud optical depth (which depends on both surface density and metallicity)
constant, so that, for isothermal gas in the absence of stellar feedback, all
the simulations would just be rescaled versions of one another, and would
produce identical results.

Overall we see a trend whereby metal-poor cases produce a slightly heavier
mass distribution over most of the IMF than metal-rich cases. However, this
trend disappears or reverses at the lower mass end of the IMF, so that metal-
rich cases also have a slightly narrower IMF overall. We attribute these
shifts to the contrasting ways that metallicity variation interacts with
different types of stellar feedback. This high mass end of the IMF is
sensitive primarily to the effects of protostellar outflow feedback inhibiting
the most massive objects from accreting, something that appears to happen
slightly more efficiently when the metallicity is higher and cooling of
outflow-shocked gas is more rapid, suppressing the formation of more massive
objects. By contrast the remainder of the IMF is shaped primarily by radiation
feedback suppressing fragmentation and preventing small objects from forming,
something that appears to happen slightly more efficiently at lower
metallicity. The net effect is that the IMF shifts to somewhat lower masses,
and is also somewhat narrower, when the metallicity is higher.

While these differences with metallicity provide interesting insight into how
feedback interacts with the star-forming environment, they are ultimately of
rather minor importance. This is because on an absolute scale the biggest
differences we see in the IMF are set by variations in surface density. Our
low surface density cases produce a slightly heavier IMF than our higher
surface ones, and the differences are substantially larger than the subtle
variations induced by metallicity. When we explore how these IMF variations
compare to those required to explain observed variations in the mass-to-light
ratio in early-type galaxies compared to spiral galaxies, we find that
metallicity-induced IMF changes are too small to explain the observations,
whereas surface density-induced ones are at the right level. We, therefore,
conclude that surface density is a more important factor than metallicity in
determining the stellar IMF and that the differences observed between the IMFs
of spiral and early-type galaxies are most likely an effect of interstellar
pressure and therefore surface density, not an effect of metallicity.

## Acknowledgements

TST and MRK acknowledge support from the Australian Research Council through
Laureate Fellowship award FL220100020. This research was undertaken with the
assistance of resources and services from the National Computational
Infrastructure (NCI), which is supported by the Australian Government, and by
resources provided by the Pawsey Supercomputing Research Centre with funding
from the Australian Government and the Government of Western Australia.

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Bate (2009) Bate M. R., 2009, MNRAS, 392, 1363
  * Bate (2012) Bate M. R., 2012, MNRAS, 419, 3115
  * Bate (2014) Bate M. R., 2014, MNRAS, 442, 285
  * Bate (2019) Bate M. R., 2019, MNRAS, 484, 2341
  * Bate (2023) Bate M. R., 2023, MNRAS, 519, 688
  * Bate & Keto (2015) Bate M. R., Keto E. R., 2015, MNRAS, 449, 2643
  * Cappellari et al. (2012) Cappellari M., et al., 2012, Nature, 484, 485
  * Chabrier (2005) Chabrier G., 2005, The Initial Mass Function: From Salpeter 1955 to 2005. p. 41, doi:10.1007/978-1-4020-3407-7_5
  * Chon et al. (2021) Chon S., Omukai K., Schneider R., 2021, MNRAS, 508, 4175
  * Conroy et al. (2017) Conroy C., van Dokkum P. G., Villaume A., 2017, ApJ, 837, 166
  * Cunningham et al. (2011) Cunningham A. J., Klein R. I., Krumholz M. R., McKee C. F., 2011, ApJ, 740, 107
  * Cunningham et al. (2018) Cunningham A. J., Krumholz M. R., McKee C. F., Klein R. I., 2018, MNRAS, 476, 771
  * El-Badry et al. (2017) El-Badry K., Weisz D. R., Quataert E., 2017, MNRAS, 468, 319
  * Geha et al. (2013) Geha M., et al., 2013, ApJ, 771, 29
  * Gennaro et al. (2018) Gennaro M., et al., 2018, ApJ, 855, 20
  * Goldsmith (2001) Goldsmith P. F., 2001, ApJ, 557, 736
  * Gu et al. (2022) Gu M., Greene J. E., Newman A. B., Kreisch C., Quenneville M. E., Ma C.-P., Blakeslee J. P., 2022, ApJ, 932, 103
  * Guszejnov et al. (2016) Guszejnov D., Krumholz M. R., Hopkins P. F., 2016, MNRAS, 458, 673
  * Guszejnov et al. (2020) Guszejnov D., Grudić M. Y., Hopkins P. F., Offner S. S. R., Faucher-Giguère C.-A., 2020, MNRAS, 496, 5072
  * Guszejnov et al. (2022) Guszejnov D., Grudić M. Y., Offner S. S. R., Faucher-Giguère C.-A., Hopkins P. F., Rosen A. L., 2022, MNRAS, 515, 4929
  * Hansen et al. (2012) Hansen C. E., Klein R. I., McKee C. F., Fisher R. T., 2012, ApJ, 747, 22
  * Hennebelle & Chabrier (2008) Hennebelle P., Chabrier G., 2008, ApJ, 684, 395
  * Hennebelle & Chabrier (2009) Hennebelle P., Chabrier G., 2009, ApJ, 702, 1428
  * Hennebelle et al. (2019) Hennebelle P., Lee Y.-N., Chabrier G., 2019, ApJ, 883, 140
  * Hopkins (2012) Hopkins P. F., 2012, MNRAS, 423, 2016
  * Hopkins (2013) Hopkins P. F., 2013, MNRAS, 430, 1653
  * Klein et al. (1999) Klein R. I., Fisher R. T., McKee C. F., Truelove J. K., 1999, Gravitational Collapse and Fragmentation in Molecular Clouds with Adaptive Mesh Refinement HydrodynamicsCD. p. 131, doi:10.1007/978-94-011-4780-4_44
  * Konigl & Pudritz (2000) Konigl A., Pudritz R. E., 2000, Protostars and Planets IV, pp 759–+
  * Krumholz (2011) Krumholz M. R., 2011, ApJ, 743, 110
  * Krumholz (2014) Krumholz M. R., 2014, Phys. Rep., 539, 49
  * Krumholz et al. (2004) Krumholz M. R., McKee C. F., Klein R. I., 2004, ApJ, 611, 399
  * Krumholz et al. (2007) Krumholz M. R., Klein R. I., McKee C. F., 2007, ApJ, 656, 959
  * Krumholz et al. (2012) Krumholz M. R., Klein R. I., McKee C. F., 2012, ApJ, 754, 71
  * Krumholz et al. (2015) Krumholz M. R., Fumagalli M., da Silva R. L., Rendahl T., Parra J., 2015, MNRAS, 452, 1447
  * Krumholz et al. (2016) Krumholz M. R., Myers A. T., Klein R. I., McKee C. F., 2016, MNRAS, 460, 3272
  * La Barbera et al. (2013) La Barbera F., Ferreras I., Vazdekis A., de la Rosa I. G., de Carvalho R. R., Trevisan M., Falcón-Barroso J., Ricciardelli E., 2013, MNRAS, 433, 3017
  * Lee & Hennebelle (2018) Lee Y.-N., Hennebelle P., 2018, A&A, 611, A89
  * Li et al. (2012) Li P. S., Martin D. F., Klein R. I., McKee C. F., 2012, ApJ, 745, 139
  * Li et al. (2021) Li P., et al., 2021, The Journal of Open Source Software, 6, 3771
  * Madden et al. (2013) Madden S. C., et al., 2013, PASP, 125, 600
  * Myers et al. (2011) Myers A. T., Krumholz M. R., Klein R. I., McKee C. F., 2011, ApJ, 735, 49
  * Nam et al. (2021) Nam D. G., Federrath C., Krumholz M. R., 2021, MNRAS, 503, 1138
  * Newman et al. (2017) Newman A. B., Smith R. J., Conroy C., Villaume A., van Dokkum P., 2017, ApJ, 845, 157
  * Offner et al. (2009) Offner S. S. R., Klein R. I., McKee C. F., Krumholz M. R., 2009, ApJ, 703, 131
  * Offner et al. (2014) Offner S. S. R., Clark P. C., Hennebelle P., Bastian N., Bate M. R., Hopkins P. F., Moraux E., Whitworth A. P., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI. p. 53 (arXiv:1312.5326), doi:10.2458/azu_uapress_9780816531240-ch003
  * Oldham & Auger (2018) Oldham L., Auger M., 2018, MNRAS, 474, 4169
  * Padoan & Nordlund (2002) Padoan P., Nordlund Å., 2002, ApJ, 576, 870
  * Padoan et al. (1997) Padoan P., Nordlund A., Jones B. J. T., 1997, MNRAS, 288, 145
  * Richer et al. (2000) Richer J. S., Shepherd D. S., Cabrit S., Bachiller R., Churchwell E., 2000, in Mannings V., Boss A. P., Russell S. S., eds, Protostars and Planets IV. p. 867 (arXiv:astro-ph/9904097)
  * Semenov et al. (2003) Semenov D., Henning T., Helling C., Ilgner M., Sedlmayr E., 2003, A&A, 410, 611
  * Sharda & Krumholz (2022) Sharda P., Krumholz M. R., 2022, MNRAS, 509, 1959
  * Smith (2020) Smith R. J., 2020, ARA&A, 58, 577
  * Spiniello et al. (2012) Spiniello C., Trager S. C., Koopmans L. V. E., Chen Y. P., 2012, ApJ, 753, L32
  * Spiniello et al. (2015) Spiniello C., Koopmans L. V. E., Trager S. C., Barnabè M., Treu T., Czoske O., Vegetti S., Bolton A., 2015, MNRAS, 452, 2434
  * Tanvir et al. (2022) Tanvir T. S., Krumholz M. R., Federrath C., 2022, MNRAS, 516, 5712
  * Treu et al. (2010) Treu T., Auger M. W., Koopmans L. V. E., Gavazzi R., Marshall P. J., Bolton A. S., 2010, ApJ, 709, 1195
  * Truelove et al. (1998) Truelove J. K., Klein R. I., McKee C. F., Holliman John H. I., Howell L. H., Greenough J. A., Woods D. T., 1998, ApJ, 495, 821
  * da Silva et al. (2012) da Silva R. L., Fumagalli M., Krumholz M., 2012, ApJ, 745, 145
  * van Dokkum & Conroy (2010) van Dokkum P. G., Conroy C., 2010, Nature, 468, 940