Noncircular Motion and Outflows in the Galactic Bulge

Stochastic Noncircular Motion and Outflows Driven by Magnetic Activity in the Galactic Bulge Region


By performing a global magneto-hydrodynamical simulation for the Milky Way with an axisymmetric gravitational potential, we propose that spatially dependent amplification of magnetic fields possibly explains the observed noncircular motion of the gas in the Galactic center region. The radial distribution of the rotation frequency in the bulge region is not monotonic in general. The amplification of the magnetic field is enhanced in regions with stronger differential rotation, because magnetorotational instability and field-line stretching are more effective. The strength of the amplified magnetic field reaches mG, and radial flows of the gas are excited by the inhomogeneous transport of angular momentum through turbulent magnetic field that is amplified in a spatially dependent manner. In addition, the magnetic pressure-gradient force also drives radial flows in a similar manner. As a result, the simulated position-velocity diagram exhibits a time-dependent asymmetric parallelogram-shape owing to the intermittency of the magnetic turbulence; the present model provides a viable alternative to the bar-potential-driven model for the parallelogram-shape of the central molecular zone. This is a natural extension into the central few 100 pc of the magnetic activity, which is observed as molecular loops at radii from a few 100 pc to 1 kpc. Furthermore, the time-averaged net gas flow is directed outward, whereas the flows are highly time-dependent, which we discuss from a viewpoint of the outflow from the bulge.

accretion, accretion disks — Galaxy: bulge — Galaxy: centre — Galaxy: kinematics and dynamics — magnetohydrodynamics (MHD) — turbulence

1 Introduction

Observations of the atomic and molecular interstellar medium in the inner Milky Way exhibit large noncircular motion of the gas (Rougoor & Oort, 1960; Scoville, 1972; Liszt & Burton, 1978; Bally et al., 1987; Oka et al., 1998; Tsuboi, Handa & Ukita, 1999; Sawada et al., 2004; Oka et al., 2005; Takeuchi et al., 2010). Among various possible explanations for the noncircular motion, de Vaucouleurs (1964) raised a possibility that the gas responds to a stellar bar potential of the Milky Way. Peters (1975) adopted a model that takes into accout gas flowing along elliptic streamlines caused by a bar-like potential and showed that the general trend of the neutral hydrogen (HI) distribution is explained by the model. Liszt & Burton (1980) developed a tilted bar model that accounts for the observed HI distributions at different latitudes.

Observationally, Blitz & Spergel (1991) identified a bar-like distribution of stars in the near-infrared distribution by Matsumoto et al. (1982, see also ())3. Based on the above consideration and finding, Binney et al. (1991) showed that a bar-like potential naturally reproduces the noncircular motion by comparing observed and calculated (Galactic longitude – line-of-sight velocity from the Local Standard of Rest; LSR, hereafter) diagrams. The streaming motion driven by the bar potential provides an only viable explanation on the parallelogram, which has been explicitly discussed in literatures (Blitz, 1994; Koda & Wada, 2002; Rodriguez-Fernandez & Combes, 2008; Baba, Saitoh & Wada, 2010; Molinari et al., 2011).

We however find that the observed parallelogram in the 2.6mm CO emission shows asymmetry in the positive and negative velocities and in , including CO features only at positive longitudes, 3 degrees (clump 2) and 5.5 degrees (e.g., Bania, 1977; Torii et al., 2010b). We also see that the central molecular zone (CMZ, hereafter) shows vertical CO features up to vertical distance pc (Enokiya et al., 2014; Torii et al., 2014a, b), which we discuss further later in this section. We note that these outstanding observed characteristics remain unexplored by the bar potential model.

We note that few efforts have actually been made to elaborate detailed observed features of the parallelogram until recently by numerical simulations in the bar potential. Rodriguez-Fernandez & Combes (2008) presented the first numerical simulations for the parallelogram. Their Figure 14 shows a simulated parallelogram to be compared to the CMZ, while their result is yet crude at best with no detailed fitting to the gas distribution in the - diagram. These authors concluded that the observed asymmetry of the CMZ cannot be due to lopsidedness of the bar and suggested that the asymmetry may be of an ad-hoc origin like infalling gas into the CMZ at =1.3. Their conclusion indicates that the bar potential alone cannot reproduce the asymmetry of the CMZ or the high-z CO features including the 1.3 complex. The difficulty in explaining the CMZ makes a contrast with the large-scale simulations in the bar potential over several kpc, which successfully reproduce the 3-kpc arm and other major features outside the central kpc (e.g., Rodriguez-Fernandez & Combes, 2008).

On the other hand, magnetic field is also supposed to play an important role in the Galactic bulge. Complicated structures, e.g., nonthermal filaments including Radio Arc, are observed in the Galactic center region, which probably attribute to magnetic fields (Yusef-Zadeh, Morris & Chance, 1984; Tsuboi et al., 1986; Chuss et al., 2003; Nishiyama et al., 2010; Morris, 2014). Crocker et al. (2010) gave a lower limit on the field strength, G, over the central 400 pc region from a non-thermal radio spectrum. By other considerations, it is argued that inferred magnetic field strength there is as strong as mG (Morris et al., 1992; Ferrière, 2009).

