Implementing Dust Shielding as a Criteria for Star Formation

TL;DR

This paper introduces a dust shielding-based star formation model for galaxy simulations, showing it produces similar galaxy properties to H2-based models and aligns with observed star formation laws.

Contribution

The paper presents a novel sub-grid model linking star formation to dust shielding, tested against other criteria in galaxy simulations.

Findings

Shielding model allows star formation at higher temperatures and lower densities.

Both shielding and H2-based models produce similar galaxy properties.

Models reproduce the observed Kennicutt-Schmidt law.

Abstract

Star formation is observed to be strongly correlated to dense regions of molecular gas. Although the exact nature of the link between star formation and molecular hydrogen is still unclear, some have suggested that shielding of dense gas by dust grains is the key factor enabling the presence of both. We present a sub-grid model for use in galaxy formation simulations in which star formation is linked explicitly to local dust shielding. We developed and tested our shielding and star formation models using smoothed particle hydrodynamic simulations of solar and sub-solar metallicity isolated Milky Way-mass disk galaxies. We compared our dust shielding-based star formation model to two other star formation recipes that used gas temperature and H$_2$ fraction as star formation criteria. We further followed the evolution of a dwarf galaxy within a cosmological context using both the shielding and H$_2$-based star formation models. We find that the shielding-based model allows for star formation at higher temperatures and lower densities than a model in which star formation is tied directly to H$_2$ abundance, as requiring H$_2$ formation leads the gas to undergo additional gravitational collapse before star formation. However, the resulting galaxies are very similar for both the shielding and H$_2$-based star formation models, and both models reproduce the resolved Kennicutt-Schmidt law. Therefore, both star formation models appear viable in the context of galaxy formation simulations.

Tables

Table 1. Table 1: Star formation and H 2 formation parameters of the models analyzed in this paper. L s is the shielding length used to calculate H 2 abundance and (if applicable) star formation. λ J subscript 𝜆 𝐽 \lambda_{J} is the Jeans length of a gas particle, while h ℎ h is the smoothing length. T max is the maximum temperature at which gas particles can form stars. c ∗ superscript 𝑐 ∗ c^{\ast} is the star formation efficiency factor. For all models, the minimum density at which stars can form is ρ m ​ i ​ n = .1 subscript 𝜌 𝑚 𝑖 𝑛 .1 \rho_{min}=.1 amu cm -3 .

Model	SF recipe	Ls	Tmax [K]	$c^{\ast }$
TC	Temperature cap	${\lambda }_J$ (T $\le 40$)	1000	$c_0^{\ast }$
H2	H2 fraction	${\lambda }_J$ (T $\le 40$)	15000	$c_0^{\ast }\left(X_{H_2}\right)$
S	Shielding fraction	${\lambda }_J$ (T $\le 40$)	15000	$c_0^{\ast }\left(1-S_d\right)$

Equations

$$S_d=e^{-{\sigma }_{d_{\text{eff}}}Z/Z_{\odot }(N_{HI}+2N_{H2})}$$

$$S_{H_2}=\frac{1-{\omega }_{H_2}}{(1+x^2)}+\frac{{\omega }_{H_2}}{(1+x{)}^{(1/2)}}{e}^{-0.00085(1+x{)}^{1/2}}$$

$$L_{shield}=\sqrt{\frac{15k_{B}T}{4\pi G{m}_H\rho }}$$

$$p=\frac{m_{gas}}{m_{star}}(1-e^{-c^{\ast }\Delta t/t_{\text{form}}}),$$

$$c^{\ast }=c_0^{\ast }{X}_{H_2}$$

$$c^{\ast }=c_0^{\ast }(1-S_d)$$

$$t_{chem}=\frac{1}{{\bar{n}}_{H}RC}.$$

$$t_{ff}=\sqrt{\frac{3\pi }{32G{\bar{n}}_H{\mu }_H{m}_H}}$$

$$\frac{t_{chem}}{t_{ff}}=24Z^{-1}{C}^{-1}{n}_o^{-1/2}$$

Full text

hydrodynamics — galaxies: dwarf — galaxies: evolution — methods: numerical — stars: formation

1 Introduction

Processes determining star formation drive the evolution of the stellar content in galaxies. Mimicking these processes in galaxy evolution simulations is critical to correctly modeling the stellar mass and distribution. However, since star formation happens on scales too small to be resolved in simulations of galaxy evolution, stochastic sub-grid models must be used to model it instead. It is therefore imperative to determine which models best mimic the actual process of star formation and produce star formation rates and locations that agree with observations from real galaxies. Additionally, comparing the effects of different star formation sub-grid models can reveal the primary drivers of star formation on 100 parsec scales.

Considerable observational evidence links the formation of stars to the surface density of gas in the star-forming region. Schmidt (1959) suggested that this link takes the form of $\rho_{SFR}\propto\rho_{gas}\,^{N}$ and suggested that $N\approx 2$. Kennicutt (1998) refined this relationship through observations of spiral and starburst galaxies, yielding the power law $\Sigma_{SFR}\propto\Sigma_{gas}\,^{1.4}$. However, more recent observational evidence suggests that star formation is linked specifically to the density of molecular gas. Wong & Blitz (2002) found that star formation rate surface density is correlated more closely with the surface density of molecular hydrogen than to the total gas surface density, a result which has been replicated by others (Boissier et al., 2003; Heyer et al., 2004). Similarly, Kennicutt et al. (2007) observed a strong correlation between $\Sigma_{SFR}$ and $\Sigma_{H2}$, but no significant correlation between $\Sigma_{SFR}$ and $\Sigma_{HI}$. Bigiel et al. (2008) also showed that the relationship between $\Sigma_{SFR}$ and $\Sigma_{gas}$ varied dramatically, while $\Sigma_{SFR}$ and $\Sigma_{H2}$ were consistently related by a Schmidt-type power law with $N=1.0\pm 0.2$. Various other studies have also found a strong and approximately linear relationship between $\Sigma_{SFR}$ and $\Sigma_{H2}$ (Bigiel et al., 2011; Schruba et al., 2011; Leroy et al., 2013), confirming the link between star formation and H2, rather than total gas.

