MHD Simulation in Galactic Center Region with Radiative Cooling and Heating

We investigate the role of magnetic field on the gas dynamics in a galactic bulge region by three-dimensional simulations with radiative cooling and heating. While a high-temperature corona with T > 106 K is formed in the halo regions, the temperature near the midplane is ≲104 K following the thermal equilibrium curve determined by the radiative cooling and heating. Although the thermal energy of the interstellar gas is lost by radiative cooling, the saturation level of the magnetic field strength does not significantly depend on the radiative cooling and heating. The magnetic field strength is amplified to 10 μG on average and reaches several hundred microgauss locally. We find the formation of magnetically dominated regions at midlatitudes in the case with the radiative cooling and heating, which is not seen in the case without radiative effect. The vertical thickness of the midlatitude regions is 50–150 pc at the radial location of 0.4–0.8 kpc from the Galactic center, which is comparable to the observed vertical distribution of neutral atomic gas. When we take the average of different components of energy density integrated over the galactic bulge region, the magnetic energy is comparable to the thermal energy. We conclude that the magnetic field plays a substantial role in controlling the dynamical and thermal properties of the galactic bulge region.


INTRODUCTION
The neutral atomic hydrogen (H I ) gas in the Galactic bulge (GB hereafter) region (R < 1.0 kpc) has a thickness of over 70 pc in the vertical direction (Sofue 2022a,b).On the other hand, if we assume the force balance, the vertical scale height can be estimated to be h = c s /Ω ≈ 20 pc 1 , where c s ≈ 10 km s −1 is the sound velocity in the gas with temperature ≈ 10 4 K and g z = − ∂ Φ ∂z , is the observed rotational frequency (e.g., So- Corresponding author: Kensuke Kakiuchi kakiuchi@g.ecc.u-tokyo.ac.jp, kakiuchi.mw@gmail.com 1 In a strict sense ∂ gz ∂z should be used for Ω where the gz is vertical component of gravity (see Section 4.1).fue 2013).c s could be smaller if the temperature of H I gas is lower, which gives further smaller h.This comparison indicates that the the vertical distribution of the H I gas cannot be supported only by the gas pressure, suggesting the presence of additional mechanisms required to maintain the observed vertical thickness.
The magnetic field is one of such plausible mechanisms; the vertically extended gas can be explained by the contribution from magnetic pressure in combination with gas pressure.In addition, there have been a variety of complex structures observed that are expected to reflect magnetic-field configurations, such as molecu-lar loops2 (Fukui et al. 2006;Fujishita et al. 2009;Torii et al. 2010a,b), spiral-shaped molecular gas (the Double Helix Nebula; Morris et al. 2006;Enokiya et al. 2014) and non-thermal filamentary structures (Yusef-Zadeh et al. 1984;Morris & Yusef-Zadeh 1985;Anantharamaiah et al. 1991;Lang et al. 1999;LaRosa et al. 2004;Paré et al. 2019;Heywood et al. 2022;Yusef-Zadeh et al. 2022).These observations illustrate that the magnetic field plays a significant role in GB region.Similar magnetic properties are obtained in other spiral galaxies (Golla & Hummel 1994;Sofue et al. 1994;Konishi et al. 2022).
The magnetic field strength in the GB region of the Milky way has been estimated by multiple independent observational techniques.For instance, if the shape of the Radio Arc, the most extended filamentary structures in the GB region, is maintained by magnetic fields in surrounding turbulent medium, the field strength of B = 0.1 − 1 mG is required (Morris & Serabyn 1996).Pillai et al. (2015) observed the polarization of the dark nebula G0.253+0.016(known as 'Brick') and estimated the magnetic field strength to be a few mG using the Davis-Chandrasekhar-Fermi (DCF) method (Davis 1951;Chandrasekhar & Fermi 1953).An upper limit of the nondetected γ-rays by EGRET gives a lower limit on the average magnetic field strength ≥ 0.05 mG within 400 pc of the Galactic center region (Crocker et al. 2010); the magnetic field needs to be sufficiently strong in order that cosmic-ray electrons effectively lose their energy via Synchrotron radiation to keep the inverse Compton radiation less than the detection limit.On the basis of these estimates, the energy density of the magnetic field equals to or even exceeds the thermal energy density of the gas in the GB region.Such strong magnetic fields are also reported in the central region of other spiral galaxies (e.g., Beck et al. 2005;Tabatabaei et al. 2018).Thus, it is essential to take into account magnetohydrodynamical (MHD) effects in order to understand the properties of the gas in the GB region of spiral galaxies.
It is also inferred that magnetic fields may affect star formation activity.While there is molecular cloud complex called the Central Molecular Zone (CMZ) in the Milky way, the star formation rate, ∼ 0.02-0.1 M ⊙ yr −1 (e.g.Yusef-Zadeh et al. 2009;Immer et al. 2012), integrated over the GB region is much smaller than that estimated by extrapolation from the Galactic disk (Longmore et al. 2013;Kruijssen et al. 2014).The reduced star formation may be a consequence of strong magnetic activity (Morris 1993;Chuss et al. 2003;Henshaw et al. 2022) Numerical investigations on the effect of magnetic fields in the dynamics of the GB region have also been carried out.For example, global MHD simulations by Machida et al. (2009) demonstrated that gas is lifted up from the mid-plane by buoyantly rasing magnetic loops, which is a possible mechanism for the observed molecular loops (Fukui et al. 2006).Suzuki et al. (2015) reported with a global MHD simulation that the observed parallelogram on the position-velocity diagram can be reproduced by radial flows excited by a non-smooth radial profile of rotational velocity.In the same simulation, "sliding slopes" associated with magnetic loops are also ubiquitously found in a time-dependent manner, which drive high-speed downflows in excess of 100 km s −1 (Kakiuchi et al. 2018).These numerical simulations treat only high-temperature gas with T ∼ 10 6 K to avoid numerical difficulty in the description of lowtemperature fine-scale clouds.However, it is essentially important to handle lower-temperature gas with radiative cooling in order to model the realistic interstellar condition.In this paper, we investigate the magnetic activity in the GB region of spiral galaxies in detail, taking into account the thermal evolution including radiative cooling and heating in the numerical simulation code of Suzuki et al. (2015) 3 .

NUMERICAL METHOD
We perform three dimensional (3D) global MHD simulations to study the time evolution of the gas disk in the GB region with an initially weak vertical magnetic field.Our simulations are performed in a cylindrical coordinate system (R, ϕ, z).In this work, we adopt an axisymmetric gravitational potential.We will present an explanation of the numerical method in the following order: basic equations (Section 2.1), the gravitational potential (Section 2.2), radiative cooling and heating (Section 2.3), and numerical setup including initial and boundary conditions (Section 2.4).

