KINEMATIC EVOLUTION OF SIMULATED STAR-FORMING GALAXIES

Over the last ∼8 billion years since a redshift of one, the population of star-forming galaxies of Milky Way mass has settled kinematically into flat, rotationally supported, disk galaxies (Kassin et al. 2007, 2012). In the past, these galaxies had larger integrated gas velocity dispersions (σ_g) and less ordered rotation (V_rot) than they do today (Flores et al. 2006; Kassin et al. 2007, 2012; Vergani et al. 2012). This σ_g is likely due to disordered motions in these galaxies (Covington et al. 2010; Kassin et al. 2012). Over 0.1 < z < 1.2, the median σ_g of star-forming galaxies of Milky Way mass has progressively decreased while V_rot has increased (Kassin et al. 2012). Observations at even higher redshifts (z ∼ 1.5–3), albeit where galaxy samples are much less representative, also find star-forming galaxies with large amounts of disordered motions, as measured through σ_g (e.g., Förster-Schreiber et al. 2009; Law et al. 2009; Wright et al. 2009; Lemoine-Busserolle & Lamareille 2010; Lemoine-Busserolle et al. 2010; Gnerucci et al. 2011).

These findings are an important benchmark for any theory or simulation of disk galaxy formation. Constraints on the internal kinematics of galaxies are significantly more stringent than constraints on, e.g., luminosity, stellar mass, or star formation rate. Internal kinematics directly dictate how the gas in galaxies is arranged, which is heavily influenced by the physical processes involved in galaxy evolution. Furthermore, if theory is able to successfully reproduce observations of kinematic evolution, then it becomes a useful tool to help interpret the observations. In particular, we would like to know the physical processes behind disk settling.

Until now, robust measurements of the internal kinematics of simulated galaxies were lacking because the simulations suffered from low resolution, overproduction of stars, and overly massive bulges. These problems were ameliorated by increased resolution (e.g., Mayer et al. 2001), better treatment of feedback (Stinson et al. 2006; Hopkins et al. 2012; Wise et al. 2012), and modeling of the ultraviolet background (Quinn et al. 1996). Among other things, the current generation of simulations are better able to model star formation, inflows, and outflows than previous generations. This makes them better able to match the rotation velocities, sizes, stellar masses, and baryon fractions of local disk galaxies (e.g., Brooks et al. 2011; Brook et al. 2012; Munshi et al. 2013; Christensen et al. 2013; Vogelsberger et al. 2013), as well as some properties of higher redshift galaxies (e.g., Hirschmann et al. 2013; Zemp et al. 2013; Hopkins et al. 2013). These simulations are the first generation which are able to make quantitative predictions for the internal kinematics of galaxies.

A few studies using these current generation simulations have measured kinematics, and most of them have only measured rotation velocities for local disk galaxies (e.g., Guedes et al. 2011; Stinson et al. 2013; McCarthy et al. 2012; Brook et al. 2012). Fewer simulations have measured rotation velocities at higher redshifts, and fewer still have measured σ_g at any redshift: e.g., Rôskar et al. (2013) and Agertz et al. (2013) for z = 0, Croft et al. (2009) for z = 0 and z = 1, Ceverino et al. (2010) and Anglés-Alcázar et al. (2014) for z = 2, and Robertson & Bullock (2008) which studied gas-rich mergers which should be common at z = 2. Bird et al. (2013) studied the evolution of V_rot/σ_g over the lifetime of a single star-forming galaxy. No studies have yet compared simulated galaxies to the progressive disk settling found at 0.1 < z < 1.2 by Kassin et al. (2012). In this paper, for the first time, we determine whether a suite of fully cosmological simulations of star-forming galaxies undergoes the kinematic disk settling found in observations of real galaxies.