One possible explanation for this link is that molecular hydrogen is required for star formation to occur. For example, Krumholz et al. (2009) demonstrated that the observed relationship between galaxies’ star formation rate and their atomic and molecular gas content can be explained by a model in which stars formed only out of molecular gas, at a rate given by  1% of the mass per free-fall time. Moreover, simulations in which star formation is explicitly linked to local H2 abundance have produced galaxies consistent with the observed Kennicutt-Schmidt and Tully-Fisher relations (Gnedin et al., 2009; Kuhlen et al., 2012; Christensen et al., 2012), as well as the observed relations between ISM pressure and molecular fraction (Robertson & Kravtsov, 2008).

However, although H2 traces star formation, it need not be a direct prerequisite for star formation to occur, as cooling in dense clouds of gas occurs primarily by C${}_{\text{II}}$ or CO lines. Indeed, by modeling star formation within a dense cloud of gas, Glover & Clark (2012) found that the presence of molecules in a gas cloud does not determine the cloud’s ability to form stars. Instead, they suggested that as long as a cloud of gas is dense and gravitationally unstable, cooling by H2 is not necessary for star formation to occur. Therefore, it is probable that molecular hydrogen is not required for star formation but instead simply correlated with it. For example, Mac Low & Glover (2012) have suggested that star formation and the presence of molecular gas both occur due to gravitational instability rather than being causally linked. Assuming such a situation, Krumholz (2012) demonstrated that in low-metallicity environments, below a few percent of solar metallicity, stars should be able to form without the presence of molecular hydrogen. Since the free-fall time of the gas will be longer than the thermal equilibrium timescale but shorter than the chemical equilibrium timescale, gas collapse will occur before the formation of H2. This theory is supported by the high atomic, low molecular content of gamma-ray burst host galaxies, which may arise when low-metallicity gas accretes onto a galaxy, then cools and collapses into stars before H2 can be formed (Michałowski et al., 2015).

The possibility that molecular hydrogen is not necessary for star formation has implications both for our understanding of galaxy formation and for our ability to accurately model star formation in low-density and low-metallicity environments. In particular, the low metallicity of the high-redshift universe makes it an extreme environment for understanding star-formation. Even after enrichment from population III stars, requiring H2 for star formation could have profound effects on the high-$z$ luminosity function, especially at the low-mass end (Kuhlen et al., 2012). Low-H2 star-forming environments also exist at low redshifts. For example, although the outer disks of spiral galaxies lack significant amounts of molecular gas, star formation can still occur, albeit inefficiently (Ferguson et al., 1998; Thilker et al., 2005; Bigiel et al., 2010). The same is true for the low-surface density, low-metallicity environments of dwarf galaxies (e.g. van Zee et al., 1997; Leroy et al., 2006; Melena et al., 2009; Roychowdhury et al., 2009). In order to accurately model these environments, we must understand the process of star formation in environments with low amounts of H2.

If molecular hydrogen is not a prerequisite for star formation, other factors must be considered. One potential explanation is that shielding by dust leads to both the formation of molecular hydrogen and the formation of cold gas necessary for gravitational collapse and star formation. Schaye (2004) argued that the transition of gas to the cold phase at high column densities both allows star formation to occur and coincides with an increase in the molecular fraction, as H2 forms more readily in dense gas and contributes to further cooling. Using detailed simulations of chemical evolution within collapsing gas clouds, Glover & Clark (2012) found that the column density of gas clouds and the associated dust shielding are the determining factors for both the creation of molecules and star formation. Dust shielding both reduces photodissociation of H2 by Lyman-Werner (LW) radiation and blocks photoheating by the Interstellar Radiation Field (ISRF); this allows cold, dense gas to form in the same regions where molecular hydrogen is able to form and remain. They argue that shielding of the ISRF, not the presence of H2 or CO, is the most important factor enabling star formation in gas clouds, a conclusion which has also been reached by others (e.g. Krumholz et al., 2011; Clark & Glover, 2014).

The evidence for a strong connection between dust shielding and star formation suggests that the accuracy and realism of numerical simulations may be increased by making star formation explicitly dependent on dust shielding. Such a model might allow for star formation to occur in areas or on timescales that a H2-based recipe would preclude. Here we implement a dust shielding-based star formation recipe and compare ISM and star formation properties of the resulting simulated galaxies to those produced using a temperature ceiling model and a H2-based model. We further analyze the effects of shielding versus H2-based star formation on a cosmological simulation of a dwarf galaxy by comparing disk size, star formation history, and the location of star formation. The comparison of shielding to other criteria enables a greater understanding of the link between dust shielding and star formation, and therefore how to more accurately model galaxy formation.

2 Methodology

We implement a novel method for dust shielding-based star formation for use within smoothed particle hydrodynamic galaxy formation simulations, including incorporating an updated method for calculating shielding. Below, we describe the isolated disk galaxy and cosmological galaxy formation simulations used to test our star formation and shielding models (§2.1). We then introduce the smoothed particle hydrodynamic code used (§2.2) and describe our implementation of shielding (§2.3) and star formation (§2.4).