Basic equations
We solve the time evolution of the gas with the following ideal MHD equations with the second-order Godunov CMoC-CT method (Evans & Hawley 1988;Clarke 1996;Sano & Miyama 1999;Suzuki & Inutsuka 2014;Suzuki et al. 2015) including radiative cooling and heating: The cooling and heating rates against temperature.The solid black line denotes the cooling rate for n = 1.0 cm −3 and the dash black lines denote the heating function with G(R) = 10 3 G0 (in the GB region R < 1.0 kpc) and G0 (at the solar neighborhood) in our simulation.The black dots plot raw data by Sutherland & Dopita (1993) and the gray dotted line plots the fitting function by Koyama & Inutsuka (2002).(b) The mean molecular weight adopted in our simulation.
We assume an equation of state for ideal gas, where γ = 5/3 is the ratio of specific heats.The thermal energy per mass e is expressed as e = k b T (γ−1)µmH where T is temperature, k b is the Bolzman constant, m H is the unified atomic mass unit.We adopt the following approximated formula for the mean molecular weight, µ, that depends on temperature, in order to take into account the effect of the ionization, which is presented in Figure 1(b).Φ in equation ( 2) is the gravitational potential, which will be described in Section 2.2.We do not consider the self-gravity of the gas.L in equation ( 4) represents the net cooling rate [erg s −1 ] as a function of number density n = ρ/µm H and temperature T : where Λ [erg s −1 cm 3 ] is cooling efficiency and Γ [erg s −1 ] is heating rate.Figure 1 (a) shows the cooling and heating functions for temperature, which will be described in Section 2.3.

The gravitational potential
We perform our MHD simulations under timeindependent and axisymmetric Galactic gravitational potential.This potential consists of four components of gravitational sources, the super-massive black hole (SMBH) at the Galactic center, the stellar bulge, the stellar disk, and the dark matter halo; We assume that the SMBH as a point mass of 4.4 × 10 6 M ⊙ (Genzel et al. 2010).The gravitational potential of the stellar bulge (i = 1) and disk (i = 2) is adopted from Miyamoto & Nagai (1975); where G is the gravitational constant.The fit parameters are set to be M 1 = 2.05 × 10 10 M ⊙ , M 2 = 2.57 × 10 11 M ⊙ , a 1 = 0.0, a 2 = 7.258, b 1 = 0.495 and b 2 = 0.52, respectively.For the dark matter halo, we adopt the NFW model (Navarro et al. 1996;Sofue 2013): where r = √ R 2 + z 2 is the spherical radius.We assume r h = 10.7 kpc and ρ h,0 = 1.82 × 10 −2 M ⊙ pc −3 .

Cooling and heating
We adopt the radiative cooling of the solar-metallicity gas for simplicity, whereas the metallicity gradually increases toward the center to give a few times greater value in the GB region than in the solar neighborhood (Shields & Ferland 1994;Maeda et al. 2002).We cover the temperature range of 10 3 K < T < 7 × 10 7 K.For the non/weakly ionized gas with T < 10 4 K, we adopt the cooling function Λ l for the interstellar medium by Koyama & Inutsuka (2002), which is shown by the gray dashed line in Figure 1(a).For the (almost) fully ionized gas with T > 10 4 K, we prescribe the cooling function Λ h for optically thin plasma by Sutherland & Dopita (1993).A fitting formula for Λ h is presented in Appendix A. We smoothly connect Λ l and Λ h via with where T b = 10 4 K, log 10 (∆T b ) = 0.1.The first term of equation ( 7), nΛ, for n = 1 cm −3 is shown as the black solid line in Figure 1.
We do not solve low-temperature (T < 10 3 K) clouds, because our simulations focus on the global properties of the GB region.However, these should be taken into account in more elaborated simulations to study star formation activity because they contain large mass in a small volume.Their scale height is typically ≲ several pc, which is comparable to the minimum resolution, a few pc, in the current simulation setup; higherresolution is required to capture such dense cool clouds.We discuss this topic in Section 4.2.
The heating rate from the interstellar radiation field (ISRF) depends on the surrounding stellar environment, such as the OB association.Therefore, the radiation heating rate is relatively high near the Galactic center where the stellar density is higher than in the solar neighborhood (Wolfire et al. 2003).The ISRF first heats up dust grains.The gas is eventually heated through collisions with the dust grains.We here assume that the dust and gas components are thermally well coupled via sufficiently frequent collisions.The temperature is set by the balance between the heating due to the absorption of photons from the ISRF and the cooling due to thermal radiation from the dust particles.
The ISRF intensity is given as a local profile, G(R), relative to the standard value, G 0 , calibrated at the solar neighborhood (Draine 1978).The ISRF is empirically estimated to be 1000 G 0 (e.g., Clark et al. 2013); this value reproduces the observed thermal properties  The heating rate due to cosmic-rays is scaled to the ionization rate by cosmic-rays ξ CR .In the solar neighborhood, ξ CR = 3.0 × 10 −17 s −1 .As with the ISRF, the ionization rate near the Galactic center is higher than one at the solar neighborhood because of the increase in the flux of cosmic-rays.The ionization rate is estimated to be ξ CR ∼ 2.0 × 10 −14 s −1 in the inner several 100 parsec (Oka et al. 2019), higher than in the solar neighborhood by a factors of ∼ 10 3 .Then, we assume that the ionization rate, and accordingly the heating rate by cosmic-rays, is also simply proportional to G(R) as in the ISRF.
We take account of the heating from ISRF and cosmicrays, and model the function of heating rate in our simulations as follows: where Γ 0 = 2.0 × 10 −26 erg s −1 (Koyama & Inutsuka 2002) is the heating rate at the solar neighborhood.The heating from the ISRF is ineffective above a certain temperature because the intensity of the UV radiation is mainly determined by massive stars.We set T cut = 5.0 × 10 4 K, referring to the effective temperature of O-type stars.In reality, the cosmic-ray heating is still active in the higher-temperature regime.However, we ignore the heating above T cut for simplicity, because the bulk of the heating is dominated by the UV radiation from massive stars.Heating functions, Γ, for G(R) = 10 3 G 0 in R < 1 kpc (long dashed line) and G(R) = G 0 at the solar neighborhood (short dashed line) are compared in Figure 1.
The temperature of the gas in thermal equilibrium can be derived from L = 0 when the heating is balanced with the cooling.Figure 3 shows thermal equilibrium curves in an n − T diagram in R < 1.0 kpc (solid) and at the solar neighborhood (dashed).