If such strong magnetic field is distributed in the bulge, it affects the global gas dynamics there (e.g., Sofue, 2007). In fact, Fukui et al. (2006) discovered two large molecular loops, loop 1 and loop 2, at a radius of 700-pc from the center, which are triggered by magnetic buoyancy (Parker instability; Parker, 1966); while the main focus in Parker (1966) was aimed at the Galactic disk, he already suggested that the magnetic activity might become even more important in the central few 100 pc. Subsequently, Fujishita et al. (2009); Torii et al. (2010a, b); Riquelme et al. (2010); Kudo et al. (2011) presented analyses of more detailed observational properties of these molecular loops, and it was shown that the molecular loops are also well reproduced by magneto-hydrodynamical (MHD, hereafter) simulations (Shibata & Matsumoto, 1991; Machida et al., 2009; Takahashi et al., 2009).

In spite of the broad recognition of the strong magnetic field in the Galactic center, there is little theoretical work that explores the role of the magnetic field in the dynamics of the gas component, except a limited number of attempts (e.g., Machida et al., 2009; Kim & Stone, 2012; Machida et al., 2013). We naturally expect that the magnetic field may play an essential role in the noncircular motion of the gas in the Galactic center region, which is the main aim of the present paper.

In this paper, we investigate the evolution and the role of magnetic field in the central region of the Milky Way by a three-dimensional (3D) MHD simulation. We do not take into account a bar potential but an axisymmetric potential to focus on the role of the magnetic field in exciting the noncircular motion of gas.

2 Simulation Setup

Figure 1: Radial profile of the equilibrium initial condition at the Galactic plane. top: Sound speed (dashed) and initial rotation speed (solid) at the Galactic plane. The rotation speed of the stellar component, (dotted) is also shown for comparison. Bottom: Initial density (dashed; left axis), gas pressure (solid; right axis), and magnetic pressure, (dotted; right axis). is multiplied by 1000 times to fit in the panel.

We treat the evolution of gas by solving MHD equations under an external axisymmetric Galactic gravitational potential. We consider three components of gravitational sources: The supermassive blackhole (SMBH) at the Galactic center (component ), a stellar bulge (), and a stellar disk (). For the SMBH we assume a point mass, (Genzel, Eisenhauer & Gillessen, 2010), where is the solar mass. For the and 3 components, we adopt a gravitational potential introduced by Miyamoto & Nagai (1975) for the bulge and disk components. The gravitation potential that includes these three components is written as


where and are cylindrical radius and vertical distance, respectively, and the adopted values for , , and are summarized in Table 1; since we assume a point mass for the SMBH, and for the bulge and disk components we use the same values for , , and in Miyamoto & Nagai (1975) (see also Machida et al. (2009)).

(kpc) (kpc)
SMBH 1 0 0
Bulge 2 2.05 0 0.495
Disk 3 25.47 7.258 0.52
Table 1: Parameters for the gravitational potential. The parameters for the super-massive black hole (SMBH) are from Genzel, Eisenhauer & Gillessen (2010), and the parameters for the bulge and disk components are from Miyamoto & Nagai (1975).

Assuming an equation of state for ideal gas, the gas pressure,


where is gas density and is sound speed. We consider the spatial dependence of temperature () but assume that it is kept constant with time in each grid cell; the gas is assumed to be locally isothermal. Since we treat the Galactic bulge and disk by global simulation, the “sound speed” here more or less reflects the velocity dispersion of gas clouds rather than the actual temperature of each cloud. In the disk region, we assume is comparable to the observed velocity dispersion (e.g., Bovy et al., 2012),


whereas this value is much smaller than the rotation speed km s. In the bulge region, however, large velocity dispersion, which is probably because of the noncircular motion, is obtained (e.g., Kent, 1992). To mimic this larger velocity dispersion, we adopt


where is the rotation speed of the stellar component.

We smoothly connect these two components to fix the radial dependence:


where we assume the boundary between the bulge and the disk, kpc, with a smoothing width, kpc.

We start our simulation from an equilibrium configuration, in which the gravity is balanced with the pressure-gradient and centrifugal forces:




We assume the initial gas density at the Galactic plane is proportional to the stellar density that fixes the gravitational potential, namely . Then, we determine the distribution of the initial density and rotation frequency, , to satisfy Equations (6) and (7).

Figure 1 presents the radial distribution of the initial equilibrium profile at the Galactic plane. The rotation speed of the gas component shows a bump near kpc, associated with a peak of the sound speed there (top panel). In this region, because of the high temperature, the outward pressure-gradient force is not negligible, compared to the centrifugal force; the inward gravity is balanced by both the pressure-gradient force and the centrifugal force. To satisfy this radial force balance, the rotation speed should be kept considerably smaller than the rotation speed of the stellar component there to suppress the centrifugal force. This equilibrium configuration plays an important role in the amplification of the magnetic field as will be discussed later. It should also be noted that the mass distribution of stars in the inner 0.1 kpc has recently been obtained (Sofue, 2013; Feldmeier et al., 2014), which shows the existence of a certain amount of mass () there. We need to incorporate this contribution in future studies, which would give moderately faster rotation speed in the inner bulge.

The bottom panel of Figure 1 presents the initial density, the gas pressure, and the magnetic pressure. The density is expressed in units of number density, cm, with mean molecular weight, , where is the mass of a hydrogen atom. In this work we adopt the one-fluid approximation and do not distinguish ions, neutral atoms, and molecules. with approximately corresponds to number density in units of hydrogen atoms. For rough estimates of molecular number density, which is dominated by H, in dense clouds, we can use . Initially we set up weak vertical magnetic field with,


in the region of kpc. A reason why we start from the weak magnetic field is that we would like to avoid the direct effect of the initial magnetic field on the dynamics of the gas. The initial magnetic pressure, , is less than times the gas pressure as shown in the bottom panel of Figure 1, namely the plasma .