2.1 Initial conditions

To test our shielding and star formation models, we used simulations of a $1.4\times 10^{12}$ M⊙ isolated Milky Way-like disk galaxy within a Navarro et al. (1997) dark matter halo. Following the work of Kaufmann et al. (2007) and Stinson et al. (2007), disks were formed by allowing a cloud of virialized gas to collapse within a virialized dark matter halo with a concentration of $c=8$, where $c=R_{vir}/R_{s}$, and $R_{s}$ is the scale radius of the Navarro-Frenk-White (NFW) profile. The initial halo had a gas fraction of 10%, and the gas particles were given an initially uniform circular velocity such that the dimensionless spin parameter $\lambda=(j_{gas}|E|^{1/2})/(GM^{3/2})=0.039$, where $j_{gas}$ is the average specific angular velocity of the gas, and E and M are the total energy and mass of the halo. The initial masses of the star, gas, and dark matter particles were $3.4\times 10^{4}$ M⊙, $1.4\times 10^{5}$ M⊙, and $1.3\times 10^{6}$ M⊙ respectively. The gravitational softening length was 206 pc and the smoothing lengths had a minimum value of 0.1 times the softening length. After 1 Gyr, following the formation of a stable disk, we scaled the gas metallicities to create one simulation with average ISM metallicity of 1.0 Z⊙ and one with 0.1 Z⊙. These two isolated disk galaxies formed the initial conditions for testing our star formation and shielding models. We then enabled H2 physics along with our various shielding and star formation models described below, and integrated the simulations for an additional 1.0 Gyr.

We also tested the effects of our shielding implementation using a zoom-in (Katz & White, 1993) cosmological simulation of a dwarf galaxy, $M_{vir}=4\times 10^{10}$ M⊙. A dwarf galaxy was chosen to highlight the impact of our model in a low-metallicity environment. The initial conditions consisted of a 25-Mpc box selected to surround a halo chosen from a low-resolution, dark-matter-only simulation. The cmbfast code was used to calculate the initial power spectrum used for the linear density field. We assumed a $\Lambda$CDM cosmology with values from the Wilkinson Microwave Anisotropy probe (Spergel et al., 2003). The masses of the star, gas, and dark matter particles were 1000, 3000, and 16,000 solar masses respectively, and the softening lengths were 87 pc.

The initial conditions for this simulation were first used in Governato et al. (2010) and produced a bulgeless dwarf galaxy. Christensen et al. (2012) simulated the same galaxy with the addition of metal line cooling and H2-based star formation. Compared to a simulation with metal-line cooling and density-based star formation, they found that the addition of H2-based star formation produced a clumpier ISM and on-going star formation at $z=0$.

2.2 General description of code

The simulations were performed using the galaxy formation code ChaNGa, a N-body smoothed particle hydrodynamics (SPH) code (N-Body Shop, 2011; Menon et al., 2014). ChaNGa is a highly scalable successor to gasoline (Wadsley et al., 2004), itself an extension of the parallel, gravity tree code PKDGRAV (Stadel & Gerhard, 2002). ChaNGa employs a modern implementation of the SPH kernel that uses a geometric mean density in the force expression (Ritchie & Thomas, 2000; Menon et al., 2014; Wadsley et al., 2017), thus eliminating artificial gas surface densities.

ChaNGa incorporates a suite of physical modules for thermodynamics, star formation, and stellar feedback, as described below. Non-equilibrium ion abundances, including H2, are calculated by integrating over the H and He chemical networks. Heating and photoionization are determined using a cosmic ultraviolet background (Haardt & Madau 2005)111 Haardt & Madau (2005) refers to an unpublished updated version of Haardt & Madau (1996), specified in cloudy (Ferland et al., 1998) as “table HM05.”. Cooling is calculated via collisional ionization (Abel et al., 1997), radiative recombination (Black & H., 1981; Verner & Ferland, 1996), photoionization, bremsstrahlung, and H, He, and metal line cooling (Cen & Renyue, 1992; Shen et al., 2010b). To calculate the non-equilibrium abundance of molecular hydrogen, we adopt the model used by Christensen et al. (2012). Local, non-equilibrium H2 abundances are determined by the rates of dust grain H2 formation, gas-phase H2 formation, collisional dissociation, and photodissociation by Lyman-Werner radiation from nearby stellar particles. HI and H2 are both shielded by dust and H2 is additionally self-shielded. The Lyman-Werner flux used for calculating the rate of H2 photodissociation and heating for a given gas particle is approximated based on the cosmic ultraviolet background and the average Lyman-Werner flux from nearby star particles. In calculating the stellar Lyman-Werner radiation, the particle gravity tree is used to select proximate gas particles, as described in Christensen et al. (2012). This method avoids a full radiative-transfer calculation while still approximating the spatial and temporal variation in the Lyman-Werner flux. When determining the rates of H2 photodissociation and heating, we further include a sub-grid model for dust and self-shielding. This same dust shielding model is also used in determining the rates of photoionization and heating of HI by the cosmic ultraviolet background. The sub-grid shielding models and the different techniques we explored to estimate column length are described in detail in §2.3. Star formation rates are determined stochastically based on local gas properties, described for each model in section 2.4.