Numerical Setup
The initial gas distribution satisfies the equilibrium configuration; the force balance among the gravity, the pressure gradient force and the centrifugal force are fulfilled: where g R = − ∂ Φ ∂R and g z = − ∂ Φ ∂z .We present a specific method to determine the initial conditions in appendix B.
Figure 2 (a) and (b) show the radial distributions of the initial rotational speed, density, and pressure.We set up a weak vertical magnetic field initially with the R (kpc) ϕ z (kpc) 0.01 < R < 50 −π < ϕ < π −0.5 < z < 0.5 Table 1.Simulation domain following radial dependence: where R b = 0.08 kpc and set to The simulation domain covers the full 2π in azimuth.The radial and vertical domain sizes are 0.01 kpc < R < 50 kpc and −0.5 kpc < z < 0.5 kpc, respectively (Table 1).The outflow condition is employed for the R and z boundaries (e.g., Suzuki & Inutsuka 2014).At the inner R boundary, B z = 0 is set for the vertical magnetic field.
Table 2 presents four cases of numerical simulations conducted for the current paper.Case II is the fiducial case.In Case I the radiative cooling and heating are turned off but the numerical resolution and other setup are the same as those in Case II.By comparison between Cases I and II, we investigate the role of the radiative effects.We reduce (enhance) the net cooling in Cases III (IV) -l by changing the radiative heating rate.In Cases III-l and IV-l lower numerical resolution is employed to save computational time.We also perform half-resolution simulations with the same setup as in Cases I and II, which are labeled as Cases I-l and II-l, in order to directly compare to Cases III-l and IV-l.
The sizes of radial and vertical cells, ∆R and ∆z, are enlarged in proportion to R and z, respectively.In Cases I and II, the minimum cell size is ∆R ≈ ∆z ≈ 1 pc.We run the simulations after the magnetic field strength is sufficiently saturated in a quasi-steady manner.The simulations times are summarized in Table 2. 3. RESULTS

Overview
Figure 4 shows a 3D snapshot in R < 1.5 kpc at t =200.1 Myr in Case II.One can see that the overall magnetic structure is mostly dominated by the toroidal component as a result of the winding due to differential rotation.On the other hand, in the denser central region of R < 300 pc, one can see vertically extended poloidal magnetic fields.We examine in detail the numerical results of the case with the radiative cooling and heating in comparison to those with adiabatic setting below.

Density distribution
We show the time evolution of the number density profiles on the mid-plane (z = 0.0 kpc) in Figure 5   (bottom right) Myr of the two cases.We can see that the radiative cooling and heating substantially affect the density distribution on the mid-plane.In Case I without radiative effects (left of Figure 5), the gas density is considerably decreasing from the initial distribution over time.We can also confirm the formation of clear non-axisymmetric structure.By contrast, the density increases rapidly at an early stage (t = 40 Myr) in Case II with radiative effects (right of Figure 5).In particular, it increases more dramatically in the inner region.
Figure 6 shows the radial distribution of mass density, ρ, normalized by the initial value, ρ 0 , within R < 1 kpc, where ρ is spatially averaged over full azimuthal angle and over time t = 184.4− 203.3 Myr.The gray and black lines show the results of Cases I and II, respectively.The solid (dotted) line depicts the distribution at z = 0.0 (0.1) kpc.Comparison between Cases I and II clearly shows the importance of the radiative cooling and heating.In Case I without radiative effects, the density on the mid-plane is decreasing (ρ/ρ 0 < 1.0) but the density above the mid-plane is increasing (ρ/ρ 0 > 1.0).This is because the gas expands vertically as a consequence of magnetic processes (Section 4.4) and mass is supplied from the mid-plane to the upper regions.On the other hand, in Case II with radiative cooling and heating, the density on the mid-plane is increasing, while the density at z = 0.1 kpc is decreasing.This is a consequence of the downflows from upper regions to the mid-plane.The temperature near the mid-plane drops rapidly by efficient radiative cooling after the simulation starts because of the high density, which causes the decrease in the gas pressure.As a result, the dynamical equilibrium in the vertical direction is no longer satisfied, which induces the downflows.We explain the details in section 3.2.In these simulations, the initially weak vertical magnetic field is amplified by various mechanisms; one of the reliable processes is magnetorotational instability (hereafter MRI ;Velikhov 1959;Chandrasekhar 1961;Balbus & Hawley 1991).MRI occurs in inner-fast differentially rotating systems, and small radial fluctuations of magnetic field lines are amplified through the transport of the angular momentum between inner and outer rings of a disk.The differential rotation also causes winding of poloidal magnetic field lines, resulting in the amplification of B ϕ .In Figure 7, we show the time evolution of the spatially averaged root mean squared (r.m.s.) magnetic field strength (solid lines),

Amplification of the magnetic field
where z , and the averaged ranges are 0.05 kpc < R < 0.8 kpc, −π < ϕ < π and |z| < 0.45 kpc.The gray and black lines are the the results without (Case I) and with (Case II) radiative cooling and heating, respectively.In Case II, the average magnetic field strength (solid black line) is amplified rapidly to 10 µG in a short time within t =10 Myr.In Case I, the magnetic field is amplified relatively slowly, but after t ≳ 200 Myr, the saturated average field strength is almost comparable to that of Case II.The dotted lines in Figure 7 plot the maximum value of the magnetic field among all grid points in the same region for Case I (gray) and Case II (black), respectively.Although the saturation times of the two cases are different, the maximum field strengths in the saturated state reach 200 − 300 µG, which is similar in both cases.Although in the lower resolution Cases I-l and II-l the growth of the magnetic fields is slightly slower than in the high-resolution counterparts of Cases I and II, similar saturated field strengths are achieved in both low and high resolution simulations.
When the radiative cooling and heating are included (Case II), the magnetic field grows more rapidly from the beginning within the initial 10 Myr.In order to understand this situation in more detail, the volume averaged total magnetic field strength , and B ϕ (green), and their time evolutions are compared in Figure 8.It is clear from this figure that the azimuthal component, B ϕ , is most strongly amplified owing to winding by the differential rotational from the early time.In addition, particularly at the beginning, the magnetic field is also compressed by the downflows toward the mid-plane, which contributes to the amplification of B ϕ .The compressible amplification continues until the magnetic pressure-gradient force is balanced by the downward gravity force.
We also show in Figure 9 the radial distribution of the magnetic field averaged over time and the ϕ and z directions The r.m.s.

Time evolution of energy density
We show in Figure 10 the time evolutions of the volume averaged different components of pressure, for Case I (left panel) and Case II (right panel).The averaged region is the same as in Figure 7.Each line denotes magnetic pressure, P B = B 2 /8π, (red solid), thermal pressure, P g = (γ − 1)ρe, (dotted blue), and turbulent pressure, , (green dashed), respectively.Here, the turbulent component is the kinetic energy excluding the initial rotational velocity.
In Case I, the thermal energy density is kept nearly constant with a very slight increase from the beginning, while the turbulent energy density is growing to the same level as the thermal energy density.Although the magnetic energy density is amplified, it is two orders of magnitude smaller than the thermal and turbulent energies, indicating that the magnetic field has a small influence on the global energetics.
On the other hand, in Case II (right panel of Figure 10), it is found that the three different components of the energy densities finally settle down to equipartition each other.The thermal energy density, of which the initial value is about 4 × 10 −11 erg cm −3 , rapidly decreases by an order of magnitude within the initial few Myr.This is due to the rapid drop in temperature caused by the radiative cooling in the dense mid-plane region.As a result, the gas pressure around the mid-plane also drops, which causes the disruption of the equilibrium configuration and excites vertical downflows to the midplane.Thus, the turbulent (kinetic) energy density is enhanced over the thermal energy density at the early stage temporarily.After that, the magnetic energy density is amplified to eventually compensate for the loss of the thermal energy up to t ≲ 10 Myr.The amplification of this magnetic field is contributed by gravitational compression as mentioned above.Winding by differential rotation and MRI also contribute to the amplification of the magnetic fields simultaneously.The figure shows that when the gradient of the total pressure consisting of the magnetic and thermal pressures can support gravity, the turbulent energy including the accretion flow settles down to be roughly comparable to the thermal and magnetic energies.
We should note that, although the equipartition among volume-averaged magnetic, thermal, and turbulent energies are satisfied in a global sense, this is not true in a local sense.Depending on the height from the mid-plane, regions dominated by gas pressure or magnetic pressure are formed, which will be explained later in Section 3.3.

