MAGNETIC FIELDS AND GALACTIC STAR FORMATION RATES

Understanding and thus predicting the star formation rate (SFR) that results from a galactic disk of given local properties, such as gas mass surface density, ${{{\Sigma }}_{g}}$, and orbital timescale, is a necessary foundation on which to build theories of galaxy evolution. Global and kiloparsec-scale correlations between star formation activity, gas content, and galactic dynamical properties have been observed (Kennicutt & Evans 2012). In molecular-rich regions of normal disk galaxies, overall SFRs are relatively slow and inefficient, i.e., a small fraction, ${{\epsilon }_{{\rm orb}}}\simeq 0.04$, of total gas is converted to stars every local galactic orbital time (Suwannajak et al. 2014). Most star formation is concentrated within ∼1–10 pc-scale regions within giant molecular clouds (GMCs), but even here star formation appears to be inefficient, with similarly small fractions, ${{\epsilon }_{{\rm ff}}}\simeq 0.01-0.05$, of total gas forming stars every local free-fall time, ${{t}_{{\rm ff}}}\equiv {{(3\pi /[32G\rho ])}^{1/2}}$, where ρ is gas density (Krumholz & Tan 2007; Da Rio et al. 2014). Low efficiencies of star formation may be the result of some combination of turbulence (Padoan et al. 2014), magnetic fields (McKee 1989; Tassis & Mouschovias 2004), and star formation feedback (Krumholz et al. 2014). Here we consider the effects of magnetic fields on the evolution, collapse, and SFR of a turbulent, shearing, kiloparsec-scale region of a galactic disk, extracted from the global galaxy simulation of Tasker & Tan (2009, hereafter TT09).

Van Loo et al. (2013, hereafter Paper I) showed that in the limit of purely hydrodynamic evolution with no star formation feedback, GMCs forming, evolving, and collapsing in this same kiloparsec-scale environment produce too much high density gas and, hence, the SFR is too high. Additional processes are needed to provide extra support to or disruption of the GMCs. One such process may be supported by magnetic fields. Observations show that clouds such as, e.g., Taurus and ρOph, are threaded by strong B-fields (e.g., Heiles 2000; Goldsmith et al. 2008) and that the field is connected on larger scales to galactic fields (Li & Henning 2011). The average solar-neighborhood Galactic magnetic field is $6\pm 2\;\mu {\rm G}$ (Beck 2001). Crutcher et al. (2010) argue for a random distribution of field strengths from 0 to ${{B}_{{\rm max} }}$ with ${{B}_{{\rm max} }}={{B}_{0}}=10\;\mu {\rm G}$ in low density gas (${{n}_{{\rm H}}}\lt 300\;{\rm c}{{{\rm m}}^{-3}}$) and ${{B}_{{\rm max} }}={{B}_{0}}{{({{n}_{{\rm H}}}/300\;{\rm c}{{{\rm m}}^{-3}})}^{0.65}}$ at higher denisties. Dynamo amplification to equipartition values is a natural expectation, which has also been seen in numerical simulations (e.g., Wang & Abel 2009; Pakmor & Springel 2013). The magnitude and direction of the B-field may play an important role in the formation of molecular clouds that are forming from compression of more diffuse atomic gas (Heitsch et al. 2009; van Loo et al. 2010). By potentially stabilizing and supporting GMCs, they may increase cloud lifetimes to the point that other processes such as GMC collisions (Tan 2000) or spiral arm passage (Bonnell et al. 2013) become important.

Here, we extend the model presented in Paper I to include B-fields, varying their strength to investigate the effect on GMC evolution, especially the amount of high-density clump gas and SFR. Section 2 describes the numerical model, Section 3 presents the results, and Section 4 summarizes and discusses their implications.