The simulations use the blastwave supernova (SN) feedback model described in Stinson et al. (2006), in which each SN releases $10^{51}$ ergs of energy into the surrounding gas. For SNe II, radiative cooling is disabled within the blastwave radius during the theoretical lifetime of the momentum-conserving phase of the supernova remnant. This sub-grid recipe prevents artificial cooling and allows affected gas to naturally rise from the disk without the use of a momentum kick. SNe I and II also inject metals into the ISM. These simulations assume a Kroupa (2001) initial mass function, SNe rates from Raiteri (1996), and metal yields from Woosley & Weaver (1995) and Thielemann (1986). Metals are also returned to the ISM by stellar winds using mass loss rates from Weidemann (1987), and they are dispersed using a shear-dependent subgrid turbulent mixing model (Shen et al., 2010a) with a metal diffusion coefficient of 0.03.

2.3 Shielding recipes

Shielding is responsible for the preservation of H2 in dense gas as both dust and self-shielding protect the gas from photodissociation. We examine the role that dust shielding may play in promoting star formation by preventing photoheating. As part of these efforts, we improve our shielding models through the implementation of a new model for calculating column length based on Safranek-Shrader et al. (2017).

Shielding has significant effects on scales below the resolution limit of our simulations. We account for these effects with a local, sub-grid approximation of shielding, following the work of Draine & Bertoldi (1996), Glover & Mac Low (2007), and Gnedin et al. (2009). The dust shielding $S_{d}$ and self-shielding $S_{H_{2}}$ functions were implemented as in Christensen et al. (2012): We calculated the dust shielding as

[TABLE]

where N${}_{\text{HI}}$ and N${}_{\text{H2}}$ are the column densities of H${}_{\text{I}}$ and H2 respectively, $Z$ is the metallicity, and $\sigma_{d_{\text{eff}}}$ is a parameter representing the effective attenuation cross section, tuned to be $2\times 10^{-21}$ cm2. The H2 self shielding is calculated as

[TABLE]

where $\omega_{H_{2}}$ is an adjustable parameter tuned to be 0.2 and $x=N_{H2}/(5\times 10^{14})$ cm2.

The column densities are the product of the column length associated with each gas particle and the number densities of HI and H2. The molecular hydrogen recipe introduced in Christensen et al. (2012) assumed the column length to be equal to the particle smoothing length, $h$, a recipe which we will refer to as Ssmooth. However, this recipe has the potential to be susceptible to resolution effects. More recently, Safranek-Shrader et al. (2017) evaluated several common approximations for shielding length—including the Sobolev approximation, the Jeans length, a length based on the local density and its gradient, a power-law approximation, and a single-cell approximation—and compared them to a detailed ray-tracing solution to the radiative transfer problem. They found that a temperature-capped Jeans length performed well at matching the effective visual extinction and performed better than other local models at calculating mass-weighted abundances of H2 and CO.

In the S model, we consider the shielding length to be equal to the particle Jeans length using the empirically-chosen temperature ceiling of 40 K from Safranek-Shrader et al. (2017):

[TABLE]

where T is the minimum of the particle’s temperature and 40 K, $m_{H}$ is the mass of atomic hydrogen, and $\rho$ is the gas density. We also investigated the use of a Jeans length approximation for shielding length without a temperature ceiling (the SJeans model).

In order to assess the differences among our shielding recipes, we examined the resulting column lengths of gas particles in isolated disk galaxies generated using the Ssmooth, S, and SJeans models. We found that the SJeans model resulted in significantly longer column lengths than the Ssmooth or S models (Figure 1). In that model, high-temperature, low-density gas is modeled as having a long column length, increasing the column densities and dust shielding beyond what would be expected in reality. This supports the findings of Safranek-Shrader et al. (2017), who concluded that H2 and CO abundances were significantly over-predicted by a Jeans length shielding model without a temperature cap. The temperature-capped Jeans length, however, led to column lengths similar to the smoothing lengths in the Ssmooth simulation.

The extent of the impact of the different column length recipes is further analyzed in Figure 2, which compares column lengths as a function of density. Although the column lengths in the SJeans simulation are significantly longer on average than the column lengths in the other two shielding simulations, the effect is less significant when looking only at gas particles eligible for star formation, indicated by darker shades in the diagram. The S model creates the shortest column lengths at high (i.e. potentially star-forming) densities. We also note that since almost all gas has temperatures greater than 40 K, the S column lengths are generally simply proportional to $\rho^{1/2}$. Additionally, for gas particles with equal masses and T > 40 K, the column lengths for S and Ssmooth will be correlated with each other.

To further test our shielding model, we examined the HI+H2 surface densities at which hydrogen transitioned from atomic to molecular for simulations using the S model at sub-solar metallicity (SLM) and with a factor of ten lower mass resolution (SLR) (Figure 3). Surface densities were calculated per-particle as $\Sigma_{HI+H2}=\rho(N_{HI}+2N_{H2})L_{shield}$. As metallicity decreased, the transition moved to slightly higher column densities since dust shielding, and therefore the rate of H2 formation, decreased. While the metallicity of the initial conditions differed by a factor of 10, the shift in the transition was roughly a factor of 2. The reduced shift was because of a combination of factors. First, self-shielding, which is not metallicity-dependent, is significant at high densities, enabling the transition from HI to H2 to occur at lower surface densities than dust shielding alone would predict. Secondly, although the average metallicity was initially $z=0.1\;Z_{\odot}$ for the low-metallicity simulation, we found that the difference in metallicities between the two simulations had decreased after they were integrated, particularly for the molecular gas. The average ISM metallicities of each of the simulations after being integrated were $\langle Z_{ISM}\rangle=1.3\;Z_{\odot}$ for S, $\langle Z_{ISM}\rangle=.73\;Z_{\odot}$ for SLM, and $\langle Z_{ISM}\rangle=1.0\;Z_{\odot}$ for SLR.

