Magnetohydrodynamic Simulations of a Plunging Black Hole into a Molecular Cloud

Figure 2 shows the results for M_BH = 10⁴ M_⊙, ${v}_{y0}=50\,\mathrm{km}\,{{\rm{s}}}^{-1}$, n₀ = 10³ cm⁻¹, and B_x0 = 500 μG (the fiducial model). In this case, the gas is accelerated to v_y ~ 100 km s⁻¹ by the magnetic pressure before the high-B layer reaches the BH. This velocity corresponds to the plunging speed of the BH in the shooting model (hereafter, we call it the inflow velocity, v_in). The top panel shows the density map at t = 5.38 × 10³ years when the high-B layer passes $\sim 0.4\,\mathrm{pc}$ after reaching the BH.

Zoom In Zoom Out Reset image size

In the region $| x| \gtrsim 0.36\,\mathrm{pc}$, the gas flows as it is not affected by the magnetic tension force. Near the BH ($| x| \lesssim 0.36\,\mathrm{pc}$), the high-density and low-velocity region (accelerated region) appears and forms an arcuate shape at the side facing the BH (see the yellow region at −0.2 pc y 0.2 pc). At the front of the BH, gas is compressed and the velocity is decelerated (accelerated toward the negative direction). In this area, the magnetic field is enhanced to ~600–3000 μG. This is because gas is frozen into the magnetic field and the field lines are caught on the gas dammed up in r < R_BHL. At the front of the BH, the magnetic pressure strongly accelerates the gas, and at both sides of the BH, the magnetic tension force due to the curved magnetic field lines contributes to the acceleration. The following flow collides with the accelerated flow inducing shock.

The bottom panel of Figure 2 shows a position–velocity (x–v_y) map. We divide the v_y-axis into grids of the width Δ v_y = 5 km s⁻¹. The color represents the column density of the gas in each velocity grid. The column density is integrated along the y-axis over the range −0.25 pc ≤ y ≤ 0.25 pc (between the magenta dashed lines in the top panel of Figure 2) in order to focus on the velocity structure around the BH.

In the region far from the BH ($| x| \gtrsim 0.36\,\mathrm{pc}$), the bulk of the gas has a velocity of ~100 km s⁻¹, because of less influence of the magnetic field at this time. In the region $| x| \lesssim 0.36\,\mathrm{pc}$, another velocity component appears at v_y  100 km s⁻¹. The ${v}_{y}\sim 100\,\mathrm{km}\,{{\rm{s}}}^{-1}$ component is contributed by the less-dense gas in the upstream side (y  −0.18 pc). The low-velocity component (v_y < 100 km s⁻¹) makes a "Y" shape on the x–v_y map. The widths of the "Y" shape in the x-direction (d) and in the v_y-direction (${\rm{\Delta }}{v}_{y}$) are ~0.72 pc and ~115 km s⁻¹, respectively. This component is contributed by the dense gas in the accelerated region around the BH. The v_y < 0 component is contributed by the gas in the vicinity of the BHL radius, which falls into the BH along the magnetic field.

Figure 3 shows the time evolution of the size of the accelerated region, which corresponds to the size of the "Y" shape on the x–v_y map, d. Here, Δt is the time measured from the moment that the high-B layer reached the BH. Using the speed of the flow (~100 km s⁻¹), the time is converted to the width of the layer that passed through the BH, L. The size increases in proportion to Δt and L. The best-fitting line (the solid black line) is d/pc = 2.05 × 10⁻⁴ Δt yr⁻¹ (d = 1.8L). This shows that the accelerated region expands with the Alfvén speed at the vicinity of R_BHL, ~200 km s⁻¹, which is roughly 5.8 times the Alfvén speed in the high-B layer.

Zoom In Zoom Out Reset image size

Each panel of Figure 4 is the same as that of Figure 1, but the parameters are different from the fiducial model. Figures 4(a) and (b) show the results in the case of M_BH = 10^3.5 M_⊙ and 10^4.5 M_⊙. In both figures, the overall structures are quite similar to the result of the fiducial model. In Figure 4(a), the size and the velocity width of the "Y" shape in the x–v_y map are d ~ 0.70 pc and Δv_y ~ 100 km s⁻¹, respectively. In Figure 4(b), these values are d ~ 0.78 pc and Δv_y ~ 130 km s⁻¹. These results show that the position–velocity structure is almost independent of the BH mass (see also Figure 2).