(0.01,60) () () 384 128 256
Table 2: Simulation Box & Numerical Resolution
Figure 2: Time evolution of the magnetic field strength averaged over the azimuthal and vertical directions (see text) at different radial locations, kpc (solid), 0.5 kpc (dashed), and 2.0 kpc (dotted).
Figure 3: Snapshot views of the bulge region ( kpc) of the simulation at Myr (upper panel). Colours indicate density in units of number density, cm, in logarithmic scale, white lines denote magnetic field lines. In the lower panel the central region is zoomed in and seen from a nearly edge-on angle to show -shaped field lines which are typical for Parker instability. The movies are available (see text and footnote).

We update , , and with time by solving ideal MHD equations,




under the fixed gravitational potential (Equation 1). The numerical scheme we adopt is the second-order Godunov-CMoCCT method (Sano, Inutsuka & Miyama, 1999), in which we solve nonlinear Riemann problems with magnetic pressure at cell boundaries for compressive waves and adopt the consistent method of characteristics (CMoC) for the evolution of magnetic fields (Clarke, 1996; Stone & Norman, 1992) under the constrained transport (CT) scheme (Evans & Hawley, 1988). Since we assume the locally isothermal equation of state (Equation 2), we do not solve an energy equation. We perform our simulation in spherical coordinates, , instead of cylindrical coordinate, , to resolve the central region by fine-scale grids. The size of the simulation box and the number of the grid points are summarized in Table 2. In the direction, the simulation domain covers from the Galactic plane; it does not cover the regions near the poles, and we need to set up an appropriate condition for the boundaries, which is described below. The size of the radial grid, , is enlarged in proportion to . We analyze the numerical data mainly in the cylindrical coordinates by converting the primitive data in the spherical coordinates.

We adopt the same boundary condition as in Suzuki & Inutsuka (2014). We prescribe the outgoing condition for both mass and waves at the boundaries, which was originally developed for simulations for the solar wind (Suzuki & Inutsuka, 2005, 2006), and the accretion condition at both the outer and inner boundaries. Because of the outgoing boundaries at the surfaces, the mass rapidly streams out by vertical outflows driven by effective magneto-turbulent pressure (Suzuki & Inutsuka, 2009, 2014). Therefore, the total mass in the simulation box does not conserve but decreases with time. This process probably overestimates the mass loss especially in the bulge region because our simulation box does not cover up to the sufficiently higher altitude (Suzuki, Muto & Inutsuka, 2010; Fromang et al., 2013). In order to compensate the rapid mass decrease in the bulge, we use very large “density floor”, of 0.8 times the local initial value, ; if density becomes less than the floor value, , at each location, we recover it to . The density floor is utilized in the only regions near the surfaces where the mass is lost by vertical outflows. The density there is lower than that near the Galactic plane, and the effect of mass gain by this treatment is small and it only partially compensates the mass lost by vertical outflows.

When we carry out the simulation, all the physical variables are non-dimensionalized by unit mass, , unit length, kpc, and unit time, Myr. We perform the simulation up to Myr.

3 Results

Figure 4: Face-on views of density in units of cm (colour) and velocity field (arrows) at the Galactic plane at different times, Myr (left) and 441.61 Myr (right).
Figure 5: diagrams at different times, Myr (left) and 441.61 Myr (right). The grid points with high-density regions, cm, in kpc are used to derive the column density in units of cm (colour) by integrating along the direction of line of sight. In the horizontal axis, both Galactic longitude, degree, (bottom axis) and the corresponding transverse distance, kpc, at kpc (top axis) are shown, where the solar system is located at kpc.
Figure 6: Same as Figure 5 but for a large range. Regions with smaller density, cm, are considered to derive the column density than for Figure 5.

3.1 Overview

After the simulation starts, the amplification of the weak vertical seed magnetic field is initially triggered by magneto-rotational instability (MRI, hereafter; Velikhov, 1959; Chandrasekhar, 1961; Balbus & Hawley, 1991) and the vertical shear of the initial rotation profile (Suzuki & Inutsuka, 2014; McNally & Pessah, 2014) to generate the radial and toroidal components. The wavelength of the most unstable mode of the MRI, which is captured by one grid scale, is slightly underresolved at first because the initial field strength is weak. However, the wavelength of the most unstable mode, which is proportional to the field strength, becomes longer with the amplification of magnetic field, and it can be well resolved in the saturated state.The toroidal component is further amplified by the winding due to the radial differential rotation. Furthermore, magnetic buoyancy (Parker, 1966) also plays a role in amplifying the vertical component (Nishikori, Machida & Matsumoto, 2006; Machida et al., 2013).

Figure 2 presents the time evolution of the magnetic field strength averaged over azimuthal () and vertical () directions,


at different radial locations, kpc (solid), 0.5 kpc (dashed), and 2.0 kpc (dotted), where the vertical integration is taken from to kpc. The growth of magnetic field is faster at smaller because the growth rate of the MRI and the winding is scaled by rotation frequency, (see §3.3), which decreases with . Figure 2 shows that the field strength saturates after Myr at and 0.5 kpc in the bulge region, while it is still gradually growing at kpc until the end of the simulation, Myr.

Figure 3 presents 3D snapshots 4 of the bulge region, kpc, after the MHD turbulence is well developed at Myr. The figure exhibits turbulent gas and complicatedly tangled magnetic field lines in the bulge region. One can see turbulent poloidal field lines excited by MRI and Parker instability, whereas the toroidal component slightly dominates the poloidal (radial and vertical) component as will be discussed in §3.3. In the lower panel that zooms in the central region, we pick up -shaped field lines excited by Parker instability, which are actually observed in the Galactic center region (Fukui et al., 2006). This turbulent state is sustained in a quasi-steady manner after the field strength is saturated. In Figure 2 we only show the evolution of the average field strength, and then one cannot see evolutionary properties of the magnetic configuration and turbulence. The movie of Figure 3 (see footnote for the link) demonstrates that the MRI and Parker instability continuously drive non-axysmmetric complex magnetic configuration rather than developing coherent axisymmetric structure. Therefore, the system is in a quasi-steady saturated state in terms of magnetic turbulence, after the saturation of the field strength.