The low-resolution simulation was generated from a similar set of initial conditions as the other isolated, Milky Way-like simulations. The only differences were a factor of 10 lower mass resolution and $10^{1/3}$ lower spatial resolution. Specifically, the low-resolution simulation had $10^{5}$ gas particles with initial masses $1.4\times 10^{6}$ M⊙, instead of $10^{6}$ gas particles with initial masses $1.4\times 10^{5}$ M⊙. Decreasing the mass resolution of the simulation did not significantly affect the surface density at which the transition occurred, indicating that that the transition is consistent over this range of mass resolution.

2.4 Star formation models

We compared three star formation models to determine the impact of tying star formation to H2, shielded gas, or cold gas (Table 1). Our reference model (H2) was the H2-based star formation model outlined in Christensen et al. (2012). We compare that model to one in which the star formation probability depends on the amount of dust shielding (S). We also examine the case in which the star formation is limited to cold gas (TC), as the motivation for shielding-dependent star formation is that shielding enables gas to cool to low temperatures.

The probability of a sufficiently cold and dense gas particle spawning a star particle is adopted from Stinson et al. (2006):

[TABLE]

where $m_{gas}$ and $m_{star}$ are the masses of the gas particle and the potential star particle respectively, $t_{form}$ is the local dynamical time, and $c^{\ast}$ is a star-forming efficiency factor, whose functionality depends on the model, as described below.

In the H2 recipe, star-forming efficiency was linked to the local abundance of molecular hydrogen:

[TABLE]

where $X_{H_{2}}$ is the number fraction of H2, and $c_{0}^{\ast}$ is an adjustable parameter tuned to 0.1 for the isolated galaxies and 0.01 for the cosmological simulations, values chosen to produce reasonable Kennicutt-Schmidt relations for the H2 model. In the shielding recipes, we modified the efficiency factor to be dependent on dust shielding:

[TABLE]

where $S_{d}$, as defined in equation 1, is the fractional reduction in the intensity (so $S_{d}$ = 1 corresponds to no shielding by dust). We used the same values of $c_{0}^{\ast}$ as in the H2 recipe. In the temperature-based recipe (TC), only those particles passing a more stringent temperature cut were allowed to form stars, and we used a constant value of $c^{\ast}=c_{0}^{\ast}$. We only tested the TC recipe on the isolated galaxies and used the same value of $c_{0}^{\ast}=0.1$ as in the H2 and S cases for those galaxies.

All models required gas particles to meet temperature and density criteria ($\rho\geq\rho_{min}$ and $T\leq T_{max}$) to be eligible for star formation. The same density threshold was used across all models: $\rho_{min}\geq 0.1$ amu cm*-3*, a value chosen to be low enough to have no significant effect on star formation. For H2 and shielding-based star formation models, $T_{max}$ was set to 15000 K, a temperature cut high enough to have no significant effect on star formation. For TC, $T_{max}$ was lowered to 1000 K, a temperature chosen to be similar to the upper temperature at which particles in the simulation contained significant amounts of H2. This temperature ceiling must be higher than the observed temperature of molecular clouds to reflect the limited resolution of the simulations.

In order to test the accuracy of our models, we compared the HI+H2 surface densities at which hydrogen transitioned from atomic to molecular for isolated disk simulations at solar metallicity. The surface densities are plotted in Figure 4 and compared to Milky Way observational data from Gillmon et al. (2006) and Wolfire et al. (2008) determined from absorption spectroscopy. In addition to showing the results of the three different star formation models, we also include shielding-based star formation using the previously discussed Ssmooth and SJeans shielding models. In all other simulations, shielding was calculated using the S model.

In order to ensure an accurate comparison to observations, the surface densities and H2 fractions were calculated from mock observations generated by post-processing the galaxy to create velocity cubes of HI and H2. The cubes were produced by using the smoothing kernel to spatially distribute the gas particles, and then calculating the expected amount of 21-cm emissions and binning it in x, y, and velocity space. Data was generated in 128 velocity channels with 1.3 km s*-1* bins and a 5 pc pixel size.

All models were able to successfully reproduce the observed transition from atomic to molecular hydrogen, indicating that all shielding models produced reasonable behavior. The different star formation criteria allowed for different amounts of H2, as will be discussed further in §3. The H2 and S models produced the greatest amounts of H2-rich gas, as in the other models, gas particles frequently formed stars before a large amount of H2 could accumulate.

Linking star formation to molecular or shielded gas relegates star formation to generally dense and cold gas for which the Jeans mass may no longer be resolved. Robertson & Kravtsov (2008), therefore, advocated the adoption of a pressure floor to prevent numerical Jeans fragmentation of dense, cold gas. However, the use of a pressure floor has also been rejected as artificial or unnecessary. For example, Christensen et al. (2012) found that the addition of a Jeans pressure floor did not significantly affect the number or distribution of stars formed for simulations of this resolution, as gas particles for which the Jeans mass and length are unresolved necessarily meet the star formation criteria. Hopkins et al. (2017) similarly found artificial pressure floors in SPH simulations to be unnecessary and unphysical. We tested the impact of adding a Jeans pressure floor to the S shielding model and found that, due to the high resolution of the simulation, it did not create any significant differences in star formation history, gas density, or gas distribution. The average density at which gas formed into stars was not affected by the presence of a Jeans floor, as shown in Figure 5. All models discussed in this paper therefore do not have an artificial pressure floor.

3 Results

We analyzed the effects of including shielding-based star formation and compared it to the more traditional star formation recipes of a temperature-capped model and a H2 abundance-based model, as shown in §3.1. We used solar and sub-solar metallicity isolated disk simulations to examine the properties of the star-forming gas for each of the three models. By examining the temperatures and densities at which gas formed into stars, we were able to assess how the different star formation recipes affected the environments within which star formation occurred.