Zoom In Zoom Out Reset image size

Figure 4(c) shows the results for the weak magnetic field, B_x0 = 158 μG. In order to adjust the inflow velocity to v_in ~ 100 km s⁻¹, which is the same as the other cases, we set the initial velocity to ${v}_{y0}=90\,\mathrm{km}\,{{\rm{s}}}^{-1}$. The other parameters are the same as those of the fiducial model. The density map shows that the accelerated region around the BH is small, $| x| \lesssim 0.15\,\mathrm{pc}$. For this parameter set, the magnetic field lines are sharply bent by the ram pressure of the flow, as the magnetic field strength is weak. As a result, in the x–v_y map, the size of the "Y" shape is d ~ 0.3 pc. The velocity width of the "Y" shape is almost the same as that of the fiducial model, but the matter located at x ~ 0 is concentrated in the narrow velocity range −15 to +15 km s⁻¹. This is because the gradient of the magnetic field strength is large at the front of the BH, and the flow is rapidly decelerated.

Figure 4(d) is the same as Figure 2 except that ${B}_{x0}=158\,\mu {\rm{G}}$ and n₀ = 100 cm⁻³. The size of the accelerated region is larger than that of Figure 4(c) and comparable to that of Figure 2. In this case, the magnetic field is weak, but the ram pressure of the flow is also small due to the low density. As a consequence, the magnetic field lines draw gentle curves near the BH. In the x–v_y map, the column density in each cell is smaller than that of the fiducial model, but the position–velocity structure is very similar to that of Figure 2.

Figures 5(a) and (b) show the size of the "Y" shape on the x–v_y map, d, as a function of the Alfvén speed in the high-B layer, ${v}_{{\rm{A}}}={B}_{x0}/\sqrt{4\pi {\rho }_{0}}$, when L = 0.2 and 0.4 pc. In order to survey the v_A-dependence, we employ n₀ = 10⁵ cm⁻³, 10⁴ cm⁻³, 10³ cm⁻³, and 100 cm⁻³ while the magnetic field strength is set to ${B}_{x0}=500\,\mu {\rm{G}}$ (red squares). In addition, we investigate the position–velocity structures for B_x0 = 50 μG, 158 μG, 500 μG, and 1.58 mG when the density is n₀ = 10³ cm⁻³ (blue triangles). Here, the BH mass is kept constant ${M}_{\mathrm{BH}}={10}^{4}\,{M}_{\odot }$. The initial velocity v_y0 is adjusted to set the inflow velocity to v_in ~ 100 km s⁻¹ before the high-B layer reaches the BH. Both Figures 5(a) and (b) show that the sizes are similar to each other if the Alfvén speeds have the same value, regardless of the combination of the density and the magnetic field strength. The size increases in proportion to the Alfvén speed and the best-fitting lines are $d/\mathrm{pc}\,=0.012{v}_{{\rm{A}}}/\mathrm{km}\,{{\rm{s}}}^{-1}$ and d/pc = 0.021v_A/km s⁻¹ for L = 0.2 pc and 0.4 pc (the solid black lines).

Zoom In Zoom Out Reset image size

Figure 6 shows the size, d, as a function of the inverse of the inflow velocity, 1/v_in, when L = 0.4 pc. We employ four different inflow velocities, v_in = 50, 100, 150, and 200 km s⁻¹, whose corresponding initial velocities are ${v}_{y0}=0$, 50, 100, and 150 km s⁻¹ (black squares). In these calculations, ${M}_{\mathrm{BH}}$, n₀, and B_x0 are the same as those of the fiducial model. The size increases as the inflow velocity decreases and this relation is fitted by $d/\mathrm{pc}=72/({v}_{\mathrm{in}}/\mathrm{km}\,{{\rm{s}}}^{-1})$ (the solid black line). The decrease of the inflow velocity suppresses the ram pressure of the flow. This leads the magnetic field lines drawing gentle curves and the large accelerated region.

