Formation of Galactic Prominence in the Galactic Central Region

We adopted the in situ formation model of prominences reported by Kaneko & Yokoyama (2015). Figure 1 schematically shows the formation mechanism of dense, cold filaments. We consider static magnetic arches at the initial state (Figure 1(a)). Subsequently, we impose converging and shearing motion at the footpoints of the arch to induce the rising motion of the arcade. Note that both the magnetic buoyancy around the loop top of the arch and the magnetic pressure gradient enhanced by the footpoint motion contribute to the rising motion of the arch. Current sheets can be formed inside the expanding magnetic arcade. Magnetic reconnection taking place in the current sheet forms rising dense plasmoids (Figure 1(b)). Since the plasmoid is confined by the helical magnetic fields, the warm gas accumulates around the bottom of the closed poloidal magnetic field lines (red arrows), where the cooling instability forms dense, cold filaments supported by magnetic tension force (Figure 1(c)). In this paper, since we neglect the variation of physical quantities perpendicular to the poloidal plane, the dense filaments are not arches but straight filaments. In three dimensions, the filaments can bend and form an arch, along which the molecular gas slides down.

Zoom In Zoom Out Reset image size

The basic equations of resistive magnetohydrodynamics (MHD) are as follows:

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}+{\rm{\nabla }}\cdot (\rho {\boldsymbol{v}})=0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial (\rho {\boldsymbol{v}})}{\partial t}+{\rm{\nabla }}\cdot (\rho {\boldsymbol{v}}{\boldsymbol{v}})=-{\rm{\nabla }}P-\rho {\boldsymbol{g}}+\displaystyle \frac{({\rm{\nabla }}\times {\boldsymbol{B}})\times {\boldsymbol{B}}}{4\pi },\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}={\rm{\nabla }}\times ({\boldsymbol{v}}\times {\boldsymbol{B}}-\eta {\rm{\nabla }}\times {\boldsymbol{B}}),\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}\displaystyle \frac{\partial E}{\partial t} & + & {\rm{\nabla }}\cdot \left[\left(E+P+\displaystyle \frac{{B}^{2}}{8\pi }\right){\boldsymbol{v}}\right.\\ & - & \left.\displaystyle \frac{{\boldsymbol{B}}({\boldsymbol{B}}\cdot {\boldsymbol{v}})-\eta ({\rm{\nabla }}\times {\boldsymbol{B}})\times {\boldsymbol{B}}}{4\pi }\right]=\rho {\boldsymbol{v}}\cdot {\boldsymbol{g}}-\rho { \mathcal L },\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&E=\displaystyle \frac{P}{\gamma -1}+\displaystyle \frac{\rho {v}^{2}}{2}+\displaystyle \frac{{B}^{2}}{8\pi }.\end{eqnarray} \tag{ 5 }$$

where the anomalous resistivity η is assumed to depend on the current density J (see Kaneko & Yokoyama 2015) as

$$\begin{eqnarray}\eta =\left\{\begin{array}{ll}0 & J\lt {J}_{c},\\ {\eta }_{0}{(J/{J}_{c}-1)}^{2} & J\geqslant {J}_{c}.\end{array}\right.\end{eqnarray} \tag{ 6 }$$

Here ${\eta }_{0}$ = 3 × 10²³ cm² s⁻¹ and J_c = 4.0 × 10⁻¹⁷ dyn${}^{1/2}$ cm s⁻¹. We adopt Cartesian coordinates $(x,y,z)$ and assume translational symmetry with respect to y. The gravitational acceleration is assumed to be ${\boldsymbol{g}}=(0,0,{g}_{z})$ and g_z = 3 × 10⁻⁹ cm s⁻². The cooling function $\rho { \mathcal L }$ includes cooling and heating terms, and other symbols have their normal meaning. In the interstellar medium, cold neutral medium (CNM), with temperature 10–100 K, and warm medium (WM), with temperature 10⁴ K, coexist by the balance of cooling and heating (Field et al. 1969; Wolfire et al. 1995, 2003). We adopted the cooling/heating function summarized by Inoue et al. (2006):

$$\begin{eqnarray}&&\rho { \mathcal L }=(-{\rm{\Gamma }}+n{\rm{\Lambda }})n,\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{\rm{\Gamma }}=2\times {10}^{-26}\ \mathrm{erg}\ {{\rm{s}}}^{-1},\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}{\rm{\Lambda }} & = & 7.3\times {10}^{-21}\ \exp \left(\displaystyle \frac{-118400}{T+1500}\right)\\ & & +7.9\times {10}^{-27}\ \exp (-92/T)\ \mathrm{erg}\ {\mathrm{cm}}^{3}\ {{\rm{s}}}^{-1}.\end{eqnarray} \tag{ 9 }$$