The kinematics of star-forming galaxies in the real universe are measured from nebular emission lines which trace gas heated in star-forming regions. In order to compare hydrodynamic simulations of galaxies most directly with observations, the simulations need to make predictions for the gas in galaxies. Treatment of the interstellar medium is crucial for this. Therefore, we look to simulations of four star-forming galaxies run with the N-Body + smooth particle hydrodynamic code GASOLINE (Wadsley et al. 2004; Stinson et al. 2006) which are described in detail in Christensen et al. (2012, 2013) and Munshi et al. (2013). They are able to match the sizes, stellar masses, bulge masses, and baryon fractions of local star-forming galaxies (Governato et al. 2007; Brook et al. 2012; Munshi et al. 2013; Christensen et al. 2013), and the size evolution of star-forming galaxies since z = 1 (Brooks et al. 2011). The simulated galaxies have z = 0 stellar masses which range from 4.2–4.5 × 10¹⁰ M_☉ for a Kroupa et al. (1993) initial mass function (IMF). They are referred to as galaxies h239, h277, h258, and h285 in the references above.

These simulations resolve high density peaks comparable in size to giant molecular clouds. In addition, they follow the non-equilibrium formation and destruction of H₂, using a gas-phase and a dust- (and hence metallicity-) dependent scheme that traces the Lyman–Werner radiation field and allows for self-shielding by H₂ gas (Gnedin et al. 2009; Christensen et al. 2013). Therefore, the simulations are able to tie star formation to the presence of molecular gas, as is indicated by observations (e.g., Leroy et al. 2008; Bigiel et al. 2008, 2010; Blanc et al. 2009; Schruba et al. 2011). The efficiency of star formation is tied to the H₂ fraction, as described in Christensen et al. (2012). With the inclusion of H₂-based star formation, stars predominantly form at high densities and low temperatures (T <1000 K; see also Krumholz et al. 2011; Kuhlen et al. 2012). Finally, to take into account the effects of the reionization of the universe, an ultraviolet background is implemented at z = 9, following a modified version of the formulation by Haardt & Madau (2001).

Our simulated galaxies were chosen from lower resolution simulations to be resimulated at high resolution based on their z = 0 virial masses and because they span a representative range of halo spin values and accretion histories. The resimulations follow a 50 Mpc co-moving box around the galaxies at lower resolution than the galaxy itself. In this manner, cosmological effects are fully implemented. The resimulations are run from z = 150 to z = 0.

The spline force softening in the high resolution region is 174 pc, and is kept fixed in physical parsecs since z = 10. In the high resolution region, dark matter particles have masses of 1.3 × 10⁵M_☉, gas particles start at 2.7 × 10⁴M_☉, and star particles are born with 30% of the mass of their parent gas particle, which corresponds to a maximum initial mass of 8100 M_☉. Each simulated galaxy has ∼5 million dark matter particles within its virial radius at z = 0 and more than 14 million total particles (dark matter, gas, and stars). Our simulations have the same mass resolution as the Eris simulation (Guedes et al. 2011), but include metal line cooling (Shen et al. 2010) and the physics of molecular hydrogen (Christensen et al. 2012). The simulations are run for a Wilkinson Microwave Anisotropy Probe 3 yr Cosmology: Ω_m = 0.24, Λ = 0.76, H_o = 73 km s⁻¹,σ₈ = 0.77 (Spergel et al. 2007).

2.1. Warm and Cold Gas in the Simulations

For each of the four simulated galaxies, we measure σ_g and V_rot at discrete redshifts sampling the range of observations in Kassin et al. (2012). We measure the intrinsic values of these quantities which are unaffected by observational effects such as seeing, slit width, and pixel scale. The quantity σ_g is measured as the average integrated velocity dispersion of the gas in the direction perpendicular to the disk, similar to observations. To measure σ_g, we step across the galaxies in a face-on position in 0.5 kpc radial bins and measure the line-of-sight velocity dispersion in each. The value of σ_g is taken as the mean value of the line-of-sight dispersions among the bins. The rotation velocity V_rot is taken as the maximum line-of-sight rotation velocity of the gas. To measure V_rot, we step across the galaxies in an edge-on position in 0.5 kpc radial bins and measure the line-of-sight rotation velocity in each. We use the actual gas particle velocities rather than circular velocities since this provides the most direct comparison with observations. From these velocity measurements, we create a rotation curve (rotation velocity versus radius). The maximum value of the rotation curve is adopted as V_rot, similar to the observations which use the rotation velocity on the flat part of the rotation curve.

