{'fraction_non_alphanumeric': 0.048116389645222905, 'fraction_numerical': 0.033080017881090745, 'mean_word_length': 3.7708414113997675, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Magnetohydrodynamic Simulations of a Plunging Black Hole into a Molecular Cloud

Using two-dimensional magnetohydrodynamic simulations, we investigated the gas dynamics around a black hole plunging into a molecular cloud. In these calculations, we assumed a parallel-magnetic-field layer in the cloud. The size of the accelerated region is far larger than the Bondi-Hoyle-Lyttleton radius, being approximately inversely proportional to the Alfv\'en Mach number for the plunging black hole. Our results successfully reproduce the"Y"shape in position velocity maps of the"Bullet"in the W44 molecular cloud. The size of the Bullet is also reproduced within an order of magnitude using a reasonable parameter set. This consistency supports the shooting model of the Bullet, according to which an isolated black hole plunged into a molecular cloud to form a compact broad-velocity-width feature.


INTRODUCTION
To date, ∼ 60 stellar mass black hole (BH) candidates have been detected in our Galaxy by X-ray observations (Corral-Santana et al. 2016), while their total number is estimated to be ∼ 10 8 -10 9 (Agol & Kamionkowski 2002;Caputo et al. 2017). This sharp discrepancy is due to the extremely low percentage of BHs in close binary systems, in which abundant mass accretion from companion stars activate them. Therefore, almost all BHs in our Galaxy still remain undetected.
An isolated, inactive BH would pull up ambient material, leaving a trace in the interstellar medium as a spatially compact, broad-velocity-width feature. The "Bullet" in the W44 molecular cloud is a candidate for such a BH trace (Sashida et al. 2013). W44 is a supernova remnant (SNR) interacting with an adjacent giant molecular cloud (Claussen et al. 1997;Seta et al. 1998Seta et al. , 2004Sashida et al. 2013). In the process of investigating the gas kinematics of the W44 molecular cloud, we noticed an extraordinary broad-velocity-width feature (Sashida et al. 2013), the Bullet. Follow-up observations by Yamada et al. (2017, hereafter Y17) have revealed its compact appearance (0.5 × 0.8 pc 2 ), broadvelocity width nature (∆V ∼ 100), and unique "Y" shape in the position-velocity maps.
The Bullet is intense in CO J=4-3 and HCN J=1-0 emissions, suggesting that it consists of warm and dense molecular gas. The total kinetic energy of the Bullet is ∼ 10 48 erg. This is approximately 1.5 orders of magnitude greater than the kinetic energy of a supernova sharing the small solid angle of the Bullet with respect to the W44 center. Y17 proposed two scenarios of Bullet formation: (1) the expansion model and (2) the shooting model. Both scenarios assume an isolated BH that contributes to the formation of the Bullet. In the expansion model, an additional explosive event triggered by mass accretion onto an isolated BH accounts for the kinematics of the Bullet, but the conversion process from the gravitational energy to the kinetic energy is unclear. The shooting model seems to be more plausible. In this model, the plunge of a 30 M ⊙ BH into the high-density layer toward us successfully explains the broad velocity width as well as the enormous kinetic energy of the Bullet (see Figure 3(b) in Y17). The gas dragged by the plunging BH might correspond to the "Y" shape on the position-velocity map.
One problem is that a native shooting model cannot reproduce the spatial size of the Bullet. This is because the size of the accretion zone, which may correspond to the Bullet size, is described by the Bondi-Hoyle-Lyttleton (BHL) radius (Hoyle & Lyttleton 1939;Bondi & Hoyle 1944;Edgar 2004), R BHL = 2GM BH /(c 2 s +v 2 ), where G, M BH , c s , and v are the gravitational constant, the BH mass, the sound speed, and the velocity of the plunging BH. When M BH ∼ 30 M ⊙ , c s ∼ 1 km s −1 , and v ∼ 100 km s −1 , the BHL radius becomes ∼ 3 × 10 −5 pc, which is too small to reproduce the Bullet size (∼ 0.5pc).
In this paper, we examine the effect of the magnetic field on the size of the accretion zone around a plunging BH. The magnetic field is frozen in the partially ionized interstellar gas, and the bent shape of the field lines propagates with the Alfvén speed. If the BH plunges into the parallel-magnetic-field layer and the magnetic field lines are caught on the gas in the vicinity of the BH, the gas outside of R BHL would be dragged by the plunging BH via the magnetic tension force. This effect may enlarge the size of the accretion zone if the magnetic field is strong enough, which is the case in the W44 expanding shell (B ∼ 500 µG; Hoffman et al. 2005). In this case, the lower part of the "Y" shape in the position-velocity map may correspond to the accelerated gas in front of the BH, which is moving with the plunging speed of BH. The upper part of the "Y" shape can be interpreted as the gas around the BH accelerated from the magnetic tension force, whose velocity may decrease with the distance from the BH.
In order to quantitatively investigate the magnetohydrodynamic (MHD) effect in the shooting model, we simulated the plunging of a BH into a parallel-magneticfield layer by considering a large-scale gas flow around a stationary BH. A two-dimensional MHD code was employed. In Section 2, we explain the calculation method. The results of the simulations are shown in Section 3. Sections 4 and 5 are devoted to discussions and conclusions.

BASIC EQUATIONS AND MODELS
We calculate the gas dynamics around a plunging BH using MHD simulations in Cartesian coordinates (x, y, z). In our simulations, a BH is located at the center of the coordinate system and ambient gas flows in with high velocity. The basic equations of ideal MHD are ∂ρ ∂t where, ρ, v, p, B, r, and e are the mass density, velocity, pressure, magnetic field, distance from the BH, and energy density of the gas written by e = p/(γ − 1) + ρv 2 /2 with γ = 5/3. In our calculations, we assume that the energy dissipation in shock has a time scale comparable to that of the cooling. At each time step, we reset the pressure and the energy density of the gas using an isothermal equation of state, p = ρc 2 s , where we set a sound speed to c s = 0.91 km s −1 , which is consistent with the observed temperature of the Bullet (T = 100 K, Y17) and a mean molecular weight of 1.0. The last term of the equation of motion is the gravitational force due to the BH located at the origin.
The computational domain is −0.5 pc ≤ x ≤ 0.5 pc and −0.5 pc ≤ y ≤ 0.5 pc. We prepare 720 × 720 grids in the domain for M BH = 10 4 M ⊙ so as to set 7-8 grids in the BHL radius, R BHL . The resolution is ∼ 1.4 × 10 −3 pc. We also employ M BH = 10 3 M ⊙ , 10 3.5 M ⊙ , and 10 4.5 M ⊙ . For these cases, we divide the domain into 7200 × 7200, 2160 × 2160, and 240 × 240 grids, corresponding to resolutions of ∼ 1.4 × 10 −4 pc, ∼ 4.6 × 10 −4 pc, and ∼ 4.2 × 10 −3 pc. The choices of M BH larger than 10 3 M ⊙ are for resolving the BHL radii with the realistic number of computational grids. After performing simulations with four M BH , we check the M BH -dependence on the size of the accelerated region (Section 4.1, Figure 7) We show the initial condition of the fiducial model in Figure 1. In the region r ≥ R BHL , the number density and the velocity are set to n = n 0 and v = (0, v y0 , 0). We suppose the situation that a BH plunges into a magnetized layer (high-B layer) with a relative velocity of ∼ 100 km s −1 . We reproduce this situation by shedding high-B layer to a stationary BH, because it is difficult to treat a BH moving in the computational box. Since the high-B layer is accelerated by the magnetic pressure before reaching the BH (see also Section 3.1), we choose the initial velocity (v y0 ) so that the plunge velocity becomes 100 km s −1 . For the fiducial model, the initial velocity is set to v y0 = 50 km s −1 . We employ n 0 = 10 3 cm −3 to simulate the high-B (density) layer after the supernova blast wave. This value is between the typical density in molecular clouds and that in the Bullet (Y17). We also perform the simulations for different n 0 and v y0 (see Section 3.2 and Section 4.1).
In the region r < R BHL , the density is set to n = 0.1n 0 and the velocity to zero. On the upstream side (y ≤ −0.2 pc), the magnetic field is set to B = (B x0 , 0, 0), which describes the high-B layer. We employ B x0 = 500 µG based on the previous measurement in the W44 expanding shell (Hoffman et al. 2005). Sections 3.2 and 4.1 describe the other cases. In the downstream side, the magnetic field is set to B = (50 µG, 0, 0). Figure 1. Initial conditions for MBH = 10 4 M⊙, vy0 = 50 km s −1 , n0 = 10 3 cm −3 , and Bx0 = 500 µG. Color map shows the density in the x-y plane. Black lines and blue arrows show the magnetic field lines and the velocity field. The BH is located at the center of the coordinate system.
We apply the free boundaries at x = −0.5 pc, x = 0.5 pc, and y = 0.5 pc. At the boundary y = −0.5 pc, the density, x-component of the magnetic field, and ycomponent of the velocity are kept constant n = n 0 , B x = B x0 , and v y = v y0 . We impose the free boundary conditions for the other variables. This means that the gas is initially flowing, maintaining the initial density and magnetic field strength in y ≤ −0.5 pc. This situation corresponds to a high-velocity plunging of a BH into a uniform gas cloud. In r < R BHL , the velocity is fixed to zero so as to reproduce the condition that the gas is trapped by the BH and it traps the magnetic field lines. The density is also kept constant at n = 0.1n 0 so that the gas flowing into this region accretes onto the BH and does not appear in the computational domain again.
Numerical simulations are carried out by using the MHD code CANS+ (Matsumoto et al. 2016). The code employs the HLLD approximate Riemann solver (Miyoshi & Kusano 2005). We apply second-order spatial accuracy by using monotone upstream-centered schemes for conservation laws (MUSCL, van Leer 1979). A third order TVD Runge-Kutta scheme is used for solving the time integration. A hyperbolic divergence cleaning method is employed (Dedner et al. 2002). Figure 2 shows the results for M BH = 10 4 M ⊙ , v y0 = 50 km s −1 , n 0 = 10 3 cm −1 , and B x0 = 500 µG (the fiducial model). In this case, the gas is accelerated to v y ∼ 100 km s −1 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 3 yr when the high-B layer passes ∼ 0.4 pc after reaching the BH.

Fiducial Model
In the region |x| 0.36 pc, the gas flows as it is not affected by the magnetic tension force. Near the BH (|x| 0.36 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 positionvelocity (x-v y ) map. We divide the v y -axis into grids of the width ∆v y = 5 km s −1 . 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| 0.36 pc), the bulk of the gas has a velocity of ∼ 100 km s −1 , because of less influence of the magnetic field at this time. In the region |x| 0.36 pc, another velocity component appears at v y 100 km s −1 . The v y ∼ 100 km 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 −1 ) makes a "Y" shape on the x-v y map. The widths of the "Y" shape in the xdirection (d) and in the v y -direction (∆v y ) are ∼ 0.72 pc and ∼ 115 km s −1 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 −1 ), the time is converted to the width of the layer that passed through the BH, L. The size increases in proportion to and ∆v y ∼ 100 km s −1 respectively. In Figure 4(b), these values are d ∼ 0.78 pc and ∆v y ∼ 130 km s −1 .
These results show that the position-velocity structure is almost independent of the BH mass (see also Figure  2). 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 −1 , which is the same as the other cases, we set the initial velocity to v y0 = 90 km 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| 0.15 pc. For this parameter set, the magnetic field lines are sharply bent by the ram pressure of the flow, since 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 km s −1 to +15 km s −1 . 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 µG and n 0 = 100 cm −3 . 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. 4. DISCUSSIONS

