Magnetohydrodynamic Simulations of Disk Galaxy Formation:
the Magnetization of The Cold and Warm Medium
Using magnetohydrodynamic (MHD) adaptive mesh refinement simulations, we study the formation and early evolution of disk galaxies with a magnetized interstellar medium. For a M halo with initial NFW dark matter and gas profiles, we impose a uniform G magnetic field and follow its collapse, disk formation and evolution up to 1 Gyr. Comparing to a purely hydrodynamic simulation with the same initial condition, we find that a protogalactic field of this strength does not significantly influence the global disk properties. At the same time, the initial magnetic fields are quickly amplified by the differentially rotating turbulent disk. After the initial rapid amplification lasting Myr, subsequent field amplification appears self-regulated. As a result, highly magnetized material begin to form above and below the disk. Interestingly, the field strengths in the self-regulated regime agrees well with the observed fields in the Milky Way galaxy both in the warm and the cold HI phase and do not change appreciably with time. Most of the cold phase shows a dispersion of order ten in the magnetic field strength. The global azimuthal magnetic fields reverse at different radii and the amplitude declines as a function of radius of the disk. By comparing the estimated star formation rate (SFR) in hydrodynamic and MHD simulations, we find that after the magnetic field strength saturates, magnetic forces provide further support in the cold gas and lead to a decline of the SFR.
Subject headings:galaxies: formation — galaxies: ISM
Rapid developments in observational cosmology recently established a standard model of cosmology, the CDM model (Spergel et al., 2007). In principle, this makes the problem of galaxy formation essentially an initial value problem. Indeed, Gpc-scale N-body simulations of CDM cosmology are revealing the details of large-scale structure formation and clustering, matching observations in many aspects (Springel et al., 2005). In this framework, it should be possible to understand the formation and evolution processes of individual galaxies. For this purpose, it is crucial to understand the structure of interstellar medium (ISM) and the dominant physics controlling star formation.
Analytical and semi-analytical studies of disk galaxy formation are useful in understanding the global properties (Mo et al., 1998; Efstathiou, 2000; Dutton et al., 2007; Kampakoglou & Silk, 2007; Stringer & Benson, 2007). Because of the complex nature of galaxy formation and the wide range of physical processes involved, numerical simulations also prove particularly useful. Hydrodynamical simulations of galaxy formation in both cosmological and isolated setups have taught us a lot about the detailed structure of the galactic ISM and star formation (see e.g. Kravtsov 2003, Li et al. 2005, Springel & Hernquist 2005, Okamoto et al. 2005, Tasker & Bryan 2006, Tassis et al. 2006, Governato et al. 2007, Wada & Norman 2007, Kaufmann et al. 2007, Robertson & Kravtsov 2007 for recent studies). However, there are still many physical processes we need to include to make the simulations realistic.
In this work, we focus on magnetic fields which have not yet been included in previous galaxy formation simulations. Magnetic fields may significantly influence the structure and evolution of ISM and has been studied extensively previously by local magnetohydrodynamic (MHD) simulations of the ISM (Mac Low et al., 2005; Balsara et al., 2004; de Avillez & Breitschwerdt, 2005; Piontek & Ostriker, 2007; Hennebelle & Inutsuka, 2006) and by observations (e.g. Beck 2007, Crutcher 1999). On a larger scale, the effect of magnetic fields is less clear because cosmological MHD simulations are still at an early stage. Miniati et al. (2001) performed cosmological hydrodynamic simulations for a passive magnetic field. Cosmological MHD simulations have been done in SPH (e.g. Dolag 1999) in the context of galaxy cluster formation and Eulerian-code simulations are also beginning to be explored (e.g. Li et al. 2006). None of those calculations have as yet resolved the internal structures of disk galaxies.
Besides the potential importance in the dynamics of galaxy formation, galactic magnetic fields are a key ingredient in many other astrophysical problems. For example, they play a dominant role in the origin and propagation of cosmic rays (Fermi, 1949; Strong et al., 2007), which may also be dynamically important in galaxy formation (Fermi, 1954; Enßlin et al., 2007). Furthermore, magnetic fields may be important in the dynamics of supernova remnants (Balsara et al., 2001), in regulating the internal turbulence of giant molecular clouds by driving outflows (Li & Nakamura, 2006; Banerjee & Pudritz, 2007) and in the collapse of individual molecular cloud cores (Mouschovias, 1987; Balsara et al., 2001; Tilley & Pudritz, 2007; Mellon & Li, 2007).
While the importance on the ISM is generally accepted, the origin of the galactic magnetic field is still an unsolved problem (see e.g. Kulsrud & Zweibel 2007 for a recent review). The traditional theory of a galactic dynamo, the alpha-Omega dynamo, is based on the mean field dynamo theory introduced by Steenbeck et al. (1966) and later developed by Parker (1971) and Vainshtein & Ruzmaikin (1971). Ferriére (1992) has extended this galactic dynamo theory to include the expansion of supernova bubbles. Turbulent amplification in protogalaxies have also been studied (Pudritz & Silk, 1989; Kulsrud et al., 1997). Those theories are built on kinematic dynamo. A different paradigm for disk dynamo comes with the discovery of the magneto-rotational instability (MRI) (Chandrasekhar, 1961; Balbus & Hawley, 1991). In MRI, field amplification is not only a consequence of the turbulent velocity field, it also produces the turbulent velocity field itself. It is still unclear what is the dominant mechanism that amplifies the galactic magnetic field. In contrast to the case of dynamos in planets, because of the enormous inductance of the ISM (Fermi, 1949), there is no decaying problem for the galactic field. So the essential part of the problem here is indeed only the initial amplification.
In the following, we first describe the numerical algorithms and initial conditions. Then in section 3 we present the results. Sections 4 and 5 are devoted to discussion and conclusions, respectively.
2. Numerical Methods and Models
2.1. MHD Algorithm
We have developed a 3D adaptive mesh refinement (AMR) MHD code based on the cosmological AMR hydrodynamics code Enzo (Bryan & Norman, 1997; O’Shea et al., 2004). We use the Dedner formulation of MHD equations (Dedner et al., 2002). This conservative formulation of the MHD equations uses hyperbolic divergence cleaning. It has also been used in several other recent AMR MHD codes (Matsumoto, 2006; Anderson et al., 2006). The conservative system is discretized by method of lines. The Riemann solver and reconstruction schemes have been tested extensively in Wang et al. (2007) for relativistic systems. Since our Riemann solvers and reconstruction schemes are designed to work for any conservative systems, when applying it to the Dedner MHD equations, we only redefine the conservative variables. Using these high-resolution-shock-capturing (HRSC) schemes designed for conservative systems ensures the conservation of mass, momentum and energy, as well as obtaining correct shock positions and velocities (LeVeque, 2002). Interested readers can find more detailed descriptions of our numerical algorithms in Wang et al. (2007) and references therein. In this work, reconstruction is done using the piecewise linear method (PLM) (Van Leer, 1979). Fluxes at cell interfaces are calculated using the local Lax-Friedrichs Riemann solver (LLF, Kurganov & Tadmor 2000), which is free of the carbuncle artifacts present in some other contact-capturing Riemann solvers (Quirk, 1994). Time integration is performed using the total variation diminishing (TVD) second order Runge-Kutta (RK) scheme (Shu & Osher, 1988). The code inherits all of Enzo’s parallel AMR capability as well as self-gravity, chemistry, cooling and radiative transfer routines. We have tested the code using various standard hydrodynamics and MHD test problems including the 1D Brio-Wu problem, 2D MHD rotor, 2D Hydro and MHD Rayleigh-Taylor problem, 3D Sedov-Taylor problem, 3D Larson-Penston isothermal collapse problem, etc.
In global simulations of galactic disks, current computational power does not allow to resolve the Jeans length down to molecular core scales where star formation actually occurs. This may gives rise to artificial fragmentation if one underresolves the Jeans length (Truelove et al., 1997). This potentially worry is avoided using a density dependent temperature floor ensuring that a cell always resolves half of the local Jeans length (e.g. Wada & Norman 2007):
Another common approach would be to create a collisionless star particle if a cell fails to resolve Jeans length (e.g. Tasker & Bryan 2006).
2.3. Initial Conditions
From pure N-body simulations it has been found that cold dark matter halos have a universal density profile well fitted by the NFW form (Navarro et al., 1996). Here we consider a small halo formed at redshift , with a mass of Mand correspondingly, a virial radius kpc. The spin parameter is assumed to be and concentration is . Such a small halo allows for higher spatial resolution. By using an isolated model, we are implicitly assuming that the effect of major and minor merger on the growth of galaxy disk is subdominant, which may be a reasonable assumption if the halo is smaller than the Milky Way halo (Guo & White, 2007). However, our model galaxy is still larger than the small galaxies with rotation velocities km/s in which star formation may be strongly suppressed by the intergalactic UV background (Thoul & Weinberg, 1996). The gravitational interaction between dark matter and baryons is a secondary effect. Consequently in this work we model the dark matter halo as a static external potential.
The initial gas density profile follows also the NFW profile with its amplitude determined from the assumed baryon mass fraction . The initial gas temperature is set by solving for local hydrostatic equilibrium. The rotation velocity is set to be the same as that of the dark matter following Springel & White (1999). In addition, we also add to the gas random velocities of the same order as the dark matter virial velocity to model crudely the turbulent motions in hierarchical formation of protogalaxies. We put the NFW halo at the center of the simulation box with virial radius (21 kpc) of the simulation box length. This makes sure that the gas evolution inside halos will be minimally affected by numerical boundary effects. The gas outside the halo virial radius is set to be the cosmic mean value at redshift with a temperature of K. The topgrid resolution is , and there are three static nested grids refined by a factor of two around the halo. Then we let the code dynamically refine to higher levels according to the Jeans criterion with Jeans number 4. After the disk forms, almost all the disk is refined to this highest level because it is strongly self-gravitating. The simulations were run with a maximum refinement level of six and seven, corresponding to resolutions of and pc, respectively. We did not find noticeable differences in the results. So in the following, only the results of the pc resolution run are discussed.
The initial magnitude and topology of the magnetic field is uncertain. Cosmological simulations for the generation of magnetic seed field will be required to address this important question. A tiny seed magnetic field of the order G is always guaranteed by the Biermann battery effect (Biermann, 1950). However, there are both observational and theoretical arguments suggesting larger pregalactic field strengths. For example, Faraday rotation has been detected in high redshift damped Lyman alpha system (Oren & Wolfe, 1995). Also, the abundance of beryllium and boron in old Galactic halo stars is directly proportional to their iron abundance. Big bang nucleosynthesis does not produce beryllium nor boron. Hence, they may dominately be produced by spallation of low energy (tens of MeV) carbon and oxygen cosmic rays, suggesting that magnetic fields and cosmic rays are already present at early times (Zweibel, 2003). On a more theoretical ground, in hierarchical structure formation, any halo formed at late times must have had progenitors that hosted prior generations of stars. So any magnetic field produced by those stars, their supernovae, or their pulsar remnants would be amplified by the turbulent ISM in the galaxies containing them. Turbulence driven by mergers may also amplify magnetic fields (Pudritz & Silk, 1989; Kulsrud & Anderson, 1992). Furthermore, Rees (2006) argued that the magnetic fields from supernova ejecta and extended radio lobes may build a field in excess of G seed field for galactic disk forming in the late Universe. In this work, we adopt the simple assumptions that the initial magnetic field is weak, uniform and directed along the rotation axis. We will discuss the results of two simulations, with (Hydro run) and G (MHD run) throughout. Initially, the magnetic field is dynamically unimportant, it will be twisted and amplified by the turbulent velocity field.
This first study of MHD models of disk galaxy formation only include a minimal set of physical processes relevant for the rich interplay between magnetic fields, the dynamics of disk formation and a multiphase ISM. We do not yet include supernova feedback, radiative heating, cosmic rays, molecular cloud chemistry, ambipolar diffusion nor other physical processes that may be relevant. With those limitations in mind, our model should be viewed as a numerical experiment rather than a realistic simulation of current day disk galaxies. Some of those neglected feedback processes, e.g., supernova feedback, may significantly change the conclusions drawing from the current calculations.
Fig. 1 shows the face-on and edge-on projections of the density, temperature and the thermal pressure. Fig. 2 shows face-on and edge on slices of velocity field and magnetic field lines overplotted on density, plasma beta and . Both figures correspond to Gyr. In some of the plots below, we will consider only the gas in the disk defined by a density threshold criterion g cm.
3.1. Disk Structures
As can be seen in Fig. 1, the global disk structures are qualitatively similar in both the Hydro and MHD cases. The density projections show filamentary structures and high density “blobs” located along those filaments. Were we to Include molecular hydrogen chemistry and cooling, these blobs would presumably resemble giant molecular cloud complexes. This is similar to the density distributions of HI gas and GMCs in observations of nearby galaxies such as M33 (Engargiola et al., 2003), which also show good correspondence between the filaments and the locations of the GMCs. The temperature projections show that those filaments are all cold, with K. The lower density disk material has a temperature of K. In the MHD run, Fig. 2 shows that much of the cold gas is already strongly magnetized.
In a self-gravitating disk, the vortex mode is the fastest unstable mode with growth rates times larger than the spiral density wave mode (Mamatsashvili & Chagelishvili, 2007). In an isolated non-rotating fluid, the perturbation mode is divided into vortical and divergent types according to the Helmholz decomposition. And it is mainly the energy in divergent flows that will be dissipated in shocks. In a differentially rotating disk, the vortex mode is a combination of vortical and divergent motions. Thus the morphology and statistical properties of large scale shocks which form the filaments are different in those two cases. However, there seems to be no quantitative comparative studies on this problem that we are aware of.
To see how self-gravity affects the growth of vortex modes, we have carried out experiments with only the external NFW potential, neglecting the self-gravity of the gas. Analytical calculations in the thin sheet approximation predict that the growth rate of vortex mode is orders of magnitude smaller without self-gravity (Mamatsashvili & Chagelishvili, 2007). Indeed in this case we saw that the disk shows much less filamentary structures and instead develops tightly wound spiral density wave patterns. Thus the vortical turbulent motion in the fiducial cases is clearly driven by self-gravity. Another interesting implication of this experiment is that many local galaxies do show tightly wound spiral patterns. Hence our experiment suggests that the gravitational force in the gas disks of those galaxies is dominated by the stars. This implies that for disk galaxies at high redshift that are still self-gravitating, the vortex modes may be the dominant structures regulating molecular cloud and star formation within them. In the following discussion we will call the turbulent velocity field seen in our simulations gravity-driven turbulence (Gammie, 2001; Wada et al., 2002).
As also can be seen in Fig. 1, the MHD run has a somewhat thicker disk than the Hydro run. This shows at Gry, magnetic pressure is already large enough to provide additional support to the vertical force balance, which can also be seen from Fig. 2.
The pressure projections, show that the cold filaments have higher pressure than the warm gas and the cores of the filaments have even higher pressure. The overpressure in the cores is expected as they are self-gravitating (Larson, 1981). The higher pressure in non-selfgravitating cold gas is consistent with previous simulations of isolated disk galaxies (Tasker & Bryan, 2006). This result is different from the idea that cold gas is formed by isobaric thermal instability in the galactic disk (Field, 1965). One important difference with the original thermal instability analysis is that instead of cosmic ray heating, the heating sources in our simulation are work, numerical viscous heating and numerical resistivity in the MHD run, all controlled by the gravity-driven turbulence velocity field in the disk. Thus in our simulations show that gravity-driven turbulent heating and radiative cooling will create a non-isobaric two phase medium. In the old isobaric theory, there is no prediction for the morphology and fraction of the cold phase. But in the turbulent heating scenario, one would expect the cold gas to have filamentary structures, which are naturally produced by a supersonic turbulent velocity field.
Fig. 3 shows the azimuthally averaged radial profiles of mass, surface density, velocity dispersion and rotational velocity for both simulations at Gyr. The Hydro and MHD calculations have similar density profiles and remarkably similar mass profiles. This is because magnetic field amplification happens after the disk is formed so the initial weak magnetic field will not affect significantly the dynamics of disk formation. In both cases, there are no sharp boundaries at the disk edge in surface density profile. From the mass profile one sees the disk radius to be kpc for both simulations. We overplotted the threshold gas surface density for star formation Mpc derived from observations (Martin & Kennicutt, 2001). It is interesting to see that it is just near the disk radius. In the mass profile plot, the dark matter halo mass profile is also shown. It shows that the dark matter mass dominates baryon at kpc. This explains why the rotational velocity profile is roughly flat at kpc with km/s, appropriate for the chosen NFW mass profile. In the rotation velocity plot, there is a big trough at kpc in both the Hydro and MHD curves. This is because of the fact that there are some big clumps at this radius whose spin provides a large cancellation to the mass-averaged rotation velocity. This also gives the trough in Toomre Q parameter discussed below.
The bottom left panel of Fig. 3 shows the velocity dispersions calculated by where represents mass-weighted azimuthal average. The velocity dispersion has been argued to remain roughly constant in self-regulated regions of disks (Silk, 1997), but this is still controversial (Ferguson et al., 1998). We can see that across the disk the velocity dispersions have about an order of magnitude fluctuations centered around km/s. So taking it to be constant is only a rough approximation. However, Silk (1997) argued that the hot phase is the key in creating the self-regulated stage and determining the global velocity dispersions. Since hot phase is completely missing in our simulations, the situation might be different once supernova feedback is included.
Fig. 4 shows the Toomre Q parameter (Toomre, 1964), where is the angular velocity, is the velocity dispersion, is the surface density. In the original Toomre’s formula, is replaced by the epicyclic frequency . We used because our disk is highly turbulent, different averaging methods for calculating would give different . In practice holds to within a factor of two. As can be seen from Fig. 4, throughout the disk Q fluctuates around , implying that the disk is Toomre unstable. Such a roughly constant suggests that the gravity-driven turbulence in the disk is in a quasi-stationary state at this time (Gammie, 2001; Wada et al., 2002). This is also supported by the quasi-stationary density PDF (see section 3.3). Beyond the disk radius, rises sharply to .
3.2. The Growth and Structure of the Magnetic Field
Comparing the density projections of the Hydro and MHD cases in Fig. 1, we can see that the low filaments are more coherent in the MHD case. Correspondingly, the MHD disk has less low density regions. By comparing the velocity field to magnetic field in Fig. 2, we can see that the large-scale velocity field and magnetic field are aligned in many places of the disk. This may be a result of self-organization in MHD turbulence by dynamic alignment of velocity and magnetic fields (Biskamp, 1993). This implies that flows along field lines may be important for the formation of clouds 111We thank Ralph Pudritz for alerting us to this point..
The edge-on slice of magnetic field lines shows that at large scale the poloidal field looks like a split monopole. This is the result of an initial uniform magnetic field along the rotation axis being dragged in by the collapse. The physics of magnetic field amplification in MRI, i.e., the stretching of magnetic field lines increases their energy density at the expense of differential rotation, must also work here. From the face-on slices in Fig. 2, it can be seen that in the horizontal direction the magnetic field amplification is restricted mainly to the disk. But from the edge-on slice, it is clear that the field amplification extends significantly above the disk plane. As a result, there are highly magnetized bubbles with around the galactic disk. Those low beta bubbles form because density decreases sharply above the disk plane but the magnetic field strength decreases much slower. Similar phenomena are well-known in the context of black hole accretion disks (e.g. Stone & Pringle 2001). The existence of such low beta bubbles implies that the magnetic pressure may significantly influence the dynamics above the disk plane. Thus halo fields built in this way may be crucial for halo-disk matter interchange (McKee & Ostriker, 1977; Norman & Ikeuchi, 1989; Gomez de Castro & Pudritz, 1992), in mergers with other galaxies (Toomre & Toomre, 1972; White, 1978), ram-pressure stripping (Gunn & Gott, 1972), galaxy harassment (Moore et al., 1996) and galaxy strangulation (Larson et al., 1980).
Fig. 5 shows the mass-weighted joint temperature-B probability distribution function (PDF). From it we can see that the cold gas has a typical magnetic field strength G. This is remarkably consistent with observations in the Milky Way (Crutcher, 1999). Since , this implies that in cold gas, fluctuates by two order of magnitudes. Thus the average gives only limited information for the magnetization of cold gas. Similar result has also been found in local MHD simulations of molecular cloud (Tilley & Pudritz, 2007). Furthermore, Fig. 5 also shows that magnetic field strength is similar in warm and cold medium (see also Fig. 2). Then since thermal pressure is different in those two phases, this results in the vast difference of plasma beta in those two phases.
Fig. 6 shows the mass-weighted joint density-B PDF. It shows that the density-B relation roughly follows power law relation . For the high density disk gas, . This implies that the Alfven speed is similar in those density ranges. This is also consistent with the observations of relation in dense molecular gas (Crutcher, 1999). The scatter in both observational data and our result are quite large. For the low density halo gas, . This results from the rise-up of the low beta bubble.
To see how the magnetic field components scale with density, Fig. 7 shows the mass-weighted joint density- and density- plots. It can be seen that both components have similar density-dependence. In the high density gas, azimuthal magnetic component dominates over the vertical component, similar to what was found in shearing box simulations of MRI (Hawley et al., 1996). In the low density bubbles, vertical magnetic field dominates over the azimuthal component. Interestingly, this is also consistent with radio polarization observations of nearby edge-on galaxies (Krause et al., 2006).
To examine the amplification history of the magnetic field, the top panel of Fig. 8 shows the time evolution of the total kinetic energy, thermal energy and magnetic energy for the disk gas. The magnetic energy increases exponentially by about five orders of magnitude in the first 600 Myr. Afterwards it remains invariant with a saturated magnetic field strength in the disk. Some of the newly amplified field after saturation will rise above the disk, magnetizing the halo and another fraction of it will be lost by dissipation. Furthermore, after saturation, the magnetic energy is still three orders of magnitude smaller than the total kinetic energy. This explains in part why in our simulation magnetic fields did not influence significantly the global disk properties. It could happen that at higher resolution magnetic fields energy can be amplified to even larger value due to unresolved small-scale dynamo effects and may influence the disk more significantly. Nevertheless, it is encouraging that with the current resolution, the magnetic field strength in the cold gas already is very similar to the observed values in the Milky Way (Crutcher, 1999). This is different from the result of dynamo simulations in MRI-driven turbulence where the total magnetic energy in the saturated state is generally larger than the kinetic energy (Hawley et al., 1996). The primary differences are that turbulence in galactic disk is mainly driven by self-gravity. Furthermore, the simulated galactic disk has a multi-phase medium. The cold phase controls the saturation. And the warm phase still has a subdominant magnetic field in the saturated state. This also contributes to the subdominance of magnetic energy.
The bottom panel of Fig. 8 shows the time evolution of the mass-weighted averaged logarithm of plasma beta for all gas, cold gas and warm gas in the disk. For Myr, the plasma beta in the cold gas quickly deceases from to one and then remains around unity. The magnetic field in the warm phase is always highly supercritical. This evolution behavior suggests that magnetic field amplification in disk galaxies is a self-regulated process. After the plasma beta reaches one in the cold gas, the magnetic fields become dynamically important and buoyant (Parker, 1966). Further growth of magnetic field is quenched in the disk. This could be because of insufficient resolution leading to too large dissipation, yet our pc resolution run gives very similar result. On the other hand, we can see clearly that magnetic fields begin to rise above the disk after saturation, forming the low beta bubbles. This suggests that at least up to Gyr, buoyancy is a relevant mechanism of removing magnetic flux from the disk. To further confirm this picture, we have run simulations that start with a ten times higher initial magnetic field, i.e. G. They reach the self-regulating phase earlier and then drive more low beta bubbles magnetizing the halo in agreement with the preceding interpretation.
The thick line in Fig. 9 shows the mass-weigthed toroidal magnetic field at Gyr. There are regular patterns of field reversal at late times. This is a consequence of field amplification by differentially rotating disks (Parker, 1979). Fig. 9 also shows the time evolution of the mass-weigthed toroidal magnetic field. It can be seen that at early times the averaged toroidal fields are small, before they are amplified by disk rotation. The amplitude generally decreases with radius because the rotation frequency decreases. A perhaps more surprising phenomenon occurred between Myr to Gry, when the toroidal field reversed almost exactly symmetrically.
3.3. Statistical Properties of the ISM
The density PDF is of interest with regards to star formation since star formation is thought to be determined by the high-density tail of the PDF (Elmegreen, 2002; Krumholz & McKee, 2005; Wada & Norman, 2007). A higher order (though not necessarily small!) effect is the spatial distribution of mass. When calculating the density PDF, we considered only gas in the disk. Fig. 10 shows the time evolution of this density PDF for both runs. It can be seen that the PDFs for the two cases are similar. The MHD calculation has slightly less volume in dense gas. In both cases the PDF reaches a quasi-stationary state within less than Myr. This quick establishment of a quasi-stationary PDF suggests that the gravity-driven turbulence in the disk reached a quasi-stationary and self-regulated state. This phenomenon, if robust even including supernova feedback, may play a crucial role in explaining the the observed universal relation between star formation rate and gas density (Kennicutt, 1998; Wong & Blitz, 2002). Our temperature floor is a feedback on the flow that in reality will be dominated by HII regions and supernova feedback.
Our density PDFs cannot be fitted by log-normal distributions. At least, they do not show the typical turn-over of a log-normal distribution. This is different from previous simulations of disk galaxies which claimed that log-normal is a good fit to the density PDF (Tasker & Bryan, 2006, 2007; Wada & Norman, 2007). The crucial difference could be that those studies start from an isolated disk while in our case, the disk is built inside-out dynamically. Thus anytime in the evolution of the disk galaxy ( Gyr) we have simulated, there is low density gas reaching the disk from outside. It is possible that if we continue to evolve the disk long enough that after the initial halo gas all fell down to the disk it may also obtain a log-normal PDF. However, for a realistic cosmological environment, infall should never stop completely.
Fig. 11 shows the joint temperature-density PDF. Above g cm, the temperature is determined by the floor value which was introduced to ensure the calculation resolves the Jeans length. So the temperature of the gas can at most be trusted below this density. From this figure, it is also clear that the ISM develops a two-phase structure with a warm K and a cold K phase. However, unlike the originally quasi-stationary picture of a two-phase ISM (Field, 1965), where the phase is characterized by a density and temperature determined by the balance of heating and cooling, in our case there are three order of magnitude density fluctuations in both cold and warm gas. This is due to the turbulent nature of heating in our simulation as discussed above. However, in both runs, one can still identify a good minimum density g cmfor the cold phase.
Fig. 12 shows the evolution of the volume and mass fractions of cold gas in both runs. It can be seen that the amount of cold gas is monotonically increasing in both volume and mass fractions. After Myr, both runs reach a quasi-stationary volume fraction of and mass fraction .
3.4. Estimated Star Formation History
where is the total gas mass with density which can be computed from the density PDF, is the dynamical time at and is a constant specifying star formation efficiency (SFE).
Using g cmand , Fig. 13 shows the star formation history (SFH) of the two simulations. It can be seen that while the SFR in the Hydro run is constantly increasing, in the MHD case,it reaches a peak at Myr and then decreases. Since Myr is shortly after the time when the average beta in the cold gas reached unity, this suggests that magnetic forces provide further support in the cold gas after the amplification saturated.
Illustrating how sensitive the estimated SFR depends on , Fig. 14 shows as a function of . It can be seen that in both cases it does depend on . Thus to get the same SFR, the SFE should decrease in those densities ranges. However, it is interesting that the effect of magnetic field is to make the dependence weaker. From cm to cm, in the MHD case increases by only a factor of .
The current analysis of the SFR may be a poor predictor of what happens when molecular cloud and star formation and feedback are included more realistically. However, the analysis in this section does suggest that after the magnetic growth saturates in the disk, magnetic fields can significantly influence GMC formation which occurs in the strongly magnetized cold and dense gas.
Using the Spitzer values for viscosity and resistivity, the magnetic Prandtl number of the warm ISM is . Thus, small-scale dynamos are expected to develop. However, those small-scale dynamo effects are not captured by our resolution. Nevertheless, previous simulations (Schekochihin et al., 2004) have shown that for very high magnetic Prandtl number MHD turbulence, the growth of the small-scale magnetic fields is controlled by the rate of strain at the viscous scale and resulted in magnetic energy tending to pile up at the much smaller resistive scale. This is because for small scale dynamos, the magnetic field seems to retain its folded structure in saturation with direction reversals at the resistive scale (Schekochihin et al., 2002). So small-scale dynamo processes may not be important for the large-scale growth of galactic magnetic field. At least, the situation is uncertain. Even if small-scale dynamos would contribute some fraction to the large-scale growth of the galactic magnetic field, our simulations of magnetic field growth would be a lower limit. Since we find that the field can be amplified to the observed value quickly, for the phenomenology of galaxy formation, it seems sufficient not resolve small-scale dynamos at the current level of details.
de Avillez & Breitschwerdt (2005) have discussed the results of 3D AMR MHD simulations of a local patch of ISM with sufficient vertical height to include the effect of halo-gas circulation. They find that gas transport to the halo is not significantly influenced by the magnetic field and they do not find the low beta bubble as seen in our simulations. We think the crucial difference of our simulation from theirs is the presence of global differential rotation, which dominates magnetic field amplification in our simulations.
An important limitation of the presented calculations is that we neglected supernova feedback while their dynamics could be important for amplifying the magnetic field in the ISM (Ferriére, 1992). More importantly, supernova feedback will create the hot phase of the ISM (McKee & Ostriker, 1977), which may play a crucial role in regulating star formation (Silk, 1997). Future work needs to include this supernova feedback to study the importance of superbubbles on the structure of ISM, global star formation and magnetic field amplification. Another limitation is that we neglected cosmic ray pressure, which may be dynamically important. For example, cosmic ray pressure arising from the Alfven wave instability (Kulsrud & Pearce, 1969) may prevent the downsliding of matter along the magnetic field line (Rafikov & Kulsrud, 2000).
The simulations discussed in this work have shown that the amplification of magnetic fields and its large-scale topology is a natural part of the galaxy formation process. Our main findings can be summarized as:
1. In the initial stage of the build-up of a galactic disk, the magnetic field is very weak and does not significantly affect the dynamics of disk formation.
2. In a self-gravitating disk, vortex modes driven by the self-gravity of the gas grow rapidly and trigger the first generation of GMCs formation.
3. In the cold and warm phase created by turbulent heating, the pressure is not isobaric. Temperature and density fluctuates by about three orders of magnitude in both phases.
4. The dynamical formation of a galactic disk results in a global density PDF that cannot be fitted by a log-normal function.
5. Differentially rotating galactic disks amplify a seed field of G to microgauss levels in Myr.
6. The growth of galactic magnetic fields is a self-regulated process. The saturation is reached first in the cold phase. After saturation, the field strength agrees with the observations of Milky Way magnetic fields. The cold phase, the magnetic field strength fluctuates by about an order of magnitude which gives larger than two order of magnitude fluctuations in the associated plasma beta. The saturated total magnetic field energy is three orders of magnitude smaller than the total kinetic energy in the case presented here.
7. The saturation results in highly magnetized material around galactic disks at all but the earliest cosmic epochs. The halo magnetic fields may significantly influence the halo-gas interactions, accretion of gas in disk galaxies, galaxy mergers and the interactions of galaxies with their environment.
8. After saturation, the toroidal field in the disk dominates over vertical components while in the magnetized halo, vertical components dominate over toroidal components.
9. Large scale magnetic field and velocity field are aligned at many places of the disk. This implies that cloud formation is likely channeled by flow along field lines.
10. After saturation, the magnetic field strength is similar in the cold and warm medium.
11. Global toroidal field reversals develops naturally in a differentially rotating disk.
12. Magnetic fields can suppress star formation by providing additional pressure support in the cold gas after saturation.
We have presented the first global numerical experiments studying magnetic field amplification and evolution during the formation of disk galaxies. We find it encouraging that many of the aspects discussed here are consistent with evidence from observations (e.g. compare with Beck 2007). Being able to include more of the relevant physics such as molecular cooling, radiation and supernova feedback a realistic treatment of star formation, and a consistent treatment of cosmic ray acceleration transport, will render the next generation of numerical models useful for a more comprehensive understanding of the physics of galaxies than the current state of the art.
- Anderson et al. (2006) Anderson, M., Hirschmann, E. W., Liebling, S. L., & Neilsen, D. 2006, Classical and Quantum Gravity, 23, 6503
- Balbus & Hawley (1991) Balbus, S. A., & Hawley, J. F., 1991, ApJ, 376, 214
- Balsara et al. (2001) Balsara, D., Benjamin, R. A., & Cox, D. P. 2001, ApJ, 563, 800
- Balsara et al. (2004) Balsara, D. S., Kim, J., Mac Low, M.-M., & Mathews, G. J. 2004, ApJ, 617, 339
- Balsara et al. (2001) Balsara, D., Ward-Thompson, D., & Crutcher, R. M. 2001, MNRAS, 327, 715
- Banerjee & Pudritz (2007) Banerjee, R., & Pudritz, R. E. 2007, ApJ, 660, 479
- Beck (2007) Beck, R. 2007, A&A, 470, 539
- Biermann (1950) Biermann, L., 1950, Z.Naturforsch, 5a, 65
- Biskamp (1993) Biskamp, D., 1993, Nonlinear Magnetohydrodynamics, Cambridge University Press
- Bryan & Norman (1997) Bryan, G. L. & Norman, M. L. 1997a, arXiv:astro-ph/9710187
- Chandrasekhar (1961) Chandrasekhar, S., 1961, Hydrodynamic and Hydromagnetic Stability, Oxford University Press, Oxford
- Crutcher (1999) Crutcher, R. M. 1999, ApJ, 520, 706
- de Avillez & Breitschwerdt (2005) de Avillez, M. A., & Breitschwerdt, D. 2005, A&A, 436, 585
- Dedner et al. (2002) Dedner, A., Kemm, F., Krner, D., Munz, C. -D., Schnitzer, T., & Wesenberg, M., 2002, J. Comput. Phys., 175, 645
- Dolag et al. (1999) Dolag, K., Bartelmann, M., & Lesch, H. 1999, A&A, 348, 351
- Dutton et al. (2007) Dutton, A. A., van den Bosch, F. C., Dekel, A., & Courteau, S. 2007, ApJ, 654, 27
- Efstathiou (2000) Efstathiou, G. 2000, MNRAS, 317, 697
- Elmegreen (2002) Elmegreen, B. G., 2002, ApJ, 577, 206
- Engargiola et al. (2003) Engargiola, G., Plambeck, R. L., Rosolowsky, E., & Blitz, L. 2003, ApJS, 149, 343
- Enßlin et al. (2007) Enßlin, T. A., Pfrommer, C., Springel, V., & Jubelgas, M. 2007, A&A, 473, 41
- Ferguson et al. (1998) Ferguson, A., Wyse, R. F. G., Gallagher, J. S., & Hunter, D. A., 1998, ApJ, 506, L19
- Fermi (1949) Fermi, E., 1949, Phys. Rev. 75, 1169
- Fermi (1954) Fermi, E., 1954, ApJ, 119, 1
- Ferriére (1992) Ferriére, K., 1992, ApJ, 389, 286
- Field (1965) Field, G. B. 1965, ApJ, 142, 531
- Gammie (2001) Gammie, C. F., 2001, ApJ, 553, 174
- Gomez de Castro & Pudritz (1992) Gomez de Castro, A., & Pudritz, R. E. 1992, ApJ, 395, 501
- Governato et al. (2007) Governato, F., Willman, B., Mayer, L., Brooks, A., Stinson, G., Valenzuela, O., Wadsley, J., & Quinn, T. 2007, MNRAS, 374, 1479
- Gunn & Gott (1972) Gunn, J. E., & Gott, J. R. I. 1972, ApJ, 176, 1
- Guo & White (2007) Guo, Q., & White, S. D. M., 2007, astro-ph/0708.1814
- Hawley et al. (1996) Hawley, J. F., Gammie, C. F., & Balbus, S. A., 1996, ApJ, 464, 690
- Hennebelle & Inutsuka (2006) Hennebelle, P., & Inutsuka, S.-i. 2006, ApJ, 647, 404
- Kampakoglou & Silk (2007) Kampakoglou, M., & Silk, J. 2007, MNRAS, 380, 646
- Kaufmann et al. (2007) Kaufmann, T., Wheeler, C., & Bullock, J. S., 2007, astro-ph/0706.0210
- Kennicutt (1998) Kennicutt, R., 1998, ApJ, 498, 541
- Krause et al. (2006) Krause, M., Wielebinski, R., & Dumke, M. 2006, A&A, 448, 133
- Kravtsov (2003) Kravtsov, A. V. 2003, ApJL, 590, L1
- Krumholz & McKee (2005) Krumholz, M. R., & McKee, C. F., 2005, ApJ, 630, 250
- Kulsrud & Anderson (1992) Kulsrud, R. M., & Anderson, S. W., 1992, ApJ, 396, 606
- Kulsrud et al. (1997) Kulsrud, R. M., Cen, R., Ostriker, J. P., & Ryu, D., 1997, ApJ, 480, 481
- Kulsrud & Pearce (1969) Kulsrud, R. M., & Pearce, W. P., 1969, ApJ, 480, 481
- Kulsrud & Zweibel (2007) Kulsrud, R. M., & Zweibel, E. G., 2007, astro-ph/0707.2783
- Kurganov & Tadmor (2000) Kurganov, A. & Tadmor, E. 2000, J. Comput. Phys., 160, 241
- Larson (1981) Larson, R. B. 1981, MNRAS, 194, 809
- Larson et al. (1980) Larson, R. B., Tinsley, B. M., & Caldwell, C. N. 1980, ApJ, 237, 692
- LeVeque (2002) LeVeque, R. J. 2002, Finite Volume Methods for Hyperbolic Problems, Cambridge University Press
- Li et al. (2006) Li, S., Li, H., & Cen, R., 2006, astro-ph/0611863
- Li et al. (2005) Li, Y., Mac Low, M.-M., & Klessen, R. S. 2005, ApJ, 626, 823
- Li & Nakamura (2006) Li, Z.-Y., & Nakamura, F. 2006, ApJL, 640, L187
- Mac Low et al. (2005) Mac Low, M.-M., Balsara, D. S., Kim, J., & de Avillez, M. A. 2005, ApJ, 626, 864
- Mac Low & Klessen (2004) Mac Low, M.-M., & Klessen, R. S. 2004, Reviews of Modern Physics, 76, 125
- Mamatsashvili & Chagelishvili (2007) Mamatsashvili, G. R., & Chagelishvili, G. D., 2007, MNRAS, 381, 809
- Martin & Kennicutt (2001) Martin, C. L., & Kennicutt, R. C., 2001, ApJ, 555, 301
- Matsumoto (2006) Matsumoto, T. 2006, PASJ in press, arXiv:astro-ph/0609105
- Mellon & Li (2007) Mellon, R. R., & Li, Z.-Y., 2007, astro-ph/0709.0445
- McKee & Ostriker (1977) McKee, C. F., & Ostriker, J. P. 1977, ApJ, 218, 148
- Miniati et al. (2001) Miniati, F., Jones, T. W., Kang, H., & Ryu, D. 2001, ApJ, 562, 233
- Mo et al. (1998) Mo, H. J., Mao, S., & White, S. D. M. 1998, MNRAS, 295, 319
- Moore et al. (1996) Moore, B., Katz, N., Lake, G., Dressler, A., & Oemler, A. 1996, Nature, 379, 613
- Mouschovias (1987) Mouschovias, T. C. 1987, NATO ASIC Proc. 210: Physical Processes in Interstellar Clouds, 453
- Navarro et al. (1996) Navarro, J. F., Frenk, C. S., & White, S. D. M., 1996, ApJ, 462, 563
- Norman & Ikeuchi (1989) Norman, C. A., & Ikeuchi, S. 1989, ApJ, 345, 372
- O’Shea et al. (2004) O’Shea, B. W., Bryan, G., Bordner, J., Norman, M. L., Abel, T., Harkness, R., & Kritsuk, A. 2004, In ”Adaptive Mesh Refinement - Theory and Applications”, Eds. T. Plewa, T. Linde & V. G. Weirs, Springer Lecture Notes in Computational Science and Engineering, 2004
- Okamoto et al. (2005) Okamoto, T., Eke, V. R., Frenk, C. S., & Jenkins, A. 2005, MNRAS, 363, 1299
- Oren & Wolfe (1995) Oren, A. L., & Wolfe, A. M., 1995, ApJ, 445, 624
- Parker (1966) Parker, E. N., 1966, ApJ, 145, 811
- Parker (1971) Parker, E. N., 1971, ApJ, 163, 255
- Parker (1979) Parker, E. N., 1979, Cosmical Magnetic Fields: Their Origin and Their Activity, Oxford University Press
- Piontek & Ostriker (2007) Piontek, R. A., & Ostriker, E. C. 2007, ApJ, 663, 183
- Pudritz & Silk (1989) Pudritz, R. E., & Silk, J., 1989, ApJ, 342, 650
- Quirk (1994) Quirk, J. 1994, Int. J. Numer. Methods Fluids, 18, 555
- Rafikov & Kulsrud (2000) Rafikov, R. R., & Kulsrud, R. M., 2000, MNRAS, 314, 839
- Rees (2006) Rees, M. J., 2006, Astron. Nachr., 327, 395
- Robertson & Kravtsov (2007) Robertson, B., & Kravtsov, A., 2007, astro-ph/0710.2102
- Rosen & Bregman (1995) Rosen, A., & Bregman, J. N., 1995, ApJ, 440, 634
- Sarazin & White (1987) Sarazin, C. L., & White, R. E., 1987, ApJ, 320, 32
- Schekochihin et al. (2002) Schekochihin, A. A., Cowley, C. S., Hammett, G. W., Maron, J. L., & McWilliams, J.C., 2002, New J. Phys., 4, 84
- Schekochihin et al. (2004) Schekochihin, A. A., Cowley, C. S., Taylor, S. F., Maron, J. L., McWilliams, J. C., 2004, ApJ, 612, 276
- Shu & Osher (1988) Shu, C.-W. & Osher, S., 1988, J. Comput. Phys., 77, 439
- Silk (1997) Silk, J., 1997, ApJ, 458, 703
- Stringer & Benson (2007) Stringer, M. J., & Benson, A. J. 2007, ArXiv Astrophysics e-prints, arXiv:astro-ph/0703380
- Spergel et al. (2007) Spergel, D. N., et al., 2007, ApJS, 170, 377
- Springel et al. (2005) Springel, V., et al. 2005, Nature 435, 629
- Springel & Hernquist (2005) Springel, V., & Hernquist, L. 2005, ApJL, 622, L9
- Springel & White (1999) Springel, V., & White, S. D. M., 1999, MNRAS, 307, 162
- Steenbeck et al. (1966) Steenbeck, M., Krause, F., & Rädler, Ke-H, 1966, Z. Naturforsch, 21a, 369
- Stone & Pringle (2001) Stone, J. M., & Pringle, J. E., 2001, MNRAS, 322, 461
- Strong et al. (2007) Strong, A. W., Moskalenko, I. V., & Ptuskin, V. S., 2007, Annu. Rev. Nucl. Part. Sci., 57, 285
- Tasker & Bryan (2006) Tasker, E. J., & Bryan, G. L., 2006, ApJ, 641, 878
- Tasker & Bryan (2007) Tasker, E. J., & Bryan, G. L., 2007, astro-ph/0709.1972
- Tassis et al. (2006) Tassis, K., Kravtsov, A. V., & Gnedin, N. Y. 2006, astro-ph/0609763
- Thoul & Weinberg (1996) Thoul, A. A., & Weinberg, D. H. 1996, ApJ, 465, 608
- Tilley & Pudritz (2007) Tilley, D. A., & Pudritz, R. E. 2007, MNRAS, 930
- Truelove et al. (1997) Truelove, J. K., Klein, R. I., McKee, C. F., Holliman, J. H., II, Howell, L. H., & Greenough, J. A. 1997, ApJL, 489, L179
- Toomre (1964) Toomre, A., 1964, ApJ, 139, 1217
- Toomre & Toomre (1972) Toomre, A., & Toomre, J. 1972, ApJ, 178, 623
- Van Leer (1979) Van Leer, B. 1979, J. Comput. Phys., 32, 101
- Vainshtein & Ruzmaikin (1971) Vainshtein, S. I., & Ruzmaikin, A. A., 1971, Sov. Astron. 15, 74
- Wada et al. (2002) Wada, K., Meurer, G., & Norman, C. A., 2002, ApJ, 577, 197
- Wada & Norman (2007) Wada, K., & Norman, C. A., 2007, ApJ, 660, 276
- Wang et al. (2007) Wang, P., Abel, T., & Zhang, W., 2007, astro-ph/0703742
- White (1978) White, S. D. M. 1978, MNRAS, 184, 185
- Wong & Blitz (2002) Wong, T., & Blitz, L., 2002, ApJ, 569, 157
- Zweibel (2003) Zweibel, E. G., 2003, ApJ, 587, 625