Figure 4 shows the velocity field with the density, , at the Galactic plane () at different times 5; the left panel presents the result at the same time as in Figure 3, and the right panel shows the result at 2.6 Myr later than the left panel, where this time difference ( Myr) roughly corresponds to of the rotation time at kpc. These two panels show that radial motion is stochastically excited particularly around kpc. Outward motion () appears to dominate, whereas one can also recognize inflows () in some regions. These behaviors are a result of the magnetic field as will be discussed in more detail later (§3.3).

3.2 diagrams

Figure 5 shows Galactic longitude–velocity () diagrams observed from the LSR at different times that correspond to Figure 4. The solar system is assumed to rotate with 240 km s in the clockwise direction at kpc (Honma et al., 2012) and , namely km s at kpc in the Cartesian coordinates. The contours indicate the column density in units of cm, integrated along the direction of line of sight. We use the grid points with high-density regions where cm, because the high-density regions more or less trace the cool molecular gas, although our simulation treats one-fluid gas and does not distinguish molecules, neutral atoms, or ions.

Parallelogram-like shapes are observed in the central region as a consequence of the excited radial motion shown in Figure 4. Comparing the two panels, the shape changes with time 6 on account of the time-dependent nature of the radial flows. Moreover, one can see that the shape is not symmetric because the distribution of the excited radial flows is not axisymmetric or bisymmetric. Interestingly, the parallelogram shape obtained in the Milky way is not symmetric (Dame, Hartmann & Thaddeus, 2001; Torii et al., 2010a, then, the shape is not a “parallelogram” in a strict sense), which cannot be explained by a simple bar potential.

Figure 6 presents diagrams for a larger degree. We use regions with smaller density, cm, to calculate surface density than for Figure 5 in order to take into account the outer region where the density is smaller. These diagrams also show non-symmetric structure in the plane, which indicates that non-axisymmetric gas distributions are extended to larger .

Figure 7: top: Time evolution of heliocentric velocities at Galactic longitude, degree The solid lines indicate the maximum and minimum velocities that give column density, cm. The dashed lines trace the radial motion of the front-side () and back-side () at kpc ( kpc in the Cartesian coordinates). bottom: Time evolution of the slope of the “parallelogram” near degree. The solid and dashed lines indicate the slope of the top and bottom sides, respectively. The slope is derived from the average between degrees, which corresponds to pc around the Galactic center. The shaded regions indicate the observational values, the region centered at 20 km sdeg for the top side of the parallelogram and the region centered at 55 km sdeg for the bottom side.
Figure 8: Radial distribution of the plasma (top) and the value (bottom; solid line) which is the sum of the Maxwell (dashed line) and Reynolds (dotted line) stresses. The data are averaged over the azimuthal direction in full rotation (), over the vertical direction within kpc, and the time Myr.
Figure 9: Radial distribution of various quantities. The data are averaged over the azimuthal direction in full rotation (), over the vertical direction within kpc (except for panel (a)), and the time –485.21 Myr. (see text). (a): The strength of differential rotation, at the Galactic plane (solid line); The average is not taken over the component but over the component and time. The dashed line represents the initial condition. Rigid-body rotation (), flat rotation (), and Keplerian rotation () are shown by dotted lines for reference. (b): (solid), (dashed), and (dotted) components of root-mean squared (rms) magnetic field in units of G. The dotted line with a constant slope denotes the initial vertical magnetic field. (c): Root mean squared radial velocity, (solid), and mean radial velocity, (dashed and dotted lines for positive and negative values).
Figure 10: Simplified schematic picture of the mechanism that drives radial flows. Let us suppose three rings, ‘1’, ‘2’, and ‘3’, that rotate with rotation frequency, , , and , at the outer edge of each ring. If , differential rotation is more intense in the red region (Region 2). The magnetic field is amplified more effectively there by the MRI and the field-line stretching, and the Maxwell stress, , (Equation 17) is larger than those in the neighboring zones, which gives the faster outward transport of the angular momentum in the red region. As a result, decreases and increases. From the radial force balance, the inner edge of the red region moves inward because of the decrease of the centrifugal force, while the outer edge moves outward by the increase of the centrifugal force. Although the inhomogeneity of angular momentum transport is distributed in a far more complicated fashion in our simulation, the essential point can be explained by this simple picture.
Figure 11: Additional mechanism that drives radial flows. The initial setting is the same as in Figure 10. Radial magnetic field, which is triggered by MRI in our simulation, is amplified to general toroidal field by the differential rotation selectively in the red region (Region 2). The magnetic pressure, , of the toroidal field there is larger than in the neighboring regions. Then, radial flows are generated by the magnetic pressure gradient force.

We examine the time-dependency of the parallelogram-like feature in Figure 7. The solid lines in the top panel shows the time-evolution of the maximum and minimum velocities ( & ) for column density cm at degree; they correspond to the positions where the upper and lower sides of the parallelogram cross the vertical line of degree in the diagram. The figure shows that and vary within 50-150 km s with time. Basically and change independently of each other, and then, the shape of the symmetric parallelogram is not always kept but distorted with time as seen in Figure 5.