Our measurements are of intrinsic quantities and therefore are different from the mock observations of, e.g., Covington et al. (2010) which take into account myriad observational effects. Our goal is to determine the intrinsic kinematic evolution in the simulations, not to investigate observational effects as in Covington et al. (2010).

Median values of σ_g and V_rot at discrete redshifts for the warm and cold gas in our four simulated galaxies are shown in the top panels of Figure 1. Values of these quantities for the individual simulated galaxies are shown in the bottom panels. In the top panels, the median values are compared to those for an observed mass-limited sample of 270 star-forming galaxies from Kassin et al. (2012). Qualitatively, the simulated galaxies follow the same trends as the observations: they increase in σ_g and decrease in V_rot with increasing redshift over 0.1 < z < 1.2. In other words, both the simulations and the observations decrease in σ_g and increase in V_rot with time over the last ∼8 billion years, demonstrating that the simulated galaxies undergo the disk settling found in observations. Similarly, Bird et al. (2013) find that the ratio of ordered to disordered motions (V_rot/σ_g) decreases with time for a similar mass galaxy from the Eris simulation. Furthermore, the scatter in σ_g and V_rot for the warm/cold gas in the individual simulated galaxies is large, similar to the scatter in the observations (Figure 5 in Kassin et al. 2012).

Zoom In Zoom Out Reset image size

As expected, the normalizations of the simulated and observed relations differ. The median values of V_rot for the simulated galaxies are greater since they are more massive on average than the observed galaxies, and V_rot generally scales with stellar mass for disk galaxies. The masses of the simulated galaxies range from 4.2–4.5 × 10¹⁰ M_☉ versus 6.3–50.1 × 10⁹ M_☉ for the observations. (The simulations and observations adopt Kroupa et al. 1993 and Chabrier 2003 IMFs, respectively, which result in consistent stellar masses.) The normalization of the σ_g versus z relations for the warm and cold gas are higher and lower than the observations, respectively. In the simulations, the warm gas is affected by feedback from supernovae, which results in higher σ_g and leads to hotter temperatures compared to observed H ii regions. By definition, the cold gas is cooler than observed H ii regions, which results in lower σ_g. In addition, we note that the median stellar mass of the mass-limited observational sample decreases with decreasing redshift (Figure 4 in Kassin et al. 2012). Since the observed evolution is to increasing V_rot and decreasing σ_g with decreasing redshift, this changing median mass makes these trends lower limits to the intrinsic evolution. The magnitude of this effect is unclear and we defer a more detailed analysis of it to future comparisons with significantly larger samples of simulated galaxies.

3.1. Quantitative Trends of σ_g and V_rot with Redshift

Linear relations are fit to the median points in Figure 1 by performing least-squares fits which take into account the errors in σ_g and V_rot (errors in redshift are negligible). When fitting, to avoid covariances, we zero-point the medians near the middle of the samples such that they vary around ∼zero. We obtain the following fits for σ_g for the warm gas, cold gas, and observations, respectively:

$$\begin{equation} \sigma _g - 50 = (21.5 \pm 6.2) (z - 0.6) - (4.0 \pm 1.8), \end{equation} \tag{ 1 }$$

$$\begin{equation} \sigma _g - 20 = (11.7 \pm 7.4) (z - 0.6) - (2.1 \pm 1.8), \end{equation} \tag{ 2 }$$

$$\begin{equation} \sigma _g - 35 = (26.0 \pm 5.7) (z - 0.6) - (3.1 \pm 1.4). \end{equation} \tag{ 3 }$$

These fits have χ² values of 0.8, 1.3, and 4.7, respectively. As for the observations, the median σ_g of both the warm and cold gas grows with increasing redshift to z = 1.2 (i.e., decreases with time), although the relation for the cold gas is consistent with no evolution. The slopes of the relations for the warm gas and observations are consistent within uncertainties. As mentioned above, the normalization of the relations for the warm and cold gas are higher and lower than the observations, respectively.