ϕ-averaged poloidal structure
Figure 11 shows the ϕ-averaged density (top), temperature (middle) and inverse of plasma β (bottom) on a R − z plane.The left column in the figure shows the initial distribution (t = 0.0 Myr).We note that the initial temperature ∼ 10 5 K in the mid-plane of the bulge region is higher than realistic interstellar medium in the mid-plane (see Appendix B), and therefore, the initial vertical scale height extends to about 100 pc.The middle and right columns show the results after time evolution of Case I and Case II at t = 200.1 Myr.The comparison between the top-left and top-middle panels clearly show that the gas expands vertically in Case I, as discussed Section 3.1.1.The center panel shows that the temperature at the mid-plane increases to T ≈ 10 6 K in Case I because of magnetic heating (Section 4.4), causing the vertical expansion.
On the other hand, in the Case II, the initial temperature is too high in comparison with the realistic value.Hence, once the simulation started, the temperature rapidly drops around the mid-plane where the density is high.Therefore, the vertical thickness at t = 200.1 Myr (top-right panel) is thinner than that in the initial distribution, indicating that the gas is more concentrated near the mid-plane.In particular, the density at the mid-plane in R < 0.2 kpc exceeds n = 10 3 cm −3 .It can also be seen that the inner dense region is almost flat, but the outer region (R > 0.4 kpc) bulges verti- cally with increasing R. As discussed below, the vertical component of the magnetic pressure-gradient force plays an important role in this structure.The temperature of the high-density gas near the mid-plane is found to be about T = 10 3−4 K.In contrast, low-density and high-temperature gas of T = 10 6−7 K is distributed in the high-latitude halo regions of |z|/R ≳ 1.This indicates the formation of multi-layer structure with a large temperature difference more than a few orders of magnitude from the mid-plane to the halo region like the solar atmosphere that consists of the low-temperature photosphere and chromosphere and the high-temperature corona.
Figure 12 et al. 2018, 2022), in which gas is floating up in the upper region from the mid-plane but falling back downward to the mid-plane.
As we have already mentioned, the effect of cooling and heating plays a significant role in determining these density and temperature distributions.Figure 13 shows n−P diagrams in a range of 0.05 kpc < R < 0.8 kpc and |z| < 0.45 kpc for Case II (with radiative cooling and heating).The left and right panels respectively show the initial condition and the results at t = 200.1 Myr.Note that each point in both panels is illustrated by a translucent circle; a more opaque region with higher concentration corresponds to a larger volume fraction.The three solid straight lines represent the isothermal gas with T = 10 3 , 10 4 , and 6 × 10 7 K on the n − P diagram, respectively.The red dotted line plots the thermal equilibrium curve for the heating rate inside the bulge within R < 1 kpc (Γ = 10 3 Γ 0 and the cut-off temperature of heating T cut = 5 × 10 4 K; see Section 2.3) The left panel of Figure 13 indicates that the initial temperature is higher than the thermal equilibrium temperature (red line) because the initial distribution is not in a thermally equilibrium state although it satisfies the dynamical balance.Since the radiative cooling is more effective for denser gas (equation ( 7)), the gas near the mid-plane cools down more rapidly to the thermal equilibrium temperature in a relatively short time.In fact, the results after time evolution (right panel in Figure 13) show that the high-density gas near the mid-plane is distributed along the thermal equilibrium curve.The pressure also decreases with temperature (equation ( 5)), within a short period of time, especially in the highdensity regions close to the mid-plane where the net cooling rate is large.As a result, the initial equilibrium state is broken so that the gravity cannot be supported by the gas pressure.Downflows are excited, which further increases the density near the mid-plane, resulting in a more cooling-dominated environment.
In contrast, the net cooling rate is smaller in the halo region since the gas density is relatively low.Therefore, the change in gas temperature in this region is expected to be moderate compared with that near the mid-plane.On the contrary, the temperature in the halo region increases from the initial value, as seen in both Figure 11 and Figure 13.A possible reason for the increase in the temperature is that there is the contribution of magnetic heating, which will be discussed in more detail in section 4.4.In summary, the efficiency of radiative cooling is substantially affected by the gas density in the gravitationally stratified structure, resulting in the large variety in the temperature as shown in these figures.
Next, we study the evolution of magnetic field strength by examining the inverse of plasma β in the bottom panels of Figure 11.The initial field strength is weak < 1 µG (Equation 16), and hence, 1/β is far below unity in the simulation domain with < 10 −3 on the mid-plane; the magnetic pressure is almost negligible initially in the simulation region.Although the magnetic field strengths in Cases I (middle-bottom panel) and II (right-bottom panel) are both monotonically amplified with time, these two cases exhibit considerably different distributions of 1/β.In Case I, the gas pressure still remains dominant compared to the magnetic pressure in most of the GB region, although 1/β increases with the amplification of the magnetic field (Section 3.1.2).
Conversely, in Case II, the magnetic pressure dominates the gas pressure in a large portion of the simu-lation domain.As shown in Section 3.1.2,the average saturated field strengths in Cases I and II are similar each other.Therefore, the difference in 1/β between these two cases is primarily due to the difference in the gas pressure.As explained previously, the gas pressure rapidly decreases by the radiative cooling at the early time in Case II.The drop in the temperature is more pronounced near the mid-plane where the initial density is higher.We actually found that the average thermal energy densities in Cases I and II differ over an order of magnitude in Section 3.1.3.Although the temperature in the halo region is kept high in both cases, the density decreases rapidly because of the infall of the gas toward the mid-plane, leading to the drop of the gas pressure there.Therefore, the magnetic field plays a significant role in the dynamics and energetics of the halo gas in Case II.
However, at the mid-plane (z = 0.0 kpc), the density increases locally owing to the compression by the downflows.Since the temperature has reached the thermal equilibrium temperature of the warm neutral medium phase (right panel of Figure 13), the gas pressure increases as the density increases.In fact, the plasma β at the mid-plane is slightly larger than unity ; the vertical distribution of the gas is primarily supported by the gas pressure with the subdominant contribution from the magnetic pressure.On the other hand, the magnetic pressure is more than 10 times higher than the gas pressure in the region slightly above the mid-plane, hereafter we call a "mid-latitude low-β region", which we explain in detail below.It is found that in Case II, while the gas pressure dominates the magnetic pressure on the mid-plane, the magnetic pressure is dominant in the upper layers.Figure 14 compares the vertical profiles of various physical quantities averaged over time between 184.8 Myr and 203.3 Myr at R = 0.15 (top), 0.45 (middle), and 0.75 (bottom) kpc for Case II.The left panels show the vertical profile of number density, n (yellow solid line), and temperature, T (sky-blue dotted line).The middle panels compare the vertical distributions of gas pressure, P g (blue line), magnetic pressure, P B (red line), and total pressure, P t (= P g + P B ; black line).In the right panels, the different components of the pressure are presented against n.In these n − P diagrams, the location close to the mid-plane corresponds to the larger n (right) side.