Here n is number density. Figure 2 shows a thermal equilibrium curve for the interstellar medium determined by the cooling/heating function (7)–(9). Two stable branches, WM (red branch) and CNM (blue branch), exist. In this paper we assume that cooling/heating takes place at the temperature range 400 ${\rm{K}}\lt T\lt {\rm{20,000}}$ K; otherwise, we set $\rho { \mathcal L }=0$ in order for the dense, cold region to be resolved numerically with a moderate number of grid points. For simplicity, we neglect thermal conduction.

Zoom In Zoom Out Reset image size

The basic equations are solved by applying the HLLD scheme (Miyoshi & Kusano 2005). Second-order accuracy in space is preserved by linearly interpolating the values at the cell interface and restricting them using the minmod limiter. The solenoidal condition ${\rm{\nabla }}\cdot {\boldsymbol{B}}$ = 0 is satisfied by applying the generalized Lagrange multiplier (GLM) scheme proposed by Dedner et al. (2002). The cooling/heating term is included by the time-implicit method. This simulation code has been applied to MHD simulations of the interaction of jets with interstellar gas (Asahina et al. 2014).

The size of the simulation box is −200 pc $\lt \,x\,\lt$ 200 pc, and 5 pc $\lt \,z\,\lt$ 2400 pc. The grid size is uniform with Δx = 1 pc and Δz = 0.5 pc.

We assume a gravitationally stratified layer in hydrostatic equilibrium. The initial magnetic field is assumed to be force-free and given by

$$\begin{eqnarray}&&{B}_{x}=-\left(\displaystyle \frac{2{L}_{a}}{\pi {H}_{m}}\right){B}_{a}\cos \left(\displaystyle \frac{\pi }{2{L}_{a}}x\right)\exp \left(-\displaystyle \frac{z}{{H}_{m}}\right),\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{B}_{y}=\sqrt{1-{\left(\displaystyle \frac{2{L}_{a}}{\pi {H}_{m}}\right)}^{2}}{B}_{a}\cos \left(\displaystyle \frac{\pi }{2{L}_{a}}x\right)\exp \left(-\displaystyle \frac{z}{{H}_{m}}\right),\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{B}_{z}={B}_{a}\sin \left(\displaystyle \frac{\pi }{2{L}_{a}}x\right)\exp \left(-\displaystyle \frac{z}{{H}_{m}}\right).\end{eqnarray} \tag{ 12 }$$

Initial distributions of density, temperature, and magnetic field are shown in Figure 3. We set ρ(z = 5 pc) = 2.8 × 10⁻²⁴ g cm⁻³ and P(z = 5 pc) = 1.56 × 10⁻¹² dyn cm⁻² (T(z = 5 pc) = 6700 K). Here H_m = 200 pc and L_a = 200 pc denote magnetic scale height and half width of the magnetic arch, respectively. We assume a strong magnetic field to sustain the dense filaments. The plasma β (=P_gas/P_mag) is assumed to be β = 0.2 at the bottom of the simulation area (B_a = 1.54 × 10⁻⁵ G). The warm gas with temperature 6700–20,000 K is assumed to satisfy mechanical equilibrium and thermal equilibrium before converging and shearing motions are imposed. In the region where $z\,\gt$ 900 pc, we assume a hot corona with T_corona = 2 × 10⁶ K.

Zoom In Zoom Out Reset image size

Boundary conditions in the x-direction are symmetric for ρ, P, B_z, v_z and antisymmetric for B_x, B_y, v_x, v_y. Absorption boundary conditions are applied at the upper boundary. At the lower boundary, density and pressure are fixed to the initial value. Converging (v_x) and shearing (v_y) motions are imposed at footpoints of the magnetic arch as follows:

$$\begin{eqnarray}{v}_{x}={v}_{y}=\left\{\begin{array}{ll}-{v}_{0}\displaystyle \frac{x}{{L}_{a}/4} & 0\leqslant x\lt {L}_{a}/4,\\ -{v}_{0}\displaystyle \frac{{L}_{a}-x}{3{L}_{a}/4} & {L}_{a}/4\leqslant x\leqslant {L}_{a},\end{array}\right.\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{v}_{z}=0.\end{eqnarray} \tag{ 14 }$$

Here v₀ is 4 km s⁻¹.

Figure 4 shows the density and temperature distribution in the x–z plane at 30, 40, and 120 Myr. The arcade field squeezed by the converging motion and stretched by the shearing motion forms a current sheet around x = 0 inside the arcade. Magnetic reconnection taking place in the current sheet creates rising flux ropes inside which warm gas is confined. In this simulation, the flux rope is lifted up to ∼250 pc. Since the warm interstellar gas slides down along the closed magnetic field lines, the density around the bottom of the flux rope exceeds the threshold for the onset of the cooling instability. Thermal instability taking place in the dense region forms cold, dense filaments. The density of the filament becomes 10–100 times that of the initial state. The dense filament is lifted up to 100–200 pc and sustained for 100 Myr. The length and thickness of the filament at 120 Myr are ∼60 pc in the vertical direction and 6–8 pc in the x-direction, respectively.

Zoom In Zoom Out Reset image size

Figure 5 enlarges the region where dense filaments are formed. Arrows show velocity vectors. The warm gas infalls along the magnetic field lines. The speed of the warm gas accumulating toward the filament is about 2 km s⁻¹. Figure 6 shows the distribution of B_y at t = 120 Myr. Since B_y is amplified by shearing motion, B_y becomes strong inside the magnetic flux rope. The flux rope corresponds to the blue-purple region in Figure 6, where ${B}_{y}\,\gt$ 5.5 μG. Magnetic field lines inside the flux rope have a helical shape, as schematically shown in Figure 1(b).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Red dots in Figure 7 show the pressure and number density at each grid point. At the initial state (Figure 7(a)), warm gas is located on the thermal equilibrium curve (solid line). The dashed line shows the isothermal line at T = 400 K, and the dot-dashed line shows that at T = 20,000 K. Figure 7(b) plots pressure P and number density n at t = 120 Myr. Since the warm gas infalling along the closed magnetic field lines compresses the gas around the bottom of the flux rope (x = 0 pc), the number density exceeds the threshold for the cooling instability and deviates from the warm gas branch. The dense, cool gas is transformed to the CNM, where temperature is less than 400 K. Since we cut off cooling when $T\lt 400$ K to numerically resolve the filament, the CNM is located at the isothermal line at T = 400 K.

Zoom In Zoom Out Reset image size

Figure 8 shows the time evolution of density, pressure, and temperature at x = 0 pc and z = 105 pc. The initially warm gas ($T\sim 10$⁴ K) is compressed during 30–50 Myr. Consequently, density and pressure increase and temperature slightly decreases. At t = 50 Myr, the gas is located in the thermally unstable region in the P–n diagram. Figure 8 shows that pressure and temperature decrease owing to the cooling instability and the density continues to increase. Finally, a cool (T = 400 K), dense filament is formed.

Zoom In Zoom Out Reset image size

Figure 9 shows the vertical distribution of temperature at x = 0 pc at 120 Myr. The dense filament locates at z = 100–160 pc. Meanwhile, a high-temperature (T = 35,000–70,000 K) region appears at z = 330–450 pc. The hot region is formed by shock waves formed by the upward flow. When the temperature of the shocked gas exceeds 20,000 K, the region stays in a high-temperature state.

Zoom In Zoom Out Reset image size

Figure 10 shows the column density of the dense gas ($T\lt 500$ K), the unstable neutral medium (UNM) in the temperature range 500 K $\lt \,T\,\lt$ 5000 K, and the warm gas in the temperature range 5000 K $\lt \,T\,\lt$ 7000 K, at 120 Myr for the line-of-sight direction parallel to the x-axis. The column density of the dense filament is ∼1.0 × 10⁻³ g cm⁻² (column number density = 6 × 10²⁰ cm⁻²) around the bottom of the dense filament. The column density of the UNM is largest at 1.0 × 10⁻⁴ g cm⁻² (column number density = 6 × 10¹⁹ cm⁻²). The column density of the warm gas is 3 × 10⁻⁴ g cm⁻² above the top of the dense filament. The dense gas can be observed by CO emission, and the UNM and the warm gas can be observed by 21 cm line emission by neutral hydrogen. Torii et al. (2010b) compared the distribution of CO emission and H i emission and found that they coincide in loop 1. In our simulation, the distribution of the UNM coincides with the dense filament, but the warm gas above the dense filament has comparable column density to the dense filament.

Zoom In Zoom Out Reset image size

The column density estimated from the CO emission is 3 × 10²¹ cm⁻² around the top of loop 1, and the H i column number density of loop 1 obtained from the 21 cm line emission is 6 × 10²⁰ cm⁻² (Torii et al. 2010b). These column densities are five times larger than that of our simulation and will be discussed in the next subsection.

In this section, we would like to estimate the mass levitated in the cold, dense filament and compare the numerical results with observations. NANTEN observations indicate that the projected lengths of two loops (loop 1 and loop 2) are ∼500 and ∼300 pc, respectively. They are lifted up to ∼220 and ∼300 pc, and the total mass of each loop was estimated to be 8 × ${10}^{4}$ ${M}_{\odot }$ (Fukui et al. 2006). If we take the average value of projected lengths of two loops and assume that the shape of the loops is a half circle, the average length of the loop is 400 pc × π/2 ∼ 600 pc. In our simulation results at 120 Myr, the dense filament ($T\lt 500$ K) is lifted up to 100–200 pc. The mass of the filament can be estimated to be 112 ${M}_{\odot }$ times the length of the loop in the y-direction in units of pc. In other words, the total mass of the dense filament is 7 × 10⁴ ${M}_{\odot }$ if we assume that the length of the filament in the y-direction is 600 pc.

The total mass of the dense filament in the simulation is consistent with the lower limit of the mass estimated from observations. However, Torii et al. (2010b) estimated by detailed analysis of the NANTEN ¹²CO and ¹³CO(J = 1–0) observations that the total mass of these two loops is 1.4 × 10⁶ ${M}_{\odot }$ and 1.9 × 10⁶ ${M}_{\odot }$, respectively. These values are an order of magnitude larger than those of our simulation. The smaller total mass of the molecular loop in our simulation can be due to the lower density of the WM, ${\rho }_{\mathrm{WM}}$ = 2.8 × 10⁻²⁴ g cm⁻³, that we assumed at the bottom of our simulation box at z = 5 pc. One possibility is that the heating rate Γ is larger in the Galactic central region. When ${{\rm{\Gamma }}}_{\mathrm{GC}}$ = 10⁻²⁵ erg s⁻¹, since the thermal equilibrium condition gives ${\rho }_{\mathrm{WM}}$(z = 5 pc, GC) ∼ 1.4 × 10⁻²³ g cm⁻³, the total mass of the filament can be 5 times that for ${\rho }_{\mathrm{WM}}$(z = 5 pc) = 2.8 × 10⁻²⁴ g cm⁻³. Another possibility is the coexistence of the molecular gas with warm gas at the initial state. According to Sofue & Nakanishi (2016), the fraction of molecular gas f_mol = ${\rho }_{{\rm{H}}2}{V}_{{\rm{H}}2}$/(${\rho }_{\mathrm{WM}}{V}_{\mathrm{WM}}$ + ${\rho }_{{\rm{H}}2}{V}_{{\rm{H}}2}$) ∼ 0.8 at z ∼ 20 pc in the Galactic central region, where ${\rho }_{{\rm{H}}2}$ is the density of molecular hydrogen and V_WM and V_H2 are the volume fractions of the warm gas and molecular hydrogen, respectively. This means that when ${\rho }_{\mathrm{WM}}$ = 2.8 × 10⁻²⁴ g cm⁻³ and ${\rho }_{{\rm{H}}2}$ = 2.8 × 10⁻²² g cm⁻³, the volume fraction of the molecular hydrogen V_H2 ∼ 0.04. Since the mean density of the two-phase medium is ρ = ${\rho }_{\mathrm{WM}}{V}_{\mathrm{WM}}$ + ${\rho }_{{\rm{H}}2}{V}_{{\rm{H}}2}$ = 1.4 × 10⁻²³ g cm⁻³, the total mass of the filament can be 5 times that for the WM with ${\rho }_{\mathrm{WM}}$ = 2.8 × 10⁻²⁴ g cm⁻³. The third possibility is the larger size of the magnetic arch. We assumed that the half width of the magnetic arch is 200 pc. When the magnetic arch is larger, it forms a larger flux rope and more mass accumulates inside the rope, so that the total mass of the dense filament can be larger.

Here we would like to point out that the mass of the dense filament increases with time. Figure 11 shows the time evolution of the total mass of the dense filament. As we have shown in Figure 5, warm gas in the magnetic flux rope continuously falls down toward the filament along magnetic field lines. The maximum mass of the filament can be estimated by assuming that all the gas inside a flux rope falls into the filament. The mass inside the flux rope corresponding to the blue-purple region in Figure 6 is 3.5 × 10⁵ ${M}_{\odot }$, assuming that the depth in the y-direction is 600 pc.

Zoom In Zoom Out Reset image size