We start with the same initial conditions as Paper I, i.e., a kiloparsec-sized box at $4.25\;{\rm kpc}$ from the galactic center extracted from the TT09 simulation. We include the same physical processes, i.e., heating/cooling functions developed in Paper I, and adopt a static background potential yielding flat rotation curve of $200\;{\rm km}\;{{{\rm s}}^{-1}}$. As the TT09 simulation is nonmagnetic, no self-consistent initial magnetic field configuration is available. Therefore, for simplicity, we thread the domain with uniform B-field along the shearing direction, i.e., $B={{B}_{0}}{{\hat{1}}_{y}}$. For B₀ we assume 0 (no B-field), 10, or 80 μG. The 10 μG case represents field strengths close to mean Galactic values expected at an inner (∼4 kpc) location. The 80 μG case exceeds observed kiloparsec-scale B-fields by $\sim 1\;$ dex, but allows for stronger fields inside GMCs and their substructures, which are already present in the initial conditions. Using the empirical relation of Crutcher et al. (2010), a field strength of ${{B}_{{\rm max} }}=80\;\mu$G corresponds to a density of ${{n}_{{\rm H}}}\simeq 7000\;{\rm c}{{{\rm m}}^{-3}}$. The mean density of the two most massive GMCs identified in Paper I is $\simeq 470$ and $330\;{\rm c}{{{\rm m}}^{-3}}$, but they contain substructures extending to more than ten times higher densities. Thus our adopted high magnetic field case may be more appropriate for predicting the dynamical evolution of these dense clumps and filaments. The above field strengths can also be compared to the critical values needed to support idealized self-gravitating clouds of a given mass surface density ${{B}_{{\rm crit}}}=21.6\;{\Sigma }/(0.01\;{\rm g}\;{\rm c}{{{\rm m}}^{-2}})$ $\mu {\rm G}=45.1\;{\Sigma }/(100\;{{M}_{\odot }}{\rm p}{{{\rm c}}^{-2}})\;\mu {\rm G}$ (e.g., Mouschovias & Spitzer 1976; McKee & Stars 1999). The TT09 simulation GMCs have typical ${\Sigma }\sim {\rm few}\times 100\;{{M}_{\odot }}{\rm p}{{{\rm c}}^{-2}}$, so only in the strongest field case do we expect significant influence on global GMC dynamics. Indeed, only in this 80 μG case are initial mass-to-flux ratios of the three least massive GMCs below the critical value. Otherwise, the GMCs are supercritical by at least an order of magnitude. A larger number of intermediate B-field strengths and different initial geometries will be explored in a future paper.

To solve the ideal magnetohydrodynamics (MHD) equations we use the Adaptive Mesh Refinement MHD code MG (Van Loo et al. 2006; Falle et al. 2012), rather than Enzo, which was used in Paper I. The basic algorithm is a second-order Godunov scheme with Local Lax–Friedrichs solver and piecewise-linear reconstruction method. Self-gravity is computed using a full approximation multigrid to solve the Poisson equation. We also solve the non-conservative internal energy equation to ensure positive pressures. We use the gas temperature value from internal energy when this is $\lt 1/10$ of total energy and from total energy otherwise. To ensure that the solenoidal constraint is met, divergence cleaning is implemented following Dedner et al. (2002).

Refinement is on a cell-by-cell basis, controlled by the Truelove et al. (1997) criterion. Note that, with the magnetic Jeans length ${{\lambda }_{B}}={{\lambda }_{J}}\sqrt{1+v_{A}^{2}/c_{s}^{2}}$ (with v_A the Alfvén speed and c_s the sound speed; e.g., Strittmatter 1966), artificial fragmentation is less likely when a magnetic field is included. We adopt the same boundary conditions (pseudo-shearing box) and effective grid resolution of ${{2048}^{2}}\;\times \;4096$ (by using a root grid resolution of ${{128}^{2}}\times 256$ with four refinement levels) as Paper I. Again, we follow evolution for 10 Myr, less than one shear-flow crossing time over the domain.

Inclusion of a uniform magnetic field in the domain introduces a numerical issue: the simulation time step is determined by the low density external medium and is $\lt 1/100$ of that in the disk midplane (region of primary interest). To speed up the simulation, we artificially reduce the magnetosonic sound speed by adopting the Boris method with a maximum magnetosonic sound speed of $2500\;{\rm km}\;{{{\rm s}}^{-1}}$ (Boris 1970; Gombosi et al. 2002). Although this procedure is, in principle, only valid for simulations of steady state problems, it is still applicable when the flow velocities are much smaller than the sound speed.

Following methods of Paper I, to model star formation we allow collisionless star cluster particles, i.e., point masses representing star clusters or sub-clusters of minimum mass ${{M}_{*}}=100\;{{M}_{\odot }}$, to form. These particles are created when the density within a cell exceeds a star formation threshold value of ${{n}_{{\rm H},{\rm sf}}}={{10}^{5}}\;{\rm c}{{{\rm m}}^{-3}}$. We assume a local SFR that converts a fraction, ${{\epsilon }_{{\rm ff}}}=0.02$, of gas above the threshold density into stars per local free-fall time (Krumholz & Tan 2007). No mass-to-flux ratio criterion is included, so magnetically-subcritical cells are also able to form stars, with necessary flux redistribution assumed to be occurring at sub-grid scales. Particle motions are calculated by using the gravitational field interpolated from the grid to the particle positions in the equation of motion. The mass of the particles is included when solving the Poisson equation.