Physical properties of multi-layer structure
First, we discuss the result at R = 0.45 kpc (middle row of Figure 14).The high-density region with n ≳ 10 cm −3 is formed across the mid-plane with the vertical thickness of |z| ≲ ± a few 10 pc by the initial infall of gas from the upper regions (Section 3.2).The force balance is sustained between the downward gravity and the gas pressure, which dominates the magnetic pressure.However, the gas pressure rapidly decreases with increasing |z| so that in the upper layer of |z| ≳ a few 10 pc, the gravity is mainly supported by the magnetic pressure.From an energetics point of view, the thermal energy reduced by the radiative cooling (Section 3.1.3)is partially compensated by the amplified magnetic energy, which forms the mid-latitude low-β region.Interestingly, similar low-β regions are obtained in a global MHD simulation for disks of spiral galaxies (Wibking & Krumholz 2023) and active galactic nuclei (Kudoh et al. 2020).The further upper halo region is occupied by low-density (n ≲ 10 −2 cm −3 ) and high-temperature (T > 10 6 K) gas.
These features can also be well captured in the n − P diagram (middle-right panel in the Figure 14).The gas pressure almost exactly follows the thermal equilibrium relation in the high density part, n ≳ 0.1 cm −3 .Although the magnetic pressure is subdominant in the Galactic-plane region (n ≳ 10 cm −3 ), it greatly dominates the gas pressure in most of the upper layers.An interesting point is that the plots of P B and P g are almost in parallel with keeping P B /P g (= 1/β) ≈ 10 in the mid-latitude low-β region with 10 −1.5 cm −3 ≲ n ≲ 10 cm −3 .Since P g almost coincides with the thermally equilibrium pressure, P TE , in this region, P B follows the line of 10P TE .In the further lower density part, n ≲ 10 −1.5 cm −3 , which corresponds to the halo region, the gas pressure, P g diviates upward from P TE .This indicates that the halo gas is heated by the dissipation of magnetic energy, which is discussed in Section 4.4 (see also Figure 11).
While the qualitative properties at the inner (R = 0.15 kpc; top row of Figure 14) and outer (R = 0.75 kpc; bottom row) locations are basically the same as those at R = 0.45 kpc, the vertical thickness moderately increases with R because of the weakend gravity.An additional remark for the inner location is that the downward gravity is not fully supported by the gas and magnetic pressures there.Therefore, gas continues to fall down from the halo region (Figure 12).

Thickness of low-β region
We have shown that the magnetically dominated regions are formed in the upper layers above the mid-plane (Section 3.3).We investigate the detailed properties and formation mechanism of these mid-latitude low-β regions by inspecting the vertical distributions of density and magnetic field.Figure 15 shows the vertical profile of density (black solid line) and plasma β (pink dotted line) in the z direction at the same Galactic center distances R = 0.15, 0.45, 0.75 kpc as in Figure 14.
When the vertical gravity is expanded around z = 0, the leading-order term gives where, we note that ω becomes the Keplerian angular velocity for point-source gravity.When the vertical gravity is supported by the gas pressure and vertical motion can be neglected, equation ( 15) gives vertical density distribution.Applying equation ( 20), we obtain where we assumed an isothermal equation of state, P g = ρc 2 s , for a constant sound velocity.The solution of equation ( 21) is where h 1 = c s /ω is the gas pressure scale height.The density, ρ 0 , and the sound velocity, c s , at the mid-plane are evaluated from the numerical data; ω, is calculated from equation 8. Following the same procedure, in the magnetically dominant condition, we derive where Here, |B h | max is the maximum horizontal field strength in the low-β  regions for a given R in the numerical data.Here, ρ ′ 0 is the normalization density that reproduces the density distribution in the upper regions.
The blue and red solid lines in Figure 15 indicate ρ 1 (z) and ρ 2 (z) for each R. The specific values used to find the scale height at each radius are shown in the Table 3. From the figure, the density distribution in the z direction can be well fitted by the sum of the two components, ρ 1 (z) and ρ 2 (z).The density distribution near the midplane (|z/h 1 | < 2) is almost consistent with ρ 1 (z).This is because the density distribution near the mid-plane is mainly supported by the gas pressure so that its thickness is determined by h 1 .In contrast, the magnetic pressure dominates in the upper regions.Thus, the gas is spread out vertically beyond z = h 1 .At R = 0.45 kpc (the middle of Figure 15), the thickness of the gas density in the upper region agrees well with the density distribution of h 2 determined by the magnetic pressure.The same tendency is obtained for the density distribution at R = 0.75 kpc.On the other hand, at R = 0.15 kpc, the density distribution extends beyond the thickness h 2 determined by magnetic pressure.Therefore, the gas located above |z| ≳ 10 pc cannot be supported by magnetic or gas pressure but still continuously falls downward The presence of the mid-latitude low-β regions is a characteristic outcome of the radiative cooling.Our treatment for the radiative cooling and heating is an approximated one (Section 2.3).Thus, we investigate the impact of the magnitude of the net cooling rate on the density distribution with 3D MHD simulations by changing the radiative heating rate.In Case III-l (IV-l) we increase (decrease) the radiative heating rate, Γ, to 10 (1/10) of the original setting in Case II.Since the increase (decrease) of the heating corresponds to the reduction (enhancement) of the net cooling, we call Case III-l (IV-l) the case with reduced (enhanced) cooling.
As described in Section 2.4, in Cases III-l and IV-l, we employ half the numerical resolution compared to the fiducial models of Cases I and II.The low resolution simulations are performed for Cases I and II4 .We compare these four cases (Cases I-l -IV-l) with the same resolution.
Figure 7 indicates the volume averaged magnetic field strengths in 0.05 kpc < R < 0.8 kpc and |z| < 0.45 kpc are very similar in Cases II-l, III-l and IV-l. Figure 16 compares the vertical distributions of Cases I-l, II-l, III-l and IV-l at R = 0.45 kpc.The qualitative vertical distributions of Cases III-l and IV-l remain unchanged even under the different net cooling rates; low-temperature and high-density gas around the midplane is sandwiched by mid-latitude low-β regions (see the lower right panel of Figure 16).However, the difference in the net cooling significantly influences in the upper region away from the mid-plane.In the reduced net cooling case (Case III-l), it is found that the density (top left), magnetic field (top right), and temperature (left bottom) structures exhibit intermediate characteristics between Case I-l and Case II-l.In the enhanced net cooling case (Case IV-l), the structure is similar to Case II-l.
In the mid-plane where the gas density is high, the cooling rate is significantly dominant over the heating rate, even though the heating rate is increased by a factor of 10 as done in Case III-l.Therefore, the net cooling rate is almost unaffected there and the system reaches a thermal equilibrium state within a short time.On the other hand, in the upper region with lower density, the cooling rate decreases, and thus, the difference in heating rate has a significant impact on the net cooling rate.
As a result, the difference in the heating rate affects the cooling time compared to the dynamical time, leading to the observed difference in the vertical distribution in the upper regions.
These results indicate that the height at which magnetic pressure dominates (β < 1.0) is controlled by the adopted net cooling rate.This is because the thermal equilibrium temperature, which depends on the net cooling rate, determines the gas pressure.In our simulations, we are using the simplified treatment for the radiative cooling and heating.In particular, we have assumed the heating rate Γ 0 for a typical interstellar cloud near the solar system in equation ( 13).In reality, photoelectric heating and heating by X-rays and cosmicrays depend on the surrounding environment in different ways.Also, we are not solving the ionization states of atoms that are crucially important in determining radiative cooling.More realistic simulations with these effects should be carried out in our future works.

