The Astrophysical Journal, 580: , 2002 November 20 # The American Astronomical Society. All rights reserved. Printed in U.S.A. THREE-DIMENSIONAL MAGNETOHYDRODYNAMIC MODELING OF THE GASEOUS

Please download to get full document.

View again

of 18
All materials on our website are shared by users. If you have any questions about copyright issues, please report us to resolve them. We are always happy to assist you.


Publish on:

Views: 45 | Pages: 18

Extension: PDF | Download: 0

The Astrophysical Journal, 580: , 2002 November 20 # The American Astronomical Society. All rights reserved. Printed in U.S.A. THREE-DIMENSIONAL MAGNETOHYDRODYNAMIC MODELING OF THE GASEOUS STRUCTURE OF THE GALAXY: SETUP AND INITIAL RESULTS Gilberto C. Gómez Department of Astronomy, University of Wisconsin Madison, 475 North Charter Street, Madison, WI 53706; and Donald P. Cox Department of Physics, University of Wisconsin Madison, 1150 University Avenue, Madison, WI; Received 2002 March 29; accepted 2002 July 25 ABSTRACT We show the initial results of our three-dimensional MHD simulations of the flow of the Galactic atmosphere as it responds to a spiral perturbation in the potential. In our standard case, as the gas approaches the arm, there is a downward converging flow that terminates in a complex of shocks just ahead of the midplane density peak. The density maximum slants forward at high z, preceded by a similarly leaning shock. The latter diverts the flow upward and over the arm, as in a hydraulic jump. Behind the gaseous arm, the flow falls again, generating further secondary shocks as it approaches the lower z material. In cases with two arms in the perturbing potential, the gaseous arms tends to lie somewhat downstream of the potential minimum. In the four-arm case, this is true at large r or early evolution times. At smaller r, the gaseous arms follow a tighter spiral, crossing the potential maximum and fragmenting into sections arranged on average to follow the potential spiral. Structures similar to the high-z part of the gaseous arms are found in the interarm region of our two-arm case, while broken arms and low column density bridges are present in the four-arm case. Greater structure is expected when we include cooling of denser regions. We present three examples of what can be learned from these models. We compare the velocity field with that of purely circular rotation and find that an observer inside the Galaxy should see radial velocity deviations typically greater than 20 km s 1.Synthetic spectra, vertical from the midplane, show features at velocities 20 km s 1, which do not correspond to actual density concentrations. Placing the simulated observer outside the Galaxy, we find velocity structure and arm corrugation similar to those observed in H in NGC Subject headings: galaxies: spiral galaxies: structure ISM: kinematics and dynamics MHD 1. INTRODUCTION 235 Even though spiral structure is one of the most prominent features of disk galaxies, details of the spiral arms in our own Galaxy remain uncertain. Georgelin & Georgelin (1976) traced the spiral structure of the Milky Way using H ii regions and developed a model with four arms. More recent attempts concluded that the Milky Way might actually have a superposition of two- and four-arm structures, each one with different pitch angles, which might arise from different components of the Galactic disk (Drimmel 2000; Lépine, Mishurov, & Dedikov 2001), suggesting that the stellar and gaseous disks might not be tightly coupled. Similar behavior has been frequently observed in external galaxies (e.g., Puerari & Dottori 1992; Grosbøl & Patsis 1998). Roberts (1969) showed that the gas must generate a largescale shock in the presence of a spiral perturbation. It was proposed that the density enhancement induced by this shock might generate a sequence of molecular clouds and star formation downstream from the shock, which itself was associated with the strong dust lane observed in the inner region of the spiral arms in external galaxies. Two-dimensional numerical models by Tubbs (1980) and Soukup & Yuan (1981) showed that the gas forms a vertical shock perpendicular to the plane of the Galactic disk. The postshock gas remained close to hydrostatic equilibrium, even with an adiabatic equation of state. Their results did not show vertical motions larger than 3 km s 1. In fact, the largest downflow that they found was due to the preshock gas readjusting its vertical structure as it flows into the arm potential. Therefore, when H i observations of face-on galaxies showed extended velocity components with dispersions of the order of 20 km s 1, they were attributed to other phenomena, such as galactic fountains, a warping of the H i disk, or intermediate-velocity clouds (Dickey, Hanson, & Helou 1990; Kamphuis & Briggs 1992; Kamphuis & Sancisi 1993). Since then, we have realized that the interstellar medium (ISM) is thicker and has a higher pressure than previously thought. The pressure scale height has been found to be larger than the density scale height, and the nonthermal pressures (turbulent, magnetic, and cosmic ray) are at least as large as the thermal component (Badhwar & Stephens 1977; Reynolds 1989; Boulares & Cox 1990). Therefore, less compressible gas needs to be considered in order to generate more realistic models of the ISM. Such a medium, with a larger effective (the ratio of the specific heats) would be more likely to display the vertical motions characteristic of a hydraulic jump. With this in mind, Martos & Cox (1998, hereafter MC) performed two-dimensional MHD simulations of the flow of the gaseous disk and found diverse structures that differed from the vertical near-hydrostatics found in previous studies. In many cases, the gas moved up ahead of the stellar arm, sped up over it, and fell behind with large bulk velocity. Frequently, there was a downstream shock at higher z as this downflow was arrested, sometimes resulting in secondary midplane density maxima. 236 GO MEZ & COX Vol. 580 The goal of our investigation is to extend calculations like those of MC to three dimensions, to a large fraction of the Galaxy, and to look for the possible observational signatures of these models. In this paper, we present the early results of these simulations. In x 2 we describe the numerical setup and the procedure to achieve the initial hydrostatic equilibrium, in x 3 we describe the results of the simulations, in x 4 we present three examples of synthetic observations that can be done with this type of simulation, and in x 5we present our conclusions. 2. THE NUMERICAL SETUP We performed three-dimensional MHD simulations in polar coordinates using the code ZEUS (Stone & Norman 1992a, 1992b; Stone, Mihalas, & Norman 1992). This code solves the ideal MHD equations for an inviscid fluid with infinite conductivity in a fixed Eulerian grid. For our standard case, the grid extends from 0 to 1 kpc in z, 3 to 11 kpc in r, and 0 to 2=N in the azimuthal angle with 50, 80, and 200=N grid points in each respective direction, N being the number of arms in each case. The boundaries in r and z are reflective, while those in are periodic Hydrostatics, Theory The azimuthally averaged gravitational potential used is model 2 from Dehnen & Binney (1998). Using this potential, we set up a hydrostatic ISM based on the scheme introduced by York et al. (1982) and described in the Appendix of R. A. Benjamin (2002, in preparation). We extended the procedure to include the effects of an azimuthal magnetic field whose magnitude depends only on the density. Given a den- Fig. 1. Initial state of the (1, 0.175, 2) case. Cases are labeled with (T=10 4 K, p mag =10 12 dyn cm 2, N), where T is the temperature, p mag is the coefficient in eq. (8), and N is the number of spiral arms. The gas is in hydrostatic and dynamical equilibrium in the vertical and radial directions with a temperature of 10 4 K. The density is presented assuming a mean particle mass l ¼ 1:27m H. No. 1, 2002 MODELING GASEOUS STRUCTURE OF GALAXY 237 sity profile in the midplane, an equation of state (isothermal, in our case), and a density magnetic pressure relation, the density and rotation velocity are uniquely defined everywhere. Define the function GðÞ as Z GðÞ ¼ d d ð p T þ p B Þ d ; ð1þ where is the density and p T and p B are the thermal and magnetic pressures, respectively. The vertical p T þ p B ¼ where is the gravitational potential, reduce via equation (1) þ Þ ¼ ð2þ ) G½ðr; zþš ¼ G½ðr; 0ÞŠ þ ðr; 0Þ ðr; zþ : ð3þ At any z, the velocity profile is given by the radial balance between the radial potential and pressure gradients, the magnetic tension, and the centrifugal force: v 2 ðr; zþ þ ðr; ð p T þ p B Þþ 2p Bðr; zþ : ð4þ r This balance can be reduced þ ¼ v2 2p B= r ¼ v2 v2 A r ; ð5þ where v A is the Alfvén velocity. Provided 2 ðg þ 2 ðg þ along with equations (2) and (5), we 2 v2 A Þ ¼ ð6þ ) v 2 ðr; zþ ¼v 2 ðr; 0Þ v 2 Aðr; 0Þþv 2 Aðr; zþ : ð7þ Fig. 2. Upper panel shows the column density of the simulation for the (1, 0.175, 2) case at 248 Myr. Lower panel shows density along a cylindrical surface at r ¼ 8 kpc. White contours show density increasing in factors of 10, starting at 10 3 cm 3. 238 GO MEZ & COX Vol. 580 Given the midplane density distribution, ðr; 0Þ, the density off the plane is obtained by solving equation (3) for ðr; zþ, while v ðr; 0Þ is obtained from equation (4) evaluated at z ¼ 0, after which the equilibrium rotation speed at all other z follows from equation (7) Hydrostatics, Implementation MC performed their two-dimensional MHD calculations of spiral arm structure using the vertical density distribution at the solar circle compiled by Boulares & Cox (1990), modified to have a slightly lower vertical scale height for the warm ionized component. We found that the vertical distribution can be reproduced fairly accurately with thermal and magnetic pressures as follows. The thermal component assumes a neutral gas with a constant temperature of 5700 K and an isothermal equation of state. The magnetic pressure is taken as p B ¼ 1: n dyn cm 2 ; ð8þ n þ n c where n c ¼ 0:04 cm 3. The form of the magnetic pressure is such that it has little gradient at high density and is proportional to at low density. The former accommodates a dense thermally supported core near the midplane, while the latter leads to a higher but constant signal speed at low density, far off the plane. With a helium abundance equal to 10% of the hydrogen abundance by number, the mean atomic mass is ð14=11þm H. In our initial work, we wanted to explore a situation that was not so heavily dominated by magnetic pressure. We therefore raised the temperature to 10 4 K and reduced the magnetic pressure by a factor of 10, keeping n c ¼ 0:04 cm 3. The midplane density distribution was taken as exponential, Fig. 3. Same as Fig. 2, with velocity structures superimposed. The solid lines in the upper panel show the clockwise velocity field of the gas in the midplane. In the lower panel, the direction and relative length of the arrows is correct, but the length of the individual components is not. See discussion in the text. No. 1, 2002 MODELING GASEOUS STRUCTURE OF GALAXY 239 with a radial scale length of 4 kpc and a density of 1.11 cm 3 at r ¼ 8 kpc. When these parameters are introduced into the above formalism and the hydrostatics found, the vertical half-disk column density at r ¼ 8 kpc is 0.12 kpc cm 3.Figure1 shows the density, rotation velocity, and magnetic field strength versus radius at the midplane and versus z at r ¼ 8 kpc. The midplane density varies by less than a factor of 10, while vertically the density drops nearly 4 orders of magnitude between the midplane and z ¼ 1 kpc, our present maximum height. At smaller radii, the vertical gravity is stronger and the density gradient in z even larger, so that both the highest and lowest densities occur at the inner boundary. A curious feature is that at high z the density increases with increasing radius: the disk flares. The rotation velocity varies only slightly with radius, by less than 15 km s 1, and by much less with z, only about 1 km s 1 increase (as per eq. [7], higher Alfvén speed requires higher rotation rate at high z in this approximation). The magnetic field strength varies only slowly with radius and appears roughly Gaussian in z, with a flat region in the inner 200 pc of the thermal core. We also report below on a case that has no magnetic field, in which the constant value of the temperature was taken as 2: K. More precisely, the thermal pressure at a given density was taken as 2.5 times that of the previous run, because the above temperature was used inconsistently with assuming that the gas was still neutral. This increase in thermal pressure was made in order to have a density distribution roughly similar to the magnetic case The Spiral Perturbation In addition to the axisymmetric potential, we used a spiral perturbation of fixed shape that rotates with P ¼ 12 km s 1 kpc 1 and has a pitch angle of 15. Details are reported in Cox & Gómez (2002). All the simulation grid is inside corotation. The depth of the perturbation varies slightly in r, weakens in z, and has a sinusoidal profile in ; in the midplane, its corresponding mass density amplitude is 52% of the disk component of the axisymmetric model at r ¼ 8 kpc, Fig. 4. Midplane perturbation potential and density at different heights for the (1, 0.175, 2) case at 248 Myr, along ¼ =2 and r ¼ 8 kpc. Dotted lines show the initial (hydrostatic) density. 240 GO MEZ & COX Vol. 580 which provides a peak-to-valley potential difference of about (30 km s 1 ) 2. This mass contrast is consistent with K- band observations performed by Rix & Rieke (1993), Rix & Zaritsky (1995), and Kranz, Slyz, & Rix (2001), who quote an arm/interarm contrast in density of old stars between 1.8 and 3 for a sample of spiral galaxies Numerical Complications In early runs, performed with outflow boundary conditions at the inner, outer, and upper boundaries, ZEUS soon reported difficulties with hot zones in which the time step became so short that the calculation terminated. This is almost certainly caused by the enormous density contrast of our hydrostatic solution. By making those three boundaries reflecting, so that the material is unable to flow off the grid, this problem was postponed or eliminated, depending on the case run. The two-arm magnetized case ran to about 270 Myr, the unmagnetized case and the four-arm magnetized case ran the full 400 Myr asked of them. In addition to changing the boundary conditions, we also changed the perturbation force field near the inner boundary after noticing that the spiral potential was pushing material against the inner boundary, causing reflected waves that propagated outward. In order to avoid splashing against the inner r boundary, this perturbation is not applied in the inner 1 kpc of the grid, while it is smoothly turned on in the subsequent kiloparsec. Thus, the useful computational grid runs from 5 to 11 kpc. Also, in order to diminish initial transient effects, the spiral perturbation is turned on gradually during the first 50 Myr. This short turnon time undoubtedly creates part of the transient behavior and may have exaggerated some of the early velocity structure. Our intention is to make runs lasting so long that such Fig. 5. Divergence of the velocity for (1, 0.175, 2) case at 248 Myr. Only those places with x v 0 are shown. The upper panel shows the midplane distribution, peaking at the inner edges of the column density arms of Fig. 2. There is a substantial forward-leaning shock (between = ¼ 0:6 and 0.7) preceding the forward-leaning density ridge of Fig. 2 (the high density at the base of the latter dominates the column density maps). But there are a variety of other shocks as well. D No. 1, 2002 MODELING GASEOUS STRUCTURE OF GALAXY 241 transients have died out and to report on that asymptotic behavior in future work. We are wary of artifacts that might be caused by our boundary conditions and will continue to experiment with alternatives, including cases with lower overall density contrast that might allow open boundaries. 3. BEHAVIOR OF THE SIMULATIONS Our primary example is the two-arm spiral with moderate (reduced from MC) magnetic pressure. Our calculation space was the positive z of half a circular disk, with periodic boundaries in. The period of time for a mass element to rotate around this half-disk, relative to the rotating pattern, is about 100, 200, and 340 Myr at r ¼ 5, 8, and 10 kpc, respectively. We have chosen a fiducial time of 248 Myr, for our initial examination of the structure. At this time, nearly all mass elements have experienced the spiral perturbation once or twice, but conditions are still transient, representative of local interaction with rather than global accommodation to the perturbation. We will compare this early structure with that of the unmagnetized case at the same time, and the four-arm magnetized case at half that time, which roughly represents the same level of maturity. Having examined those single early time characteristics, we will present features of the subsequent development of the two-arm cases and the much more mature four-arm case, as indicative of features requiring a longer time to appear. In the remainder of this paper, we will refer to each case using a three-element naming scheme, describing the tem- Fig. 6. Vertical velocity structure at z ¼ :31 kpc for the (1, 0.175, 2) case at 248 Myr. White contours in the upper panel are drawn at v z ¼ 20, 0, and 20 km s 1. Lower panel shows a cut along the vertical ( ¼ =2) line. Comparison of the upper panel with that of Fig. 2 shows that the velocity field at this height has twice the frequency of the column density; i.e., it looks like a four-arm spiral. The sawtooth pattern of the lower panel is consistent with the velocity field of Fig. 3 and x v in Fig. 5. There is a forward-leaning shock (in z; see Fig. 5) preceding the forward-leaning density ridge of the gaseous arm. Upstream from the shock, the gas is falling; downstream, it is rising as in a hydraulic jump. Immediately over the arm, the vertical velocity is close to zero; gas flows up and over the density ridge and down the other side. As found by MC, the falling gas on the downstream side again shocks. In their case, this led to secondary density maxima in the midplane. In our runs thus far (e.g., Fig. 4), there are secondary density maxima at high z, but they are too weak to show up in column density maps. D 242 GO MEZ & COX Vol. 580 perature (in units of 10 4 K), the numerical coefficient in the magnetic pressure density relation (eq. [8]; in units of dyn cm 2 ), and the number of spiral arms in the perturbation potential. Therefore, the magnetic two-arm case will be denoted (1, 0.175, 2), the nonmagnetic two-arm case (2.5, 0, 2), and the four-arm magnetic case (1, 0.175, 4) Two-Arm Magnetized Case Figures 2 and 3 show the results for the (1, 0.175, 2) case. The lower panels show the density and (in Fig. 3) velocity field along a surface of constant radius r ¼ 8 kpc. The upper panels show the column density for z 0, half the total. In the upper panels, rotation is clockwise, and in Figure 3, the lines show the integrated velocity field at the midplane. For clarity in the visualization, we modified the components of the velocity in the lower panel of this figure, so the arrows representing velocity in the inertial reference frame are parallel to the flowlines, with relative lengths proportional to the total velocity. The most important feature in Figures 2 and 3 is the presence of a simple grand design spiral. The density concentration in the midplane contains most of the column density. At higher z, this feature leans forward. The midplane gaseous arm appears slightly after the perturbation potential minimum (outward in radius), shifting to a better alignment farther above the plane (Fig. 4). Note also the strong arm-to-interarm contrast; a significant fraction of the material is located in the arms. The convergences of the velocity flow field in the upper panel of Figure 3 are consistent with this concentration. Figure 5 shows the negative of the velocity divergence at the same positions as Figure 3. Only the regions with neg
Related Search
Similar documents
View more...
We Need Your Support
Thank you for visiting our website and your interest in our free products and services. We are nonprofit website to share and download documents. To the running of this website, we need your help to support us.

Thanks to everyone for your continued support.

No, Thanks