3.1. Cloud Structure and Fragmentation

Figure 1 shows mass surface densities and temperatures within the numerical domain for models with different initial magnetic field strengths. The pure hydrodynamical model has a very fragmented structure, which becomes smoother as magnetic field strength increases. Significant differences are also seen in the temperature distributions, which are the result of different gas densities (and thus different heating/cooling rates and equilibrium temperatures) and also different amounts of shock heating from large scale flows and accretion heating from dense clump formation. Figure 2 shows zoom-ins of two selected regions, to illustrate more clearly the effects of magnetic field strength on resulting gas structures. Differences in star formation activity are discussed below (Section 3.4).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The fragmentation can be quantified by counting the number of separate clouds that form after 10 Myr of evolution. As in Paper I, we assume a threshold density of ${{n}_{{\rm H}}}=100\;{\rm c}{{{\rm m}}^{-3}}$ to define gas that is part of molecular "cloud" structures. For simplicity, all connected cells are counted as a single cloud. Using this routine, we find 287, 134, and 28 clouds in the $0,\;10$, and $80\;\mu$G models, respectively. The mean[median] cloud masses are 4.0 $\times {{10}^{4}}\;[27]\;{{M}_{\odot }}$ for $0\;\mu$G, 7.3 $\times {{10}^{4}}\;[100]\;{{M}_{\odot }}$ for $10\;\mu$G, and 3.1 $\times {{10}^{5}}\;[5.0\times {{10}^{4}}]\;{{M}_{\odot }}$ for $80\;\mu$G. Most clouds forming in the 0 and 10 μG models are thus of very low mass. As expected, increasing the magnetic field strength has a major impact on the fragmentation of the initial GMCs. The magnetic critical mass of spherical clouds, ${{M}_{B}}=16.2{{(B/10\;\mu {\rm G})}^{3}}$ ${{({{n}_{{\rm H}}}/1000\;{\rm c}{{{\rm m}}^{-3}})}^{-2}}$ (Bertoldi & McKee 1992), shows we expect a factor of ∼500 change in typical cloud mass comparing the 10 and 80 μG cases, which is seen in the median cloud masses.

In Figures 1 and 2 we also see that the magnetic field is more distorted by gas motions in Region 1 than Region 2. In the $10\;\mu$G run, filaments predominantly lie parallel or perpendicular to field orientation. On $\sim 100\;{\rm pc}$ scales, the 80 μG fields are hardly affected by the cloud motions, but do become more distorted and disordered in regions of high star formation activity.

3.2. Column and Volume Density Distributions

We also quantify ISM structure by examining the probability distribution functions (PDFs) of ${{{\Sigma }}_{g}}$ (as viewed from above the disk) and number density of H nuclei, ${{n}_{{\rm H}}}$. Figure 3 shows these quantities weighted by area/volume (top row) and by mass (bottom row). The ${{{\Sigma }}_{g}}$ PDF is a directly observable feature of clouds and galaxies, while the volume density PDF needs to be reconstructed using different techniques (e.g., Kainulainen et al. 2014).

Zoom In Zoom Out Reset image size

The ${{{\Sigma }}_{g}}$ PDFs show the typical combination of a log-normal with power-law tail at high surface densities (e.g., Kainulainen et al. 2009). The transition occurs around ${{10}^{-2}}\;{\rm g}\;{\rm c}{{{\rm m}}^{-2}}$, corresponding roughly to the minimum ${{{\Sigma }}_{g}}$ of regions containing "cloud" (GMC) gas, which is essentially all part of self-gravitating structures. There is variation from region to region, with Region 1 having its log-normal peaking at higher ${{{\Sigma }}_{g}}$. The shape of the power-law distribution differs only slightly between these individual regions, roughly scaling as ${\Sigma }_{g}^{-{{\alpha }_{{\Sigma },{\rm PDF}}}}$ with ${{\alpha }_{{\Sigma },{\rm PDF}}}\simeq 1$ for the area-weighted PDF, and thus relatively flat slopes of the mass-weighted distributions. Power law tail normalization at given ${{{\Sigma }}_{g}}$ is several times higher in Region 1 than Region 2, i.e., the former has a higher dense gas mass fraction. This also correlates with its more fragmented morphology and larger SFR (below).

There are only modest differences seen in the PDFs as a function of magnetic field strength (although note the large dynamic range in the figures). At higher ${{{\Sigma }}_{g}}$, these are more clearly seen in Region 2, which also shows the greatest effect from B-field strength on its degree of fragmentation and SFR (below).