In §3.2, we further analyzed the effect of using H2 vs shielding-based star formation on simulations of more realistic galaxies by comparing the evolution of two cosmological dwarf galaxy simulations. These simulations were chosen to highlight the effects of the models in low-metallicity environments. We focused our analysis on the star formation within the simulations, and compared the star formation histories, Kennicutt-Schmidt relation, and stellar and gas density profiles for the simulations. Due to the computational expense, we limited our analysis of cosmological simulations to the H2 and S models.

3.1 Star-Forming Gas in Isolated Simulations

Figure 6 illustrates the ISM for solar and sub-solar isolated disk galaxies simulated with the TC, H2, and S star formation models. The star formation criteria has a clear effect on the structure of the gas. The H2 and S simulations led to formation of denser gas than the TC model. This effect is especially pronounced in the H2 model at solar metallicity, in which the greater amount of dense and cold gas results in the formation of more pronounced spiral arms. For all star formation models, lower metallicity gas resulted in more compact disks, because the reduced cooling efficiency prevented gas from reaching high densities, and low-metallicity gas took longer to collapse into a disk for a given initial radius.

Narrowing our focus to the star-forming gas, Figure 7 shows the temperatures and densities at which star formation occurred for all star particles in the solar and sub-solar metallicity simulations. In the H2 simulation, no star formation occurred at densities below 1.0 atomic mass unit cm*-3*, and little occurred at temperatures above 1000 K, despite the lack of an explicit temperature cap. In the S simulation, however, star formation occurred at densities as low as 0.1 atomic mass units cm*-3* and temperatures as high as 15000 K. While explicitly limiting star formation to gas with temperatures below 1000 K, the TC model also allows for star formation at densities lower than those in the H2 model and similar to the S model. When comparing across metallicity, star formation occurred at slightly higher densities in the sub-solar metallicity simulations than the solar metallicity simulations, regardless of the star formation criteria. In the sub-solar metallicity gas, higher densities were required for the gas to cool, shield, or to maintain H2.

The S simulation resulted in star formation in gas with densities below $1$ amu cm*-3* and temperatures above $10^{3}$ K. Much of this material is low-density gas in thermal equilibrium that is continuing to cool. Those gas particles that are both higher density and hot include gas that, while clearly able to form stars, is not currently in thermal equilibrium because it recently experienced an injection of energy from nearby supernovae. Because of the use of a delayed-cooling feedback model, these particles may currently or until recently have had their cooling disabled.

Although we would not expect significant amounts of star formation at such low densities and high temperatures in reality, the densities and temperatures of gas particles in the simulation should be interpreted as average values for the entire gas cloud. Therefore, given the limited resolution of the simulations, these low-density, high-temperature star forming gas particles represent clouds of gas that may contain much denser regions. Finally, the star formation in the shielding simulation was still concentrated at higher densities and lower temperatures than the temperature capped simulation, indicating that although the S model does allow for some low-density star formation, most stars still form out of cold, dense gas as would be expected from observations.

We further compare the densities at which star formation occurred in Figure 8. As seen in Figure 7, both star formation model and metallicity had a strong effect on the density of the star forming gas. As predicted, star formation occurred at higher densities in the sub-solar metallicity simulations than the solar metallicity simulations for all models. The offset between the peaks of the histograms for the solar and sub-solar simulations was similar regardless of the star formation recipe. As previously seen, the H2 model required the highest gas densities for star formation, followed by the S model, and then the TC model.

A clear difference between the shielding-based and H2-based star formation models is the reduced density of the star forming gas in the former case. This difference happens in spite of the close causal link between H2 and dust shielding, and is due to the discrepancy between the chemical equilibrium timescale and the gas collapse timescale in some environments. Following the derivation in Krumholz (2012), we can articulate the difference in time scales as follows. Molecular hydrogen formation occurs on the surface of dust grains, at a rate per H atom of $\bar{n}_{H}RC$ where $C=\langle n_{H}^{2}\rangle/\bar{n}_{H}^{2}$ is a clumping factor and $R\approx 3\times 10^{-17}Z$ cm3 s*-1* is the rate coefficient for H2 formation on grain surfaces (Wolfire et al., 2008). Then the timescale for conversion of hydrogen from atomic to molecular is:

[TABLE]

The timescale over which gravitationally unstable gas will collapse and begin to form stars is the free-fall time:

[TABLE]

where $\mu_{H}\approx 1.4$ is the mean mass per H nucleus in units of the hydrogen mass $m_{H}$. The ratio of these timescales is

[TABLE]

where $n_{o}=\bar{n}_{H}/$1 cm*-3* (Krumholz, 2012) and C = 10 for these simulations. At low metallicities, the chemical equilibrium timescale can be longer than the free-fall timescale. For example, for metallicity equal to 0.1 Z*⊙, $t_{chem}>t_{ff}$ for number densities less than 580 cm-3*. This difference in time scales suggests that, assuming H2 is not required for star formation, a cloud of gas could begin to form stars before a substantial amount of H2 is formed. In our implementation, this resulted in gas forming stars at lower density in the S model than in the H2 model. In the latter situation, the gas had to form sufficient H2 for star formation to proceed, thereby causing the gas to collapse to higher densities first. Therefore, the low-density star formation in the shielding model suggests one possible model for the small amounts of star formation that do occur in environments that have low metallicities and low average surface densities, such as dwarf galaxies and the outer disks of spiral galaxies.