Distribution of magnetic field strength
In our MHD simulations, the magnetic field strength is amplified to > 10 µG with the maximum strength of a few to several 100 µG, which is consistent with the observationally suggested value, 10 − 1000 µG in the GB region (Morris & Serabyn 1996;Pillai et al. 2015;Crocker et al. 2010).In this subsection, we further examine the distribution of the obtained magnetic field from a statistical viewpoint, particular focusing on the dependence on density.
Figure 17 shows the histograms as a function of the magnetic field strength in 0.05 kpc < R < 0.3 kpc, which corresponds to the scale of the CMZ.The volume fraction (the top panel of Figure 17) is dominated by low-density gas above the mid-plane.Since the gravity is stronger at small R, the gas that was cooled by radiative cooling cannot maintain the dynamical equilibrium.Therefore, low-density gas falls down toward the mid-plane and forms a compressed dense region.The two peaks in the volume fraction diagram correspond to the compressed high-density gas near and on the mid-plane and low-density gas that falls from the upper layers.From the volume fraction histogram, the spatially averaged magnetic field strength is estimated to be about 10 µG.
The bottom panel of Figure 17   The mass spectrum can focus on denser regions; therefore, the mass weighted average magnetic field strength is higher than the volume weighted field strength, with the peak close to 100 µG.In particular, the mass fraction that exceeds 100 µG is about 11% (see in Figure 18).In the inner part of the GB region, there are dense clouds with strong magnetic fields where the field strength locally exceeds 1 mG.For one specific example, Pillai et al. (2015) estimated to be magnetic field strengths of 1-2 mG inside the "Brick" region, a highdensity dark nebula in the CMZ.Such strong-field regions tend to be correlated with high-density molecular clouds.The total mass of CMZ at R < 300 pc is M CMZ = 2 − 6 × 10 7 M ⊙ (Dahmen et al. 1998;Ferrière et al. 2007).The average molecular gas density in the CMZ is ρ = 10 −20 g cm −3 (Güsten & Henkel 1983;Bally et al. 1987;Mills et al. 2018).Therefore, the volume of molecular gas in the Galactic center region is estimated to be V = M CMZ /ρ ∼ 7×10 60 cm 3 , which is about 1% of the total volume in R < 300pc and |z| < 100 pc.On the other hand, the molecular gas is likely to be the major component in mass since it accounts for more than 90% of the total mass of the entire CMZ scale.The current numerical setting cannot handle such high-density and low-temperature clouds because of the the cut-off temperature for the radiative cooling (equation ( 7)) and the insufficient numerical resolution.In our future study, we plan to investigate strong-field concentrations in highdensity regions by improving the numerical treatment.

Cosmic-Rays and γ-Ray Excess in the GB region
In this paper, we show that the magnetic field supports the density structure in the z direction in the upper region of the mid-plane.On the other hand, it is suggested that cosmic-ray pressure may contribute to the same extent (Boulares & Cox 1990).It has also been reported that cosmic-rays can be a source of heating of dilute ionized gas and drive outflow (Shimoda & Inutsuka 2022).The study of dilute gas at higher altitudes (z > 1 − 10 kpc) should incorporate the effects of cosmic-rays and supernova explosions, which are not considered in this study.
The cosmic-rays around the GB region are also related to the observed ∼ 100 TeV γ-ray excess (HESS Collaboration et al. 2016).The cosmic-ray protons with an energy of ∼ PeV can be responsible for the ∼ 100 TeV γray emissions although their origin and acceleration process are still under debated (e.g., Shimoda et al. 2022).In the most accepted scenario of cosmic-ray acceleration, the diffusive shock acceleration at supernova remnant shocks, the maximum energy of accelerated particles is determined by the magnetic filed strength at the shock upstream (e.g., Lagage & Cesarsky 1983a,b;Bell 2004).The required field strength is > 100 µG for accelerating ∼PeV cosmic-rays at the supernova remnant shock.If the volume fraction of gas with such a strong field was large, the PeV cosmic-ray acceleration would occur everywhere in the GB region.Thus, further investigation of the field strength distribution is a necessary step to elucidate the origin of the high energy γ-rays and cosmic-rays.