In Figure 7, the time evolution of the motion of the gas at kpc is also shown (dashed lines) for comparison. These lines indicate that the direction of the radial motion tends to be outward ( for the front side and for the back side), whereas inward motions are also seen occasionally.

The bottom panel presents the slope of the parallelogram measured near , which is derived from the average in the region with deg, corresponding to pc. The solid and dashed lines denote the top and bottom sides of the parallelogram, which mostly reflect far and near sides with respect to the Galactic center, respectively. The top and bottom sides change independently because of the non-axisymmetric excitation of radial flows. The slopes vary from 15 to 90 km sdeg, which covers the observed slope, km sdeg, for the top side of the parallelogram and km sdeg for the bottom side (Torii et al., 2010a).

3.3 Origin of Noncircular Motion –Roles of Magnetic Field–

Our MHD simulation shows that noncircular motions are excited, even though the axisymmetric gravitational potential is adopted. Here we discuss the roles of the magnetic field in driving the radial gas flows by inspecting the numerical data. In Figures 8 and 9 we present the radial profile of various quantities averaged over azimuthal () and vertical () directions and over time . We take the time-average of the final stage of the simulation from Myr to Myr. Although the magnetic field is already amplified to achieve the saturated state in the bulge region, it is still developing outside the bulge (Figure 2). Therefore, in these figures we focus on the inner region, kpc. We take the average of a physical quantity, , as follows:


where we set and kpc. For the quantities concerning velocity, we adopt density-weighted averages; the average of flow velocity and root-mean squared (rms) velocity is taken by




while the simple average is taken for the quantities on magnetic field.

Figure 8 presents the plasma value,


which is the ratio of gas pressure to magnetic pressure (top panel), and the value (bottom panel; solid line) (Shakura & Sunyaev, 1973),


which is the sum of Reynolds (; dotted line) and Maxwell (; dashed line) stresses. The fluctuation component, , is derived by subtracting the mean rotation,


The top panel of Figure 8 shows that the plasma is below 100 in the bulge region, kpc, which is much smaller than the initial value, , as a consequence of the amplification of the initial weak magnetic field. The bottom panel shows that a large value is obtained. In the inner bulge, kpc, the Maxwell stress dominates the Reynolds stress. In the outer bulge, the Reynolds stress exceeds the Maxwell stress, and they both show bumpy structures, which reflects the complicated rotation profile as discussed from now.

Figure 9 (a) presents the time and –averaged strength of the radial differential rotation, , at the Galactic plane. Note that correspond to rigid-body, flat, and Keplerian rotations, respectively. The MRI is active when , and its maximum growth rate to generate radial magnetic field is given by


(Balbus & Hawley, 1991). The amplification of toroidal field by winding in the radial differential rotation follows (e.g., Balbus & Hawley, 1998)


Equations (19) and (20) show that stronger differential rotation amplifies magnetic field faster by both MRI and field-line stretching.

Figure 9(a) shows that the initial rotation profile is more or less kept at the later time. The local minimum of is located at kpc, where the rotation rapidly decreases with because of the contribution from the gas pressure-gradient force (§2). In this region, the magnetic field is effectively amplified on account of the strong differential rotation. The magnetic field strength is amplified to mG in the bulge region, which is consistent with an observational lower limit, G (Crocker et al., 2010) and an empirical estimate, mG, based on multiple observations (Ferrière, 2009).

The effective amplification of the magnetic field around the location for the maximum differential rotation at kpc gives a bump in the Maxwell stress, (Equation 17) as shown by the dashed line in the bottom panel of Figure 8. determines the outward transport of angular momentum (e.g., Lynden-Bell & Pringle, 1974) by magnetic field. The inhomogeneous distribution of indicates that spatially dependent gain or loss of the angular momentum takes place, which is illustrated by a simplified cartoon in Figure 10; in a region with the gain of angular momentum the rotation is accelerated, and vice versa. This leads directly to the increase or decrease of the centrifugal force, which excites radial motion. Although what is taking place in our simulation is more complicated, namely is distributed in a non-axisymmetric manner, and both outward and inward flows are generated even at the same as exhibited in Figure 4, this simple conceptual cartoon in Figure 10 explains the essential point.

In addition to the inhomogeneous transport of angular momentum, magnetic pressure also contributes to the excitation of radial flows. In the outer bulge region near kpc, the differential rotation is weak and the amplification of magnetic field is suppressed there. Figure 9(b) exhibits the rapid decrease of the magnetic field strength with in kpc. The rapid decrease of causes the radial pressure gradient, , which drives outward radial flows (Figure 11).

We have discussed the two types of the processes, the inhomogeneous angular momentum transport and the magnetic pressure-gradient force, that generate radial flows. By inspecting each term in the momentum equation (10) numerically, we found that the former dominates the latter by .

Figure 9(c) shows that the direction of the mean flow is outward (dashed line ) in kpc, which is expected from both processes. This radially outward flow also transports the angular momentum, which is a reason in part why the Reynolds stress (the bottom panel of Figure 8) shows a bump in kpc. Compared to the mean flow, the rms velocity, , (solid line) at kpc is quite large, km s, because the MRI-triggered turbulence drives both inward and outward intermittently (Sano & Inutsuka, 2001, see also movies of Figures 3 and 4), as discussed so far. This is the main reason why the simulated diagrams show a thick parallelogram-like shape (Figure 5).

3.4 Non-axisymmetric Structure

Figure 12: Density structure in the Galactic center region at Myr. (Zoomed-in view of the left panel of Figure 3.)