Comparing to observations, somewhat steeper indices, i.e., ${{\alpha }_{{\Sigma },{\rm PDF}}}\approx 1.33-3.65$, have been reported in dust-emission-derived ${{{\Sigma }}_{g}}$ PDFs by Schneider et al. (2014). We note that some extinction-based studies find PDFs that can be fit by a single log-normal up to ${\Sigma }\simeq 0.5\;{\rm g}\;{\rm c}{{{\rm m}}^{-2}}$ (Butler et al. 2014a). Note, these observations are derived from "in-plane" higher-resolution views of clouds, in contrast to the "top-down" views of the galactic plane from the simulations. Direct observational constraints of top-down views are expected to be able to achieved with high angular resolution observations of nearby galactic disks, e.g., with ALMA.

Volume density PDFs in the simulations show similar trends as the ${{{\Sigma }}_{g}}$ PDFs, with Region 1 having a larger fraction of gas at higher densities. For ${{n}_{{\rm H}}}\geqslant 100\;{\rm c}{{{\rm m}}^{-3}}$, i.e., the cloud gas, the distributions can be reasonably well approximated with single power laws up to $\sim {{10}^{5}}\;{\rm c}{{{\rm m}}^{-3}}$ with power law index of $\simeq -1$.

3.3. Evolution of Dense Gas Mass Fraction SFRs

Here we consider the time evolution of cloud (${{n}_{{\rm H}}}\gt 100\;{\rm c}{{{\rm m}}^{-3}}$) and clump (${{n}_{{\rm H}}}\gt {{10}^{5}}\;{\rm c}{{{\rm m}}^{-3}}$) mass fractions over the $10\;{\rm Myr}$ timescale of the simulations (Figure 4). For the whole kiloparsec-sized region, cloud mass fractions start from relatively high values, $\simeq 0.6$, inherited from the TT09 simulation. Over $10\;{\rm Myr}$, they show a moderate decline, partly because of the build up a mass fraction in stars, and partly in the magnetized cases from expansion of cloud envelopes due to the introduced magnetic pressure. Clump mass fractions grow to peak values of $\sim 0.05-0.1$, with smaller peak values reached more slowly as the B-field is increased. Since stars are only allowed to form from clump gas, SFR evolution closely tracks clump mass fraction evolution, although it exhibits more stochasticity. At late times the different models appear to converge in their clump mass fractions and SFRs, perhaps due to the exhaustion of most of the initially unstable gas via star formation.

Zoom In Zoom Out Reset image size

Regions 1 and 2 follow similar trends. Indeed Region 1's area of 0.16$\;{\rm kp}{{{\rm c}}^{2}}$ contains large fractions of the clump mass and star formation of the whole kiloparsec region. Region 2 starts with gas structures that are somewhat easier to stabilize with magnetic fields (at least in directions perpendicular to the initial y-direction field orientation), so larger differences are seen in clump and SFR evolution as field strength is increased. Region 1 contains $3\times$ as much cloud gas as Region 2. It contains two GMCs in the process of merging, and the resulting cloud has significant kinetic energy. The magnetic field therefore plays a minor role in this case. The overall SFR surface densities are about $10\times$ larger in Region 1 compared to Region 2.

Some caveats are in order. First, the formation of clumps is a response to simulation initial conditions, where dense, self-gravitating clouds are allowed to collapse to high densities as the resolution is suddenly increased. Second, once formed, star particles contribute gravitationally, but stellar feedback, implementation of which is much more uncertain and numerically challenging, is not yet included in the simulation.

3.4. Magnetic Fields and Average SFRs

The detailed time history of clump and star formation is sensitive to the choice of initial conditions. It is therefore useful to compare the SFR averaged over the full $10\;{\rm Myr}$ evolution to the observations since these are also averaged over $\sim 10\;{\rm Myr}$. Here we consider the absolute and relative average SFRs that are seen in the simulations on the kiloparsec, "region," and "cloud" (i.e., $100\;{\rm pc}$) scales.