3.2 Star Formation in Cosmological Dwarfs

We examine the impact of the shielding-based star formation prescription on the evolution of a dwarf galaxy. These simulations were chosen to highlight any differences that might occur between shielding- and H2-based star formation in the low-metallicity, low-surface density environments of dwarf galaxies. In particular, we investigated the effect of the star formation models on the history and distribution of star formation as well as the resulting galaxy morphology.

To assess the differences in star formation among the various models, we first compared the star formation in our simulated galaxies to the observed redshift zero resolved Kennicutt-Schmidt relation (Figure 9). To replicate the observations, we calculated surface densities of HI and H2 from mock THINGS (Walter et al., 2008) observations generated at the resolution and sensitivity for a 5 Mpc-distant dwarf galaxy. We generated data in 128 velocity channels with 1.3 km s*-1* bins and a 1.5-arcsec pixel size for the galaxy oriented at a 45°angle. We then spatially smoothed the data cubes using a Gaussian beam with a full width at half-maximum of 10 arcsec $\times$ 10 arcsec. Finally, we made a sensitivity cut and discarded all emission from cells below a $2\sigma$ noise limit, in which $\sigma=$ 0.65 mJy beam*-1*. Star formation rates were calculated using the simulated far UV and 24 $\mu$m emission generated by the radiative transfer code, sunrise (Jonsson, 2006). Surface densities of gas and star formation were then calculated for 7502 pc2 grid squares.

Both the S and H2 model produce very similar results to each other and to the observed trend for local galaxies (Bigiel et al., 2010). Both galaxies have elevated surface densities and star formation rates at their center, resulting in a data point above the observations. The scatter in the simulated data for locations beyond $r_{25}$ is also somewhat higher than in the observed case. Otherwise, however, the data from the simulated galaxies follows the same trend as observed galaxies. While the H2 model results in slightly more high-surface density data points than the S model, both models produce very similar relations to each other. This similarity is in spite of the fact that the S model results in lower density gas particles forming stars. Apparently, when the surface densities are averaged over a 7502 pc2 grid, the comparatively high densities of the star forming gas in the simulation with the H2 model were diluted by the surrounding interstellar media. Similar results were found in Christensen et al. (2012), where a density threshold-based star formation law produced a similar Kennicutt-Schmidt relation to the H2-based star formation law. We therefore conclude that despite the differences between the shielding and H2-based star formation models, both are consistent with the resolved Kennicutt-Schmidt relation.

One environment where shielding-based star formation might occur is in the low-metallicity outskirts of galaxies (Ferguson et al., 1998; Bigiel et al., 2010; Goddard et al., 2010) In these regions, star formation is observed to occur despite very low levels of H2, and it is possible that the initial gravitational collapse for star formation might occur prior to H2 formation. We therefore compare the ISM, stellar, and star formation rate surface densities as a function of radius (Figure 10). Both models produce similar gas profiles, though the H2 model does produce a higher surface density of H2. Likewise, the stellar profiles are almost identical. However, we do find that the star formation rate surface density in the S model is slightly higher than in the H2 model, particularly in the outer disk, despite the lower H2 surface densities.

Although the average gas surface densities were similar between the simulations produced with the H2 model and the S model, a detailed analysis of the particle densities reveals differences (Figure  11). We found that the H2 simulation produced more numerous and higher density clumps of gas than the S simulation. These clumps are evident in the figure as spatially localized spikes in the particle densities. This difference in gas densities follows from the tendency for star formation to occur in lower density gas in the S simulation. As a result, the galaxy was less able to maintain high density clumps of material. The similarity in the average surface densities between the models, which was also apparent in the K-S law (Figure 9) and gas profile (Figure 10), camouflages these small scale differences.

Finally, we compare the star formation across the entire history of the galaxy. The total stellar masses formed over the course of the simulations were $3.6\times 10^{9}$ M⊙ for the H2 model and $4.1\times 10^{9}$ M⊙ for the S model, indicating that the star formation recipe had a significant effect on the total mass of stars produced. This difference in stellar mass could have been due to either or both of the following possible causes. First, star formation may have been allowed in a wider range of environments in the S simulation. From Figures 7 and 8, it is clear that the S model expands the range of gas particles capable of star formation, supporting this possibility. Secondly, changes to the star formation environments caused by the different models could have had other indirect effects; in particular, the efficiency of the feedback recipe may have been affected. This possibility is supported by previous research with similar simulations. For example, Governato et al. (2010) and Christensen et al. (2012) showed that the blastwave feedback recipe is most effective at regulating star formation when star formation is confined to high-density clumps. It is likely that both explanations play a role here in the total mass formed. It is also important to understand that the amount of feedback energy injected per SNe (dESN) is a tunable parameter that must be set for each combination of star formation parameters. While we used the same value of dESN for both simulations to create a fair comparison, when using either star formation model in other research, dESN would be set to ensure appropriate z = 0 stellar masses.

Looking in more detail at the star formation histories (Figure 12), it is clear that the S model begins star formation earlier, resulting in greater amounts of star formation within the first 2 Gyrs. During this early time period, the metallicity of the ISM is extremely low. Low metallicity not only increases the surface densities at which dust shielding occurs, it also decreases the rate of H2 formation on dust grains. As a result, while low metallicity limits star formation in both the S and H2 models, it has a stronger effect with the H2 model.