As we have discussed in the previous subsection, the velocity field shows non-axisymmetric radial flows because of the intermittency of the MRI-triggered turbulence. The density distribution also shows non-axisymmetric inhomogeneity, which is typical for the MRI turbulence as well (Suzuki & Inutsuka, 2014), in the central region as exhibited in Figure 12. One can see multiple arm-like structure, and when measured at pc, the density varies largely within, cm, which roughly corresponds to cm in molecular number density.

It is well known that the Milky way possesses the central molecular zone (CMZ), which contains molecular gas with mass and occupies a volumetric filling factor within pc (Morris & Serabyn, 1996). The CMZ consists of non-axisymmetric clouds with complicated structure (Oka et al., 2005; Kruijssen, Dale & Longmore, 2015). Our simulation implies that such non-axisymmetric clouds are a natural outcome of magnetic activity in the Galactic center. The molecular number density estimated from our simulation above is smaller than typical observed values, cm (e.g., Nagayama et al., 2007). This is because our simulation assumes the one-fluid gas with the locally isothermal equation of state. In pc, the temperature is assumed to be constant K in the simulation, which is much higher than the observed typical temperature, 30–200 K, of molecular clouds in the CMZ (e.g. Morris & Serabyn, 1996). We need to treat the energy equation with heating and cooling to deal with different phases of the gas material such as molecular clouds and warm neutral media that respectively possess different filling factors. The temperature of dense regions in the simulation would be much lower if the radiative cooling was taken into account, and accordingly the density would be higher to satisfy the pressure balance.

3.5 Outflows

Figure 13: Vertical cut at of (colours) and velocity field (arrow) at Myr.

Figure 13 presents the snapshot of a vertical cut at Myr, Here the colours indicate ; redder colours correspond to regions with relatively strong magnetic pressure compared to gas pressure. The figure shows that flow patterns are quite complicated; one can see both inward flows to the bulge and outward streams from it. Although both inward and outward streams coexist, if we integrate the mass flux in space and time, the direction of the net gas flow is outward from the bulge; for example the mass loss rate measured at the spherical surface, kpc, averaged during Myr, gives 70 yr, and the kinetic energy luminosity, erg s. It should be noted that this outflow rate might be suppressed when a stellar-bar potential was taken into account (Chandran, Cowley & Morris, 2000).

In addition to radial flows, upward flows (winds) stream out of the surfaces of the simulation box. The typical speed of these flows is several tens km s, and in some regions the speed exceeds 100 km s. Such vertical outflows possibly explain the observed vertical structure of molecular gas associated with the DHN (Enokiya et al., 2014; Torii et al., 2014b). Recent observation further reveals CO emission at high Galactic latitude, , including several filamentary features in addition to diffuse extended halo-like CO gas up to 2 in above and below the CMZ (Torii et al., 2014a). A plausible interpretation is that these high- features are also driven by the Parker instability; although they are filamentary without a clear loop-like shape having two foot points, the MHD numerical simulations by Machida et al. (2009) show that magnetic-flotation loops are generally not symmetric and can be seen as an open single filament depending on the ambient density/field distribution and evolutionary effect. Further signature is associated with the double helix nebula (DHN, hereafter) toward (Morris, Uchida & Do, 2006), where a column of molecular gas of 200-pc height is found to be associated with the DHN which is nearly vertical to the plane (Enokiya et al., 2014; Torii et al., 2014b). Such a feature is explained as created by a magnetic tower formed above and below the central black hole according to preceding numerical simulations (Kato, Mineshige & Shibata, 2004; Machida et al., 2009).

Recently, the existence of outflows or winds in the Galactic center region is extensively discussed (Everett et al., 2008; Crocker et al., 2011), and it might be closely linked to the formation mechanisms of the Fermi bubbles via magnetic activity (Carretti et al., 2013). Our simulation is not directly applicable to the Fermi bubbles since the simulation does not cover the region near the polar axis. However, the above rough estimate of the energetics of the outflow in our simulation shows that magnetic activity is a reliable mechanism in driving outflows from the bulge.

4 Summary and Discussion

We have investigated the excitation of radial gas flows in the Galactic center region by the 3D global MHD simulation with the axisymmetric gravitational potential. The initial weak vertical magnetic field is amplified by the MRI at the beginning, and in addition eventually by the Parker instability and the field-line stretching owing to the differential rotation. In the bulge region, the final field strength is mG, and the plasma . Because of the gas pressure-gradient force, the equilibrium rotation frequency is non-monotonic with radial distance. The amplification of the magnetic field is systematically more effective in regions with stronger differential rotation. Then, the transport of the angular momentum is spatially dependent, and the acceleration or deceleration of the rotation speed occurs depending on regions. This breaks the radial pressure balance owing to the change of the centrifugal force, and excites radial flows of the gas in a time-dependent and non-axisymmetric manner. In addition, the radial component of the magnetic pressure gradient is produced, which also excites radial flows. As a result, the simulated diagram exhibits a time-dependent and asymmetric parallelogram shape. The rotation curve near the Galactic center exhibits complicated features (Sofue, 2013). Our simulation implies that stochastically excited radial motion might be contaminated in these features.

In order to interpret the parallelogram in the CMZ, the elliptical orbit in the stellar-bar potential has been used as the only viable model (Binney et al., 1991). The model may not be fully appropriate as the interpretation from an observational point. First, the bar potential alone is not able to create high- features up to 200 pc above the plane; these “loops” are better explained by the Parker instability (Fukui et al., 2006). There is no other driving mechanism known to expel gas up to that height by only stellar gravity. Second, the parallelogram is not symmetric in the - diagram; the velocity gradient in the positive velocity, 20 km s deg, is significantly different by a factor of 3 from that in the negative velocity, 55 km s deg (Torii et al., 2010a). The present model shows that the velocity gradient of the parallelogram is naturally produced as a time dependent manner, whereas it is not clear how such a difference is explained in the bar potential model.