Size of Accelerated Gas
Figures 5(a) and 5(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 A = B x0 / √ 4πρ 0 , when L = 0.2 pc and 0.4 pc. In order to survey the v A -dependence, we employ n 0 = 10 5 cm −3 , 10 4 cm −3 , 10 3 cm −3 , and 100 cm −3 while the magnetic field strength is set to B x0 = 500 µ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 0 = 10 3 cm −3 (blue triangles). Here, the BH mass is kept constant M BH = 10 4 M ⊙ . The initial velocity v y0 is adjusted to set the inflow velocity to v in ∼ 100 km s −1 before the high-B layer reaches the BH. Both  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 −1 , whose corresponding initial velocities are v y0 = 0, 50, 100, and 150 km s −1 (black squares). In these calculations, M BH , n 0 , 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/pc = 72/(v in /km 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. 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 3 M ⊙ , 10 3.5 M ⊙ , 10 4 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 3 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 ⊙ .
From Figures 3, 5(a), 5(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

Comparison with Observations
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 highvelocity 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 −1 and the sizes of ∼ 1 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 A L lines employing the relation d = 5.1L/M 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 A L 20, which can be rewritten as 5(L/pc) −1 (n/cm −3 ) 1/2 B/µG 9(L/pc) −1 (n/cm −3 ) 1/2 . The value of v A L can also be estimated from the observations. In the case of the Bullet, employing B = 500 µG (Hoffman et al. 2005), n = 10 4 cm −3 (Y17), and L = 0.1 pc that is the thickness of thin filaments detected in W44 (Jones et al. 1993), we obtain v A L ∼ 1.1 pc km s −1 . 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 ∼ 500 µG (e.g., Morris & Serabyn 1996) and the density of the molecular cloud is n ∼ 10 5 cm −3 (Güsten & Henkel 1983). We assume that L is comparable to the HVCC size and is smaller than the size of the molecular cloud, L ∼ 1-5 pc. As a result, we obtain v A L ∼ 3.5-17.3 pc km s −1 . 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/M 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 ∼ R BHL , because we ignore the dynamics in r < R BHL and assume the simple conditions n = 0.1n 0 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/M 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 ∼ 0.5 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 BH = 10 M ⊙ , n = 10 4 cm −3 , and v in = 100 km s −1 , the total luminosity is estimated to be L X = 0.06Ṁ BHL c 2 ∼ 2 × 10 34 erg s −1 , whereṀ 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Ṁ 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 (RIAF, 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 is 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 must be important and interesting future work.
In our simulations, we assumed a simple parallelmagnetic-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.

CONCLUSIONS
Performing the MHD simulations, we investigated the gas dynamics around a BH plunging into a molecular cloud with a parallel-magnetic-field. We found the following results: 1. The MHD effects enlarge the accelerated region compared to the native shooting model, and the acceleration region expands within |x| 0.36 pc when L = 0.4 pc for the fiducial model (M BH = 10 4 M ⊙ , B x0 = 500 µG, n 0 = 10 3 cm −1 , and v in = 100 km s −1 ).
2. When L = 0.4 pc, the accelerated gas exhibits a "Y" shape with a size of d ∼ 0.72 pc and a velocity width of ∆v ∼ 115 km s −1 on the x-v y map for the fiducial model.
3. The size of the accelerated gas increases in proportion to the time (the distance traveled in the layer, L).
4. Our simulations show that the size of the "Y" staple is almost independent of the BH mass and stellar mass BH (M BH 100 M ⊙ ) is preferred for the Bullet if we assume the BHL accretion rate and the standard accretion disk.
5. The size of the "Y" shape is approximately determined by d = 5.1L/M A . 6. Our model can reproduce the "Y" shape and the velocity width on the position-velocity map of the W44 Bullet.
7. The size of the "Y" shape expected from the model is consistent with that of the Bullet within one order of magnitude. 8. Our model can explain the velocity width and the size of the small HVCCs in the Galactic center.
The foregoing results support the shooting model of the Bullet and the small HVCCs and indicate that the MHD effects are necessary to reproduce the size of the accelerated gas. There is a possibility that the plunging of the isolated BH into the molecular cloud be responsible for the formation of the extraordinary high velocity component.