Figure 5 shows SFR surface density, ${{{\Sigma }}_{{\rm SFR}}}$, versus ${{{\Sigma }}_{g}}$. Observational results are those of ∼kiloparsec-scale regions of galactic disks from Bigiel et al. (2008) and of individual Galactic GMCs from Heiderman et al. (2010). On the kiloparsec region scale, the simulations all start with ${{{\Sigma }}_{g}}=17.0\;{{M}_{\odot }}\;{\rm p}{{{\rm c}}^{-2}}$, decreasing to $13.9,\;14.8,\;15.2\;{{M}_{\odot }}\;{\rm p}{{{\rm c}}^{-2}}$ after 10 Myr for $0,\;10,\;80\;\mu {\rm G}$ cases, respectively, i.e., ${{{\Sigma }}_{{\rm SFR}}}$ = $0.28,\;0.21,\;0.19\;{{M}_{\odot }}\;{\rm y}{{{\rm r}}^{-1}}\;{\rm kp}{{{\rm c}}^{-2}}$. On these scales, the magnetic field only has a modest impact on SFR, the values of which are much higher ($\sim 10-30\times$) than systems with equivalent ${{{\Sigma }}_{g}}$ observed by Bigiel et al. (2008). However, as noted in Section 3.3, total SFR is dominated by that occurring in Region 1. When separate results for Regions 1 and 2 are considered, we see more significant effects of increasing the magnetic field strength: a factor of 3 reduction in Region 2 (Figure 5).

Zoom In Zoom Out Reset image size

Focusing on smaller $100\;{\rm pc}$-size "Cloud" scales (see Figure 1), we find simulated SFRs overlap with observational rates in Galactic GMCs derived by Heiderman et al. (2010). Stronger magnetic fields again almost always lead to lower SFRs. One reason for the decrease in SFR with increasing B-field is that clumps form later in models with stronger fields (Section 3.3). In one case, Cloud H, stars only form in the $B=0\;\mu {\rm G}$ run.

Our simulations thus reproduce SFRs similar to some observed local Galactic GMCs. Note the Heiderman et al. GMCs do not contain especially vigorous regions of massive star formation, so may be relatively less affected by internal star formation feedback. On the kiloparsec-scale, even our strongly magnetized simulation has values of ${{{\Sigma }}_{{\rm SFR}}}$ that are too large. However, this is largely a consequence of the "starburst" of Region 1. Inclusion of star formation feedback is expected to have an impact in reducing this activity, to be investigated in a future paper.

We have studied the effects of magnetic fields on molecular clouds extracted from a global galaxy simulation, especially for cloud structure and SFRs. Magnetic fields suppress fragmentation, as expected from consideration of the magnetic critical mass. Their effects on ${{{\Sigma }}_{g}}$ and ${{n}_{{\rm H}}}$ PDFs are more modest.

In the context of models of star formation from gas above a threshold density (${{n}_{{\rm H}}}\geqslant {{10}^{5}}\;{\rm c}{{{\rm m}}^{-3}}$), we find that on the largest kiloparsec-scale, average SFRs are only modestly suppressed by magnetic fields. However, this is due to the presence of a starbursting region (Region 1) in the simulation domain. Considering other regions, including down to "GMC"-scales of $\sim 100\;$ pc, we find a variety of ${{{\Sigma }}_{{\rm SFR}}}$s, with values quite similar to those of some Galactic GMCs. Average suppression factors of ${{\epsilon }_{B}}=0.76,\;0.54$ for $B=10,\;80\;\mu$G compared to the nonmagnetized case are seen in the clouds, but in some cases there can be much more dramatic effects, including complete suppression in one cloud.

As a numerical experiment, this study highlights several important aspects, including issues of B-field initialization when overdense structures are already present in initial conditions. Stochastic effects are important, i.e., variation in behavior of individual GMCs, and kiloparsec-scale SFRs can be influenced by a single GMC complex. A sub-grid model for star formation is needed: here via a chosen threshold density and an empirically motivated efficiency per local free-fall time, ${{\epsilon }_{{\rm ff}}}=0.02$. Such low efficiencies may require the effects of local star formation feedback, like protostellar outflows, to maintain turbulence and prolong star cluster formation (Nakamura & Li 2007), and/or they may result from the influence of magnetic fields themselves. Note, we have not included magnetic field effects on the star formation sub-grid model, rather keeping its empirically based parameters fixed for simplicity. There is scope in future work for exploring magnetic sub-grid star formation models, where cell mass-to-flux ratio is also an input, which should lead to a more natural threshold criterion for star formation activity.

Other effects of local star formation feedback, such as stellar winds, ionization, and supernovae, have not yet been included in these simulations, but are expected to act to further reduce SFRs. To disentangle the relative importance of these effects and of magnetic fields will require testing many properties of the simulated clouds and young stellar populations against observed systems (e.g., Butler et al. 2014b).

S.V.L. acknowledges support from the SMA Postdoctoral Fellowship (SAO). J.C.T. acknowledges support from NASA grant ATP 09-0094. Resources supporting this work were provided by NASA High-End Computing Program, the Smithsonian Institution High Performance Cluster and the High Performance Computing facilities at the University of Leeds.