Zoom In Zoom Out Reset image size

Figure 7 shows the M_BH-dependence of the size, d, when L = 0.4 pc. We calculate the position–velocity structures for M_BH = 10³ M_⊙, 10^3.5 M_⊙, 10⁴ M_⊙, and 10^4.5 M_⊙ (black squares), while the other parameters are same as those of the fiducial model. This figure shows that the size hardly depends on the BH mass at least in the range 10³ M_⊙ ≤ M_BH ≤ 10^4.5 M_⊙. The best-fitting line is d/pc = 0.51(M_BH/M_⊙)^0.040 (the solid black line). If we extrapolate this relation in the low-mass range, a "Y"-shaped structure with a size of ~0.5 pc is produced by a BH of M_BH ~ 10 M_⊙.

Zoom In Zoom Out Reset image size

From Figures 3, 5(a) and (b), 6, and 7, we find that the size, d, depends mainly on v_A, v_in, and L, besides having a weak dependence on M_BH. The size is approximately determined by $d={aL}/{{ \mathcal M }}_{{\rm{A}}}$, where a is a proportionality constant, and ${{ \mathcal M }}_{{\rm{A}}}$ is the Alfvén Mach number in the high-B layer, ${{ \mathcal M }}_{{\rm{A}}}$ = ${v}_{\mathrm{in}}/{v}_{{\rm{A}}}$. We derive a = 5.1 from the results of M_BH = 10⁴ M_⊙. The size of $L/{{ \mathcal M }}_{{\rm{A}}}$ is the scale of the warped magnetic field lines when the high-B layer passed L after reaching the BH. This is derived from the balance between the magnetic tension force and the ram pressure of the flow. The proportionality constant, $a(\gt 1)$, shows an enhanced Alfvén speed and a decreased velocity around the BH.

Our results show that the plunging of the BH into the high-B layer reproduces the characteristic "Y" shape on the position–velocity map of the Bullet in the W44 SNR. Here, we quantitatively compare our results to the two objects, the W44 Bullet (Y17) and the small high-velocity compacts clouds (HVCCs) detected near the Galactic nucleus (HCN–0.009–0.044 and HCN–0.085–0.094; Takekawa et al. 2017, hereafter T17) The small HVCCs have the velocity widths of Δ v 60 km s⁻¹ and the sizes of $\sim 1\,\mathrm{pc}$. The high-velocity components originate from the dense molecular clouds. Although the "Y" shape is not resolved, the position–velocity structure of the small HVCCs is similar to that of the W44 Bullet.

In Figure 8, we compare the theoretical predictions based on our results and the observations on the d–v_in plane. We plot the constant v_AL lines employing the relation $d=5.1L/{{ \mathcal M }}_{{\rm{A}}}$ (the solid lines) as well as the size and the velocity width of the Bullet and the small HVCCs (gray regions). Here, on the basis of our simulations, we assume that the velocity width is comparable to the inflow velocity. We find that both the Bullet and the small HVCCs are located in the range 10  v_AL  20, which can be rewritten as $5{(L/\mathrm{pc})}^{-1}{(n/{\mathrm{cm}}^{-3})}^{1/2}\lesssim B/\mu {\rm{G}}\lesssim 9{(L/\mathrm{pc})}^{-1}{(n/{\mathrm{cm}}^{-3})}^{1/2}$.

Zoom In Zoom Out Reset image size

The value of v_AL can also be estimated from the observations. In the case of the Bullet, employing B = 500 μG (Hoffman et al. 2005), n = 10⁴ cm⁻³ (Y17), and $L=0.1\,\mathrm{pc}$, which is the thickness of thin filaments detected in W44 (Jones et al. 1993), we obtain v_AL ~ 1.1 pc km s⁻¹. This is consistent with the value predicted by our results within an order of magnitude. In the case of the small HVCCs, the typical magnetic field strength in the central region of our Galaxy is $B\sim 500\,\mu {\rm{G}}$ (e.g., Morris & Serabyn 1996) and the density of the molecular cloud is n ~ 10⁵ cm⁻³ (Guesten & Henkel 1983). We assume that L is comparable to the HVCC size and is smaller than the size of the molecular cloud, $L\sim 1$–5 pc. As a result, we obtain v_AL ~ 3.5–17.3 pc km s⁻¹. This is comparable to or slightly smaller than the theoretical prediction.