After 2 Gyrs, the star formation histories have a similar shape, though the star formation rate is consistently slightly higher in the S simulation. Both the H2 and S recipe produce bursty star formation histories at high redshift and smoother star formation histories at low redshift. Bursty star formation histories with periods of 150-200 Myr are characteristic of present-day dwarf galaxies (Weisz et al., 2012). To resolve bursts on these timescales, we examine the star formation history over a 1 Gyr time period (Figure 12). Bursts of star formation lasting 50-100 Myr appear in both the H2 and S simulations during the first Gyr but are not present in either during the latter part of the simulation, even when examining 1 kpc2 subsections of the disk. Continually bursty star formation histories for dwarf galaxies are characteristic of simulations with efficient supernova feedback and high star formation density thresholds (Governato et al., 2010; Hopkins et al., 2014; Shen et al., 2014; Sparre et al., 2017). Since both the H2 and S simulations have these qualities, it is somewhat surprising that neither maintain bursty star formation across the entire history of the galaxies. Nevertheless, the similar shapes of the star formation histories indicate that both the H2 and S star formation models are operating in a similar fashion.

4 Discussion and Conclusions

In this study we developed and analyzed an implementation for dust shielding-based star formation in galaxy evolution simulations. We examined the effect of tying star formation probability to the amount of dust shielding, rather than the fraction of molecular hydrogen, finding that the former resulted in star formation in lower-density gas while still being consistent with the Kennicutt-Schmidt relation. We summarize our implementation and these effects below.

The shielding-based star formation model consisted of a stochastic sub-grid model in which the probability of a gas particle forming a star particle was a function of the amount of dust shielding. This dependency on dust shielding replicates the theoretical model in which only shielded gas is able to achieve the low temperatures necessary for gravitational collapse. As part of the model development, different prescriptions were used for calculating the column length associated with a given gas particle. In agreement with Safranek-Shrader et al. (2017), we found that a prescription in which the column length is the Jeans length calculated with a 40 K temperature cap produced the best agreement with observations. With this prescription, we were able to replicate the observed transition from HI to H2 as a function of the HI + H2 surface densities calculated for the disk during post-processing.

We compared our shielding-based star formation model to a molecular hydrogen-based star formation model and a temperature ceiling model using simulations of isolated disk galaxies. The primary differences between the models were the densities and temperatures of the star forming gas. The shielding-based model resulted in star formation at higher temperatures and lower densities than the H2-based model, while still requiring higher densities on average than a temperature ceiling recipe. One consequence of the different star forming gas densities was the presence of denser gas in the simulation using the H2 model, which also lead to more pronounced spiral structure.

The denser star forming gas in the H2 model can be explained by the time delay introduced by requiring the gas to form H2 prior to star formation. During the time that the gas is forming H2, it also may gravitationally collapse and reach higher densities before spawning a star particle. In models that do not rely on an explicit link between H2 and star formation, the discrepancy between free-fall timescales and molecular formation timescales allows star formation to begin earlier while the gas particle is at a lower density.

We further compared the shielding and H2 star formation models by examining the evolution of cosmological simulations of dwarf galaxies. Both the shielding and H2 models produced similar fits to the observed Kennicutt-Schmidt relation when resolved on 750 pc scales. At this resolution, differences in the densities of the star forming gas were washed out, and both recipes provide realistic results. Despite similar gas surface density profiles, the star formation rate in the outer disk was slightly higher for the shielding simulation than the H2 simulation. This difference in star formation at large radius is likely the result of shielding-based star formation being able occur at lower densities. However, the resulting stellar profiles were essentially the same for both models. Likewise, both star formation models produced similarly shaped star formation histories with similar levels of burstiness. The exception to this similarity is that star formation is able to begin slightly earlier when using the S model, because extremely low metallicity has a weaker effect on the presence of dust shielding than the amount of H2. The star formation rate continued to be higher when using the S model for most of the duration of the simulation.

The strongest difference the star formation models produced in the cosmological simulations was in the total mass of stars formed. The shielding model resulted in approximately ten percent more stellar mass formed over the history of the universe than the H2 model. The additional stellar mass appears to be driven by a higher star formation rate in the outer disk. This difference appears to be caused not by differences in the amount of star formation per cold gas, as seen in the Kennicutt-Schmidt relation, but by the existence of greater amounts of cold gas in the H2 simulation. We attribute this differences to a decrease in the effective efficiency of the supernova feedback model when star formation is allowed to take place in lower density gas. Similar behavior was observed in Christensen et al. (2012), where one of the main impacts of using an H2-based star formation recipe was the enhanced feedback efficiency due to star formation happening in higher density gas. The slightly lower relative efficiency of feedback observed when using the shielding model can be modulated by increasing the amount of energy released per supernovae mass, $dESN$. We therefore stress that while both dust shielding and H2-based star formation models produce realistic behavior (e.g. matching the Kennicutt-Schmidt relation), the tuning of $dESN$ and other star formation or feedback models must be model specific.

While we found that the dust shielding-based model enable star formation in lower density gas than the H2-based model, there were only slight differences in the observable properties of the resulting galaxies. It is likely that these effects would be stronger in simulations that focused on the extremely early Universe or that resolved the internal structure of star forming clouds. On the resolution and time scales of cosmological galaxy formation simulations, though, the different models do not produce markedly different observable characteristics. Therefore, we conclude that both the dust shielding-based and H2-based star formation models are viable and realistic in the context of galaxy formation simulations.

5 Acknowledgments

The authors would like to thank the anonymous referee for insightful comments. LB acknowledges support through the Grinnell College Mentored Advanced Project program. The cosmological simulations were computed at NASA AMES. Analysis of the simulations was undertaken using the Python analysis package for astrophysical N-body and Smooth Particle Hydrodynamics simulations, pynbody (Pontzen et al., 2013).