We obtain the following fits for V_rot for the warm gas, cold gas, and observations, respectively:

$$\begin{equation} V_{\rm rot} - 175 = (-56.9 \pm 22.0) (z - 0.6) - (7.9 \pm 6.2), \end{equation} \tag{ 4 }$$

$$\begin{equation} V_{\rm rot} - 200 = (-20.7 \pm 15.1) (z - 0.6) - (4.9 \pm 4.3), \end{equation} \tag{ 5 }$$

$$\begin{equation} V_{\rm rot} - 100 = (-48 \pm 11.5) (z - 0.6) - (7.8 \pm 3.0). \end{equation} \tag{ 6 }$$

These fits have χ² values of 2.7, 1.5, and 9.2, respectively. As for the observations, the median V_rot of both the warm and cold gas decreases with increasing redshift to z = 1.2 (i.e., increases with time). The simulated galaxies have a higher normalization (i.e., have faster median rotation velocities) than the real galaxies, as expected because they are on average more massive than the observed galaxy sample, as discussed above.

3.2. Kinematics and Morphology

As is the case for real galaxies over 0.1 < z < 1.2 (e.g., Flores et al. 2006; Kassin et al. 2007, 2012; Yang et al. 2008), the kinematics of simulated galaxies are reflected in their morphologies. Images of two of our four simulated galaxies are shown in Figure 2 at discrete redshifts over the redshift range in Figure 1. Both galaxies have more disturbed morphologies at z = 1.0, and become progressively less disturbed with decreasing redshift. At z = 0 they appear as ordered disk galaxies with little disturbances or peculiarities. This morphological transformation is consistent with the kinematic settling they undergo (i.e., decrease in σ_g and increase in V_rot with time).

Zoom In Zoom Out Reset image size

The images in Figure 2 are produced by taking into account the effects of radiative transfer and dust using a Monte Carlo ray tracing program designed to pair with hydrodynamic simulations (Sunrise; Jonsson 2006; Jonsson et al. 2010). In short, Sunrise creates a spectral energy distribution for each star particle in the simulation based on its age and metallicity, using the Starburst99 stellar population modeling software (Leitherer et al. 1999). The metallicities of the gas particles determine how light from the galaxy is attenuated by dust, and a constant dust-to-gas ratio of 0.4 is assumed. The resulting images are then convolved with the Sloan Digital Sky Survey g-, r-, and i-band filters to produce those in Figure 2.

For the first time, we measure the evolution of the kinematics of a suite of simulated star-forming galaxies over the last ∼8 billion years, over nearly half of the age of the universe. Specifically, we measure the evolution of disordered motions (σ_g) and ordered rotation (V_rot) of gas in a suite of four simulated galaxies. Measurements are compared with recent observations which show the progressive settling of gas in disk galaxies over 0.1 < z < 1.2 from disordered systems into the ordered disk galaxies common today (Kassin et al. 2012). We find that the simulated galaxies follow the same trends as the observations: they progressively decrease in disordered motions (σ_g) and increase in ordered rotation (V_rot) with time. The scatter in σ_g and V_rot among the four simulated galaxies is also similar to that found for the observations. Differences in normalization between the observations and simulations can be attributed to differences in the average stellar mass of the galaxy samples and gas temperatures/regions probed. A larger sample of simulated galaxies by at least an order of magnitude is needed to confirm our findings. However, it is encouraging that with the simulations we have at hand, we find similar trends and scatter as the observations.

Reproducing these observations is an important benchmark for any theory or simulation of disk galaxy formation. The observations trace the internal kinematics of star-forming galaxies over a significant period of time. Internal kinematics directly dictate how the gas in galaxies is arranged, and are strongly affected by physical processes such as feedback from star formation, major/minor mergers, and smooth accretion of baryons. Our next step in a future work is to place constraints on the processes which have the most direct effect on disk settling.

F.G. acknowledges support from NSF grant AST-0607819. Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center. A.B. thanks Jay Gallagher for helpful conversations.