The difference between our results and the observations might be caused by the uncertainty of the coefficient, a ~ 5.1, in the relation $d={aL}/{{ \mathcal M }}_{{\rm{A}}}$. The coefficient, a, approximately corresponds to the ratio of the Alfvén speed at r ~ R_BHL to that of the high-B layer at the initial condition. It is difficult to accurately evaluate the Alfvén speed in the vicinity of $r\sim {R}_{\mathrm{BHL}}$, because we ignore the dynamics in r < R_BHL and assume the simple conditions n = 0.1n₀ and v_y = 0. The magnetic field strength and the density near the BH would change if we consider the realistic magnetic structure around the BH and the feedback from the accretion flow such as radiative heating. In that case, there is the possibility that the coefficient, a, increases by several times to ten times, and the size of the Bullet would be well explained by the relation $d={aL}/{{ \mathcal M }}_{{\rm{A}}}$.

We found that the size of the accelerated region is almost independent of the BH mass (Figure 7). This result indicates that "Y"-shaped position–velocity structure with $d\sim 0.5\,\mathrm{pc}$ can be reproduced by the plunging of a stellar mass BH. Here, we estimate the total luminosity of the BH based on the BHL accretion model. Assuming ${M}_{\mathrm{BH}}=10\,{M}_{\odot }$, $n={10}^{4}\,{\mathrm{cm}}^{-3}$, and ${v}_{\mathrm{in}}=100\,\mathrm{km}\,{{\rm{s}}}^{-1}$, the total luminosity is estimated to be ${L}_{{\rm{X}}}=0.06\,{\dot{M}}_{\mathrm{BHL}}{c}^{2}\sim 2\times {10}^{34}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$, where ${\dot{M}}_{\mathrm{BHL}}$ is the BHL accretion rate (e.g., Edgar 2004). This luminosity is consistent with the absence of the X-ray counterpart in ROSAT All Sky Survey (RASS; Haakonsen & Rutledge 2009). If we assume the BHL accretion rate and standard accretion disk (Shakura & Sunyaev 1973), non-detection of the X-ray counterpart in the RASS indicates that the mass of the plunging BH is less than ~100 M_⊙. Note that the accretion rate would be smaller than ${\dot{M}}_{\mathrm{BHL}}$ because the magnetic field suppresses the mass accretion rate (Lee et al. 2014). In addition, the accretion disk might become radiation inefficient accretion flow (Narayan & Yi 1994). In such cases, a BH mass larger than 100 M_⊙ is acceptable. In order to calculate the accurate accretion rate and corresponding X-ray luminosity, the three-dimensional MHD and/or radiation hydrodynamics simulations in the small scale around the BH are necessary. Detection of an X-ray counterpart with modern X-ray imaging telescopes (such as Chandra) will provide a strong support for our scenario.

As a first step of the theoretical approach to the Bullet, we performed the two-dimensional simulations on the x–y plane (z = 0 plane). In three-dimensional simulations, the gas and magnetic field behaviors in the z = 0 plane are expected to be similar to those in the two-dimensional simulations. This is because the gas and magnetic field lines far from the z = 0 plane (z R_BHL) do not affect those in the z = 0 plane. The three-dimensional MHD simulations of BHL accretion in a small computational domain show the bow shock similar to our results (Lee et al. 2014), supporting that our two-dimensional simulations, at least qualitatively, well reproduce the gas and magnetic field behaviors in the z = 0 plane. In the three-dimensional simulations, the magnetic pressure in the z-direction might lower the density and magnetic field strength near the BH, and thereby reduce the size of the accelerated region in z = 0 plane. Investigating these three-dimensional effects will be important and interesting future work.

In our simulations, we assumed a simple parallel-magnetic-field layer, although magnetic fields in molecular clouds are not highly ordered. The unordered fields may produce the asymmetric "Y" shape, and the size of the accelerated region might slightly change. More accurate measurements of the magnetic field configuration near the Bullet and the simulations employing realistic settings may be interesting future works.