On the other hand, we see some features that are not reproduced in the present simulation. The large broad features like the clump 2 and the 5.5 deg feature are not reproduced in the present simulation. In our simulation we have assumed one-fluid gas and neglected heating and cooling with a simplified treatment of the locally isothermal gas. In order to treat dense molecular gas in a quantitative fashion, we should replace these simplifications. If radiative cooling, as well as adiabatic heating and cooling, is included, large density contrast between cool molecular gas and warm/hot gas will be naturally reproduced, which may explain the observed fine-scale features.

We have shown that the magnetic activity possibly produces the overall trend of the observed asymmetric parallelogram-shape in the diagram even though the axisymmetric gravitational potential is considered. It is also desirable to test quantitatively the importance of the magnetic activity in contrast to role of the stellar bar potential in the noncircular motion of the gas after the cooling/heating effect is taken into account in our simulation.

Although in this paper we have focused on the bulge region of our Milky Way, results of our simulation can be applied to other galaxies. For example, Sakamoto et al. (2006) observed NGC 253 by the Submillimeter Array and reported that they found an expanding circumnuclear disk with km s. The obtained position velocity diagrams show asymmetric parallelogram features in different wavelengths. These features are qualitatively similar to but quantitatively different from those obtained in the Milky Way. A possible explanation of the difference is the time variability inhering in the MHD turbulence, based on our simulation.