Magnetic heating
In Figures 13 and 14, we explained that the lowdensity gas in the halo region is significantly heated from the thermal equilibrium temperature to a hightemperature state.We here discuss the role of the dissipation of magnetic fields in heating the halo gas.
The magnetic energy can be transformed to kinetic and thermal energy via magnetic reconnection processes that happen in the geometrically thin current sheets where we cannot neglect physical resistivity.In realistic situations, it is regarded that magnetic turbulence induces small-scale reconnections (Lazarian & Vishniac 1999;Lazarian et al. 2020), which finally converts the magnetic energy to the thermal energy of the surrounding gas.In general, such reconnections of turbulent magnetic fields take place at very fine scales, which requires extremely small sizes of numerical cells.It is not realistic to perform such high-resolution simulations in our global treatment.Instead, utilizing the energy conserved scheme we are employing (Section 2.1), we estimate the heating rate that is expected from the dissipation of magnetic energy (Suzuki & Inutsuka 2006;Matsumoto & Suzuki 2014).
In our simulations dissipative phenomena are treated by the MHD Riemann solver; the physical heating by viscosity within the sub-grid scale is handled by this shock capturing scheme.Although this method cannot describe and directly capture the details of the physical dissipation scale, we believe that it can give a reasonable estimate on the global magnetic heating rate.
The magnetic energy conservation law is expressed as where the term on the right-hand side denotes Joule heating rate using the electric current J and the electric conductivity σ c .We note that, since σ c = ∞ is assumed in the ideal MHD, the right-hand side is zero in the ideal MHD limit.However, what we are solving in our sim- t=200.1 Myr (0.05 < R < 0.8, |z| < 0.45 kpc) Figure 19.
Volume average magnetic heating rate for number density.
ulations is not equation (25) but equation (3) for the conservation of the total energy.Therefore, the sum of the three terms on the left-hand-side of equation ( 25) is not zero numerically, which gives a finite value of the Joule heating rate.This non-zero Joule heating at each numerical cell is counted in the change of the thermal energy.This increase in the thermal energy due to the numerical dissipation can be interpreted as the energy released at MHD shocks and by small-scale reconnections as a consequence of turbulent cascade that occurs at a sub-grid scale.In the following discussion, the Joule heating rate, J 2 /σ, derived numerically from equation ( 25) is referred to as the "magnetic heating rate".
Figure 19 shows magnetic heating rate (red solid line) averaged in the range of 0.05 kpc < R < 0.8 kpc and |z| < 0.45 kpc at t =200.1 Myr, in comparison to the radiative cooling rates for different temperatures (solid bluish lines).Although the magnetic heating is largely dominated by the radiative cooling in the high-density regions with n ≳ 1cm −3 , it is not negligible in lowerdensity regions.In particular, for n < 10 −2 cm −3 , the magnetic heating rate exceeds the radiative cooling rate, causing the deviation of the temperature in the halo region from the equilibrium value shown in Figures 13  and 14.
Although we are taking account for radiative cooling and heating in the energy equation, we do not consider the effect of thermal conduction (diffusion).The timescale, τ i , of the thermal conduction for fully ionized gas in the GB region can be estimated as follows: where κ i = 1.24 × 10 −6 T5/2 [erg K −1 cm −1 s −1 ] is the thermal conductivity (Parker 1953;Asai et al. 2007).This estimate indicates that thermal conduction is more effective in lower-density and higher-temperature envi-ronments, where magnetic heating is also actively working.If we substitute typical temperature, T = 10 6 − 10 7 K, and density, n < 10 −2 cm −3 , into equation (26), we obtain τ i ranging from 10 −4 Myr to 0.1 Myr.Since this is much shorter than the rotation time ∼ a few 10 Myr in the GB region, it is expected that the temperature of the hot gas should have a flatter temperature profile by thermal conduction than that shown in Figure 11.
On the other hand, magnetic heating and radiative heating are weaker against radiative cooling in the highdensity regions near the mid-plane.Thus, the averaged gas temperature is low ≲ 10 4 K there.In such low-temperature gas, neutral atoms contribute to the thermal conduction.The corresponding conductive timescale is 1953;Koyama & Inutsuka 2000).Therefore, the effect of thermal conduction is negligible on the scale of L ≳ 10 − 100 pc near the mid-plane.

Application to Milky Way -the stellar bar
Our simulations are focusing on the GB region of general spiral galaxies with axisymmetric bulge structure (e.g., M31, NGC 278, NGC 628 5 , NGC 772 and NGC 4030).If we are to compare to the Milky way, we have to take into account the effect of the bar potential.This is because the stellar bar structure with a size of about 3 kpc in our galaxy is widely known from infrared observations (Matsumoto et al. 1982;Nakada et al. 1991;Weinberg 1992;Dwek et al. 1995).How the gravitational field affects the distribution of the interstellar gas and its motion inside this bar structure has been discussed for more than 30 years.Assuming the inner bar potential based on the basic theory of Binney et al. (1991), several intersecting orbits are formed, and the orbits near the Galactic center, called X1 and X2 orbits, reproduce the elongated orbits of the gas in and around the CMZ and the gas accretion flow along the dust lanes (Athanassoula 1992).The orbits and structures of the interstellar gas in the bar potential have been examined in detail by N-body simulation (Rodriguez-Fernandez & Combes 2008).In recent years, high-resolution hydrodynamic simulations provide quantitative results, and the impact on the bar potential on the star formation history is being discussed (Sormani et al. 2018;Armillotta et al. 2019;Baba & Kawata 2020;Tress et al. 2020).Indeed, it is identified that the contribution of the stellar bar forms a stream along dust lanes with the time-averaged inflow rate of a few M ⊙ yr −1 (Sormani & Barnes 2019).
It is also discussed how the bar potential contributes to the gas distribution and motion via the magnetic field in the Galactic center region (Kim & Stone 2012).Recently, Moon et al. (2023) performed 3D MHD numerical simulations under the magnetized inflowing gas along a dust lane into the Galactic center region, which is expected from the stellar bar.They reported that the amplified magnetic field plays a substantial role in the suppression of star formation in the nuclear ring region and the excitation of inflows from the ring to the nuclear disk.In our next step, we are also planning to include the effect of the stellar bar in global MHD simulations.

CONCLUSION
We examined the impact of magnetic fields on gas dynamics in the GB region using 3D-MHD simulations with radiative cooling and heating.Although the thermal properties are greatly affected by the radiative cooling, the saturated field strengths are not so different in the cases with and without the radiative effects; the magnetic fields are amplified from the initial value < µG to the volume integrated average of ≈ 10 µG with several µG in localized field concentrations.We also found rough equipartition between the volume-averaged magnetic and thermal energies, whereas the local plasma β value varies with height.This result emphasizes the substantial role of magnetic fields in controlling dynamical and thermal properties in the GB region.
An important finding is that, by the effect of the radiative cooling in combination of magnetic fields, a vertically multi-layered structure is created.The cool dense gas that is mainly supported by gas pressure is distributed in a thin layer on and near the mid-plane.Above that, low-β plasma occupies the vertical range of z = 50 − 150 pc.The magnetic pressure, which largely dominates the gas pressure in the mid-latitude low-β regions, mainly maintains this vertical thickness comparable to the observed scale height of neutral atomic gas.The further upper regions are occupied by low-density and high-temperature halo gas.
In the cool dense mid-plane and the mid-latitude lowβ regions, the temperature is determined by the thermal equilibrium between the radiative cooling and heating.In contrast, the temperature is largely deviated from the equilibrium value in the Galactic halo.This indicates that the magnetic heating may be important in such low-density gas.Sutherland & Dopita (1993).We used the fitting function of the fourth order equation created from the table data from Sutherland & Dopita (1993) as Λ h at the high temperature side of T > 10 4 K:

B. INITIAL DISTRIBUTION
We set up the initial distributions, following Suzuki et al. (2015).First, we determine the temperature distribution (the middle left of Figure 11).We give T = 1.25 × 10 5 K at the mid-plane in the bulge, T = 5 × 10 3 K at the mid-plane in the disk, T = 2.0 × 10 6 K at z = 0.5 kpc in the bulge, and T = 2.0 × 10 5 K z = 0.5 kpc in the disk.We smoothly interpolate temperature in the boundary layers between these four regions.
Next, the density at the mid-plane (z = 0.0 kpc) is assumed to be 0.05 times the stellar mass distribution that is determined by Φ * (equation 9 and Section 2.2).We integrate equation 15 to determine P g (R, z), using g z = − ∂ Φ ∂z .As temperature is already fixed, the density distribution ρ(R, z) is derived from P g .Finally, by integrating equation ( 14) along the radial direction, we determine the initial rotation frequency.The top left panel shows that density around the mid-plane is slightly higher in the high-resolution case than in the low-resolution case.The magnetic field strength (top right) is also slightly more amplified near the mid-plane plane in the high-resolution case, compared with the low-resolution case.On the other hand, temperature (bottom left) in the mid-plane region does not depend on numerical resolution.This is because the thermal equilibrium temperature determined by radiative cooling is not affected by numerical resolution.However, temperature in the halo region is higher for higher resolution owing to the resolution dependence of magnetic heating.

Figure 2 .
Figure 2. Initial distribution for Galactic distance R. (a) Velocity and strength of differential rotation.(b) Density, gas pressure and magnetic pressure.(c) G(R).
of 'Brick'.Sormani et al. (2018) introduced a functional form for G(R) to mimic the empirically derived profile (Figure2(c)), which we adopt in our simulations.The heating rate Γ(R) is treated in our simulations as proportional to G(R).

Figure 4 .
Figure 4. 3D overview inner the GB region at t =200.1 Myr in Case II.Colors represent number density.Red lines represent magnetic field lines.

Figure 5 .Figure 6 .
Figure 5.Time evolution of the number density profile on the mid-plane (z = 0.0 kpc).The left four panels show the results for Case I (without radiative cooling and heating) and the right four panels show the results for Case II (with radiative cooling and heating).The snapshots are taken at t = 0.0, 40, 100, and 200 Myr.The arrows represent velocity vectors in the mid-plane.

Figure 7 .
Figure 7. Time evolution of the volume-averaged magnetic field strengths of various cases (top) and the maximum field strengths in the same region for two cases (bottom).

Figure 8 .
Figure 8.Time evolution of radial (cyan dot-dashed), azimuthal (dark-green dashed), and vertical (pale red dotted) components of the magnetic field averaged over the R, ϕ and z directions.The averaged range is the same as in Figure 7.

Figure 9 .
Figure 9. Radial profile of magnetic field.The lines are the same as in Figure 8, except for gray dotted line (initial condition).

Figure 10 .
Figure 10.Time evolution of thermal (blue dotted), magnetic (red solid), and turbulent (green dashed) energy densities averaged over the R, ϕ and z directions.

Figure 11 .
Figure 11.Number density n (top), temperature T (middle) and an inverse of plasma β (bottom) averaged over the ϕ direction in a R − z plane.The left panels present the initial profile.The middle and right panels respectively present the profile at t = 200.1 Myr without (Case I) and with (Case II) cooling and heating.
shows velocity field (arrows) with density (colors) on a R − z diagram.It is shown from the figure that gas flows upward in the halo region in both Cases I and II.The averaged mass outflow rate in 0.05 kpc < R < 0.8 kpc is ∼ several times 10 −2 M ⊙ yr −1 in Case I and several times 10 −4 M ⊙ yr −1 in Case II.Additionally, one can clearly see failed winds (e.g., Takasao

Figure 12 .
Figure 12.Velocity fields (purple arrows) with density distribution (colors) in Case I (top left) and II (top right).Top two panels cover the same area as in Figure 11.The bottom panel zooms in 0.05 < R [kpc] < 0.2 and |z| < 0.1 kpc of Case II.

Figure 14 .
Figure 14.Vertical distributions of various physical quantities at R = 0.15 kpc (top), 0.45 kpc (middle), and 0.75 kpc (bottom) of Case II.The left panels present number density (solid orange) and temperature (light blue dashed).The middle panels compare total pressure (black solid), gas pressure (blue solid), and magnetic pressure (red solid).The right panels show n − P diagrams; the same line types are used for pressure; the thermal equilibrium curved for PTE (gray dashed) and 10 PTE (gray dotted) are also plotted (see text for the detail).

Figure 15 .
Figure 15.Density distribution in Figure 14 zoomed near the mid-plane.The black solid line is the density of the numerical simulation of Case II.The red and green solid lines are the density profiles, supported by gas pressure, ρ1 (equation 22), and magnetic pressure, ρ2 (equation 23), respectively.The purple dashed line represents the plasma β profile.The top axis is normalized by gas scale height, h1.

Figure 16 .
Figure 16.Comparison of vertical profiles for different heating rates.Each panel represents vertical distributions at R = 0.45kpc of density (top left), magnetic field strength (top right), temperature (bottom left), and plasma β values (bottom right).The lines correspond to the results of Case I-l (adiabatic case; cyan solid), Case II-l (fiducial cooling and heating; pale red dashed), Case III-l (reduced net cooling; green dotted), and Case IV-l (enhanced net cooling; yellow dot-dashed), respectively.

Figure 17 .
Figure 17.Histograms of the volume fraction (top) and mass (bottom) against magnetic field strength.Each bin is spacing via ∆ ln B = 0.02.

Figure 18 .
Figure 18.Cumulative mass for magnetic field strength.The integration is taken from the maximum |B|max = 200 µG to smaller |B|.The horizontal solid and dashed lines indicate fraction of 0.1 and 0.5, respectively.
K. K. thanks Y. Fukui, R. Enokiya and A. Tanikawa for fruitful discussions.This work was supported by JSPS KAKENHI Grant Nos.19J11052 (K.K.), 17H01105, 21H00033, 22H01263 (T.K. Suzuki) and 18H05436, 18H05437 (S.Inutsuka), and by Program for Promoting Research on the Supercomputer Fugaku by the RIKEN Centre for Computational Science (Toward a unified view of the universe: from large-scale structures to planets, grant 20351188 -PI J. Makino) from the MEXT of Japan.This research used computational resources of Oakforest-PACS provided by Multidisciplinary Cooperative Research Program in Center for Computational Sciences, University of Tsukuba, Yukawa-21 at YITP in Kyoto University and Cray XC50 at Center for Computational Astrophysics, National Astronomical Observatory of Japan.
FUNCTION FOR RADIATIVE COOLING In the Figure 1(a), the black points show the data table of
Figure20compares the vertical distributions of various physical quantities taken from Cases I, II, I-l, and II-l.The top left panel shows that density around the mid-plane is slightly higher in the high-resolution case than in the low-resolution case.The magnetic field strength (top right) is also slightly more amplified near the mid-plane plane in the high-resolution case, compared with the low-resolution case.On the other hand, temperature (bottom left) in the mid-plane region does not depend on numerical resolution.This is because the thermal equilibrium temperature determined by radiative cooling is not affected by numerical resolution.However, temperature in the halo region is higher for higher resolution owing to the resolution dependence of magnetic heating.

Table 2 .
Case NR N ϕ Nz Radiative cooling and heating Simulation time Summary of simulation cases.

Table 3 .
The value used to estimate scale height at each radius.