This work was supported by Grants-in-Aid for Scientific Research from the MEXT of Japan, 24224005 (PI: YF). Numerical simulations in this work were carried out at the Cray XC30 (ATERUI) operated in CfCA, National Astrophysical Observatory of Japan, and the Yukawa Institute Computer Facility, SR16000. The authors thank the anonymous referee for constructive comments to improve the paper. TKS thanks Prof. Shu-ichiro Inutsuka and Dr. Kazunari Iwasaki for fruitful discussion.



  1. pagerange: Stochastic Noncircular Motion and Outflows Driven by Magnetic Activity in the Galactic Bulge RegionLABEL:lastpage
  2. pubyear: 2015
  3. Infrared observation of the Galactic center region was also independently done by Okuda et al. (1977).
  4. Movie of the simulation is available at
  5. Movie is available at
  6. Movie is available at


  1. Baba J., Saitoh T. R., Wada K., 2010, PASJ, 62, 1413
  2. Balbus S. A., Hawley J. F., 1991, ApJ, 376, 214
  3. Balbus S. A., Hawley J. F., 1998, Reviews of Modern Physics, 70, 1
  4. Bally J., Stark A. A., Wilson R. W., Henkel C., 1987, ApJS, 65, 13
  5. Bania T. M., 1977, ApJ, 216, 381
  6. Binney J., Gerhard O. E., Stark A. A., Bally J., Uchida K. I., 1991, MNRAS, 252, 210
  7. Blitz L., 1994, in Astronomical Society of the Pacific Conference Series, Vol. 66, Physics of the Gaseous and Stellar Disks of the Galaxy, King I. R., ed., p. 1
  8. Blitz L., Spergel D. N., 1991, ApJ, 379, 631
  9. Bovy J. et al., 2012, ApJ, 759, 131
  10. Carretti E. et al., 2013, Nature, 493, 66
  11. Chandran B. D. G., Cowley S. C., Morris M., 2000, ApJ, 528, 723
  12. Chandrasekhar S., 1961, Hydrodynamic and hydromagnetic stability. Oxford: Clarendon
  13. Chuss D. T., Davidson J. A., Dotson J. L., Dowell C. D., Hildebrand R. H., Novak G., Vaillancourt J. E., 2003, ApJ, 599, 1116
  14. Clarke D. A., 1996, ApJ, 457, 291
  15. Crocker R. M., Jones D. I., Aharonian F., Law C. J., Melia F., Ott J., 2011, MNRAS, 411, L11
  16. Crocker R. M., Jones D. I., Melia F., Ott J., Protheroe R. J., 2010, Nature, 463, 65
  17. Dame T. M., Hartmann D., Thaddeus P., 2001, ApJ, 547, 792
  18. de Vaucouleurs G., 1964, in IAU Symposium, Vol. 20, The Galaxy and the Magellanic Clouds, Kerr F. J., ed., p. 195
  19. Enokiya R. et al., 2014, ApJ, 780, 72
  20. Evans C. R., Hawley J. F., 1988, ApJ, 332, 659
  21. Everett J. E., Zweibel E. G., Benjamin R. A., McCammon D., Rocks L., Gallagher, III J. S., 2008, ApJ, 674, 258
  22. Feldmeier A. et al., 2014, A&A, 570, A2
  23. Ferrière K., 2009, A&A, 505, 1183
  24. Fromang S., Latter H., Lesur G., Ogilvie G. I., 2013, A&A, 552, A71
  25. Fujishita M. et al., 2009, PASJ, 61, 1039
  26. Fukui Y. et al., 2006, Science, 314, 106
  27. Genzel R., Eisenhauer F., Gillessen S., 2010, Reviews of Modern Physics, 82, 3121
  28. Hayakawa S., Matsumoto T., Murakami H., Uyama K., Thomas J. A., Yamagami T., 1981, A&A, 100, 116
  29. Honma M. et al., 2012, PASJ, 64, 136
  30. Kato Y., Mineshige S., Shibata K., 2004, ApJ, 605, 307
  31. Kent S. M., 1992, ApJ, 387, 181
  32. Kim W.-T., Stone J. M., 2012, ApJ, 751, 124
  33. Koda J., Wada K., 2002, A&A, 396, 867
  34. Kruijssen J. M. D., Dale J. E., Longmore S. N., 2015, MNRAS, 447, 1059
  35. Kudo N. et al., 2011, PASJ, 63, 171
  36. Liszt H. S., Burton W. B., 1978, ApJ, 226, 790
  37. Liszt H. S., Burton W. B., 1980, ApJ, 236, 779
  38. Lynden-Bell D., Pringle J. E., 1974, MNRAS, 168, 603
  39. Machida M. et al., 2009, PASJ, 61, 411
  40. Machida M., Nakamura K. E., Kudoh T., Akahori T., Sofue Y., Matsumoto R., 2013, ApJ, 764, 81
  41. Matsumoto T., Hayakawa S., Koizumi H., Murakami H., Uyama K., Yamagami T., Thomas J. A., 1982, in American Institute of Physics Conference Series, Vol. 83, The Galactic Center, Riegler G. R., Blandford R. D., eds., pp. 48–52
  42. McNally C. P., Pessah M. E., 2014, ArXiv e-prints
  43. Miyamoto M., Nagai R., 1975, PASJ, 27, 533
  44. Molinari S. et al., 2011, ApJ, 735, L33
  45. Morris M., Davidson J. A., Werner M., Dotson J., Figer D. F., Hildebrand R., Novak G., Platt S., 1992, ApJ, 399, L63
  46. Morris M., Serabyn E., 1996, ARA&A, 34, 645
  47. Morris M., Uchida K., Do T., 2006, Nature, 440, 308
  48. Morris M. R., 2014, ArXiv e-prints
  49. Nagayama T., Omodaka T., Handa T., Iahak H. B. H., Sawada T., Miyaji T., Koyama Y., 2007, PASJ, 59, 869
  50. Nishikori H., Machida M., Matsumoto R., 2006, ApJ, 641, 862
  51. Nishiyama S. et al., 2010, ApJ, 722, L23
  52. Oka T., Geballe T. R., Goto M., Usuda T., McCall B. J., 2005, ApJ, 632, 882
  53. Oka T., Hasegawa T., Sato F., Tsuboi M., Miyazaki A., 1998, ApJS, 118, 455
  54. Okuda H., Maihara T., Oda N., Sugiyama T., 1977, Nature, 265, 515
  55. Parker E. N., 1966, ApJ, 145, 811
  56. Peters, III W. L., 1975, ApJ, 195, 617
  57. Riquelme D., Bronfman L., Mauersberger R., May J., Wilson T. L., 2010, A&A, 523, A45
  58. Rodriguez-Fernandez N. J., Combes F., 2008, A&A, 489, 115
  59. Rougoor G. W., Oort J. H., 1960, Proceedings of the National Academy of Science, 46, 1
  60. Sakamoto K. et al., 2006, ApJ, 636, 685
  61. Sano T., Inutsuka S., Miyama S. M., 1999, in Astrophysics and Space Science Library, Vol. 240, Numerical Astrophysics, Miyama S. M., Tomisaka K., Hanawa T., eds., Boston, MA: Kluwer, p. 383
  62. Sano T., Inutsuka S.-i., 2001, ApJ, 561, L179
  63. Sawada T., Hasegawa T., Handa T., Cohen R. J., 2004, MNRAS, 349, 1167
  64. Scoville N. Z., 1972, ApJ, 175, L127
  65. Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337
  66. Shibata K., Matsumoto R., 1991, Nature, 353, 633
  67. Sofue Y., 2007, PASJ, 59, 189
  68. Sofue Y., 2013, PASJ, 65, 118
  69. Stone J. M., Norman M. L., 1992, ApJS, 80, 791
  70. Suzuki T. K., Inutsuka S.-i., 2005, ApJ, 632, L49
  71. Suzuki T. K., Inutsuka S.-i., 2006, Journal of Geophysical Research (Space Physics), 111, 6101
  72. Suzuki T. K., Inutsuka S.-i., 2009, ApJ, 691, L49
  73. Suzuki T. K., Inutsuka S.-i., 2014, ApJ, 784, 121
  74. Suzuki T. K., Muto T., Inutsuka S.-i., 2010, ApJ, 718, 1289
  75. Takahashi K. et al., 2009, PASJ, 61, 957
  76. Takeuchi T. et al., 2010, PASJ, 62, 557
  77. Torii K., Enokiya R., Fukui Y., Yamamoto H., Kawamura A., Mizuno N., Onishi T., Ogawa H., 2014a, in IAU Symposium, Vol. 303, IAU Symposium, Sjouwerman L. O., Lang C. C., Ott J., eds., pp. 106–108
  78. Torii K., Enokiya R., Morris M. R., Hasegawa K., Kudo N., Fukui Y., 2014b, ApJS, 213, 8
  79. Torii K. et al., 2010a, PASJ, 62, 675
  80. Torii K. et al., 2010b, PASJ, 62, 1307
  81. Tsuboi M., Handa T., Ukita N., 1999, ApJS, 120, 1
  82. Tsuboi M., Inoue M., Handa T., Tabara H., Kato T., Sofue Y., Kaifu N., 1986, AJ, 92, 818
  83. Velikhov E. P., 1959, Zh. Eksp. Teor. Fiz., 36, 1398
  84. Yusef-Zadeh F., Morris M., Chance D., 1984, Nature, 310, 557
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description