Structure and Instability of the Ionization Fronts around Moving Black Holes

Here, we briefly review the analytical model developed in PR13. This quasi-1D model of flow around a moving BH with radiation feedback reproduces the variation of the type of flow and the resulting accretion rate as a function of BH velocity observed in their simulations. We refer the readers to PR13 for the details of this model.

The model considers the structure of flow along the axis of the BH motion. In the BH rest frame, the gas is moving toward the BH with velocity $v_{\infty}$. Hereafter, the subscripts $\mathrm{I}$ and $\mathrm{II}$ indicate the neutral and ionized sides of the I-front, respectively. For simplicity, we assume that the sound speed (or temperature) discontinuously changes from $c_{\mathrm{I}}$ in the neutral gas upstream of the I-front to $c_{\mathrm{II}}$ in the ionized gas across the I-front. In this section, we assume that the BH velocity with respect to the neutral ambient gas is supersonic (i.e., $v_{\infty}\geq c_{\infty}$), and the ambient gas has constant temperature, i.e., $c_{\mathrm{I}}=c_{\infty}$. We consider the subsonic case ($v_{\infty}<c_{\infty}$), as well as the effect of X-ray pre-heating of the ambient gas ($c_{\mathrm{I}}\not=c_{\infty}$), in Appendix A.

Let us begin with writing down the jump conditions for density $\rho$ and velocity $v$ across the I-front,

In order for Eq. (1) to have a (real) solution, the condition $v_{\mathrm{I}}<v_{\mathrm{D}}$ or $v_{\mathrm{I}}>v_{\mathrm{R}}$ must be satisfied. Here, the D-type critical velocity is defined as

and the R-type critical velocity is

The I-fronts with $v_{\mathrm{I}}<v_{\mathrm{D}}$ and $v_{\mathrm{I}}>v_{\mathrm{R}}$ are called D-type and R-type, respectively.

In the case with $v_{\infty}>v_{\mathrm{R}}$, the gas can directly reach the I-front with $v_{\mathrm{I}}=v_{\infty}$ and $\rho_{\mathrm{I}}=\rho_{\infty}$. The I-front becomes weak R-type (i.e., the minus sign is taken in Eq. 1), and the density and velocity inside the HII region are given, respectively, as

and

Here, the density and velocity jump is generally weak because $\Delta^{(-)}$ takes its maximum value of $2$ $(c_{\mathrm{I}}\ll c_{\mathrm{II}})$ at $v_{\infty}=v_{\mathrm{R}}$ and approaches unity as $v_{\infty}$ increases. Hereafter, we name this type of flows as the R-type flows.

In the case with $(v_{\mathrm{D}}<c_{\mathrm{I}}<)\ v_{\infty}<v_{\mathrm{R}}$, however, the jump conditions in Eq. (1) cannot be satisfied without extra structures. In this case, a shock is generated in front of the I-front. The gas experiencing the shock forms a dense shell between the shock and the I-front. The I-front becomes D-type because the gas slows down to below $v_{\mathrm{D}}$ as a result of the following two effects. Firstly, the isothermal shock decreases the velocity by $(v_{\infty}/c_{\mathrm{I}})^{2}$, while increasing the density by the same factor. Secondly, the tangentially diverging motion decelerates the gas in the shell (we will discuss this process in detail in Sec. 2.2). In PR13, they further assumed that the I-front is critical D-type (i.e., $v_{\mathrm{I}}=v_{\mathrm{D}}$), and thus the velocity inside the HII region is given by

Here, the total pressure $P_{\mathrm{tot}}=\rho(v^{2}+c_{\mathrm{s}}^{2})$, defined as the sum of the ram pressure $P_{\mathrm{ram}}=\rho v^{2}$ and the thermal pressure $P_{\mathrm{th}}=\rho c_{\mathrm{s}}^{2}$, is conserved across the shock and the I-front. It is also conserved inside the shocked shell, where the flow is subsonic, because $P_{\mathrm{tot}}\approx P_{\mathrm{th}}$ and $P_{\mathrm{th}}$ is conserved due to the approximate pressure equilibrium. From Eq. (6) and $P_{\mathrm{tot}}=\rho_{\infty}(v_{\infty}^{2}+c_{\mathrm{I}}^{2})$, the density inside the HII region is given by

Hereafter, we name this type of flows as the D-type flows.

Using the density and velocity inside the HII region obtained above, the accretion rate can be estimated with the Bondi-Hoyle-Lyttleton (BHL) formula as

where the accretion rate for the D-type flows

and that for the R-type flows

Here, we have used $\rho_{\mathrm{II}}$ and $v_{\mathrm{II}}$ given by Eqs. (4) and (5) for $\dot{M}_{\mathrm{R}}$ and those given by Eqs. (6) and (7) for $\dot{M}_{\mathrm{D}}$. In Fig. 1, we plot the accretion rate given by Eq. (8), with the temperature of the neutral gas $T_{\mathrm{I}}=10^{4}\,\mathrm{K}$ ($c_{\mathrm{I}}=8.1\,\mathrm{km/s}$) and that of the ionized gas $T_{\mathrm{II}}=4\times 10^{4}\,\mathrm{K}$ ($c_{\mathrm{II}}=23\,\mathrm{km/s}$). We use the same $T_{\mathrm{I}}$ and $T_{\mathrm{II}}$ for the Park & Ricotti model in the rest of the paper unless otherwise stated. As pointed out in PR13, $\dot{M}$ takes its maximum value at $v_{\infty}=v_{\mathrm{R}}$ in Figure 1.

For $\dot{M}$ obtained above, the luminosity of the X-ray emission from the BH accretion disk can be calculated as

with radiative efficiency $\eta$. Here, we assume an $\dot{M}$-dependent efficiency modeled as

with the normalization factor $\eta_{0}=0.1$. The Eddington luminosity is defined as $L_{\mathrm{Edd}}\equiv 4\pi GM_{\mathrm{BH}}cm_{\mathrm{p}}/\sigma_{\mathrm{T}}$, where $m_{\mathrm{p}}$ is the proton mass and $\sigma_{\mathrm{T}}=6.65\times 10^{-25}\,\mathrm{cm^{2}}$ is the Thomson cross-section. This form of $\eta$ is motivated by the fact that the radiative efficiency becomes lower with lower $\dot{M}$ in the in the low accretion rate RIAF regime (Narayan & Yi, 1994).

In order to better understand the stability of the I-front, we estimate its size $R_{\mathrm{IF}}$, using the luminosity and density modeled above. As pointed out in PR13, $R_{\mathrm{IF}}$ is either determined by the balance of the ionization and recombination rates inside the ionized region, or by the balance of the ionization rate and the flow rate of neutral gas across the I-front. In the former case, $R_{\mathrm{IF}}$ is well approximated by the Strömgren radius,

where $\dot{N}_{\mathrm{ion}}=L/\left<h\nu\right>$ is the ionizing photon emission rate, with average energy of ionizing photons $\left<h\nu\right>$, $n_{\mathrm{II}}\approx\rho_{\mathrm{II}}/m_{\mathrm{p}}$ is the ionized gas number density, and $\alpha_{\mathrm{B}}\approx 8.3\times 10^{-14}\,\mathrm{cm^{3}\,s^{-1}}\,(T_{\mathrm{II}}/4\times 10^{4}\mathrm{K})^{-1}$ is the Case-B recombination coefficient (Ferland et al., 1992). Alternatively, in the latter case, $R_{\mathrm{IF}}$ is approximated by

where we have used the particle number conservation $n_{\mathrm{I}}v_{\mathrm{I}}=n_{\mathrm{II}}v_{\mathrm{II}}$ across the I-front. The size $R_{\mathrm{IF}}$ is determined by the smaller between $r_{\mathrm{Strm}}$ and $r_{\mathrm{neutral-flow}}$. For the cases explored in our simulations, generally we have $r_{\mathrm{Strm}}\leq r_{\mathrm{neutral-flow}}$, and thus $R_{\mathrm{IF}}\sim r_{\mathrm{Strm}}$.

In the RIAF regime ($\dot{M}<L_{\mathrm{Edd}}/c^{2}$; see Eq. 12), with Eqs. (5), (11), (13), and (14), we can show that the condition $r_{\mathrm{Strm}}>r_{\mathrm{neutral-flow}}$ is satisfied if

Here, we have assumed $T_{\mathrm{II}}=4\times 10^{4}\,\mathrm{K}$ and $\left<h\nu\right>\approx 41\,\mathrm{eV}$ for the power-low radiation spectrum $L_{\nu}\propto\nu^{-1.5}$ with a lower cut-off at $13.6\,\mathrm{eV}$. Therefore, for velocities higher than the critical velocity in Eq. (15), $R_{\mathrm{IF}}\approx r_{\mathrm{neutral-flow}}$, while for velocities lower than the critical value from Eqs. (8) and (11) – (13) we find:

with the BHL radius for the ionized gas defined as

The ratio of $R_{\mathrm{IF}}$ to $r_{\mathrm{BHL,II}}$ is independent of either $M_{\mathrm{BH}}$, $n_{\infty}$, or $v_{\infty}$, except for the dependence through $T_{\mathrm{II}}$. Eq. (16) implies that the BHL radius is always much smaller than the size of the I-front in the RIAF regime.

In the following, we extend the Park & Ricotti model to understand the property of the shells forming in the D-type flows. Specifically, we model the thickness of the shell $\Delta R_{\mathrm{shell}}$, which is a critical parameter for the stability of the I-front, as we will see in Sec. 4. To model $\Delta R_{\mathrm{shell}}$, we start from mass conservation in a cylindrical region of the shell with radius $b$ from the axis of the BH motion (see Fig. 2 a). Below, we first estimate the mass flux at each surface of the cylinder and then derive the expression for the shell thickness. In this section, rather than restricting the model to isothermal shocks, we allow the sound speed in the shell, $c_{\mathrm{shell}}$, to differ from that of the incident flow $c_{\mathrm{I}}$. However, for simplicity, we assume that the temperature of the shell is kept constant at $\approx 10^{4}\,\mathrm{K}$ due to the strong temperature dependence of Ly$\alpha$ cooling (see also Appendix A.2), with inefficient cooling below $10^{4}$ K expected in a low-metallicity gas. We also neglect molecular hydrogen formation and cooling in the shell.

Let us start with obtaining the incoming mass flux across the shock

and outgoing mass flux across the I-front $\dot{m}_{\mathrm{IF}}=\rho_{\mathrm{II}}v_{\mathrm{II}}$. Using Eqs. (6) and (7), we rewrite the mass flux in the HII region as

Note that the expression for $\rho_{\mathrm{II}}$ (Eq. 7), which is derived from total pressure equilibrium and conservation of mass flux, is valid even if the temperature changes across the shock.

Next, we estimate the outgoing mass flux at the lateral surface of the cylinder. Here, we make a rough approximation and regard the bow-like shell as a thin disk with radius $B\sim O(R_{\mathrm{IF}})$, neglecting the curvature of the shell and the structure along the axis (Fig. 2b). The thermal pressure of the shocked gas, $P_{\mathrm{shell}}$, which is higher than that of the ambient gas, $P_{\infty}$, causes the radial expansion of the disk. We estimate the gas velocity at the outer edge as $V_{\perp}\sim O(v_{\infty})$, considering the approximate balance between the thermal pressure of the shell $P_{\mathrm{shell}}\approx\rho_{\infty}v_{\infty}^{2}$ and the ram pressure of the ambient gas $\rho_{\infty}V_{\perp}^{2}$ in the rest frame of the edge. Since the vertical velocity $v_{\perp}$ is proportional to $b$ in a uniformly expanding disk, we model it as

where we introduce an $O(1)$ numerical factor $\alpha$ in the last expression to absorb the uncertainties of the model. The density of the shell is determined by the approximate balance between the inner thermal pressure and outer ram pressure at the shock as

Therefore, from Eqs. (20) and (21), the vertical mass flux $\dot{m}_{\perp}=\rho_{\mathrm{shell}}\,v_{\perp}$ at radius $b$ is given by

Finally, we obtain $\Delta R_{\mathrm{shell}}$ considering mass conservation. The equation for the conservation of mass inside the cylinder with radius $b$ and length $\Delta R_{\mathrm{shell}}$ is given by

where the first, second, and third terms correspond to the flows through the shock, I-front, and lateral surface of the cylinder, respectively. With Eqs. (19), (18) and (22), Eq. (23) can be solved for $\Delta R_{\mathrm{shell}}$ as

where we have made an approximation $c_{\mathrm{I}}\ll v_{\infty}$ and used the approximate expression of the critical velocity $v_{\mathrm{R}}\approx 2c_{\mathrm{II}}$ (Eq. 3) in the last equation. Later, we will determine the $O(1)$ value for $\alpha$ and the effective temperature of ionized region $T_{\mathrm{II}}$ by comparing Eq. (24) to the simulation results. Here we define $T_{\mathrm{II}}$ as an effective temperature because it absorbs the uncertainties coming from the temperature variation inside the HII region and the simplifying assumptions in modelling $\dot{m}_{\mathrm{IF}}$. This expression implies that $\Delta R_{\mathrm{shell}}$ approaches zero as $v_{\infty}$ approaches $v_{\mathrm{R}}$ from the lower side. Thus, we can consider that the transition from the D- to R-type flow occurs at $v_{\mathrm{R}}$, not only because the R-type I-fronts are allowed only for $v_{\infty}>v_{\mathrm{R}}$, but also because the shells in the D-type flows can have finite thickness only for $v_{\infty}<v_{\mathrm{R}}$.

Before closing this section, we shortly discuss the column density of the shell. From Eqs. (21) and (25), we can roughly estimate it as $m_{\mathrm{p}}\,\Delta N_{\mathrm{shell}}=\rho_{\mathrm{shell}}\,\Delta R_{\mathrm{shell}}\sim\rho_{\infty}\,R_{\mathrm{IF}}\,(1-v_{\infty}/v_{\mathrm{R}})$, omitting the $O(1)$ numerical factors. Using Eqs. (16) and (17), for velocities $v_{\infty}<v_{\mathrm{R}}$, we obtain the shell column density:

which is proportional to $n_{\infty}\,M_{\mathrm{BH}}$ with the dependence on $v_{\infty}$ through the factor $(1-v_{\infty}/v_{\mathrm{R}})$. We see that the shell is typically optically thick to ionizing radiation for $n_{\infty}\,M_{\mathrm{BH}}>10^{3}\,M_{\odot}\,\mathrm{cm^{-3}}$.

Figure 3 shows the snapshot of the steady-state flow for the run with $n_{\infty}=10^{5}\,\mathrm{cm^{-3}}$, $M_{\mathrm{BH}}=10^{2}\,M_{\odot}$, $T_{\infty}=10^{4}\,\mathrm{K}$, and $v_{\infty}/c_{\infty}=2$. We perform this run as the reference case for the stable D-type flows in the Region I of Figure 1. As expected, we see that a stable, dense shell forms between the D-type I-front and the preceding shock, alike a bow shock around a blunt body (see, e.g., Yalinewich & Sari, 2016; Keshet & Naor, 2016, for recent studies).

Let us investigate the structure of the flow in detail. In the HII region, the gas is heated to the equilibrium temperature $T_{\mathrm{II}}\approx 4\,\text{--}\,5\times 10^{4}\,\mathrm{K}$, determined by the balance of the photo-ionization heating and the Ly$\alpha$ and free-free cooling. The shock is isothermal due to the efficient Ly$\alpha$ cooling in the neutral gas, and the density jump in the shell is $(v_{\infty}/c_{\infty})^{2}\approx 4$ of the ambient value. As considered in the analytical model in Sec. 2, the gas motion is approximately plane-parallel except for inside the shell, where the tangentially diverging motion has a significant effect on the streamlines. The shell is rather thick ($\Delta R_{\mathrm{shell}}/R_{\mathrm{IF}}\sim 0.1$) and stable. The size of the I-front, $R_{\mathrm{IF}}\sim 2\times 10^{4}\,\mathrm{au}$, agrees with the analytic Strömgren radius in Eq. (13). In general, the flow structure is consistent with previous 2D simulations in PR13 and agrees with the analytical model.

To understand the properties of the shell and its stability, we investigate the dependence of the shell thickness on the BH velocity, by performing runs with various BH velocities $v_{\infty}/c_{\infty}=1.5\,\text{--}\,3$ for $n_{\infty}=10^{5}\,\mathrm{cm^{-3}}$ and $M_{\mathrm{BH}}=10^{2}\,M_{\odot}$. We observe stable D-type flows for the velocity range mentioned above, but the shell becomes unstable or disappears (the I-front becomes R-type) for velocities $v_{\infty}/c_{\infty}>3$, as we will see in the next section.

Figure 4 summarizes the main results found in this paper with regard to the stability of the I-front. The points in the figure show the ratio of the shell thickness to the size of the I-front, $\Delta R_{\mathrm{shell}}/R_{\mathrm{IF}}$, as a function of $v_{\infty}/c_{\infty}$ for a large set of simulations, as shown in the legend. We see that the ratio becomes smaller, i.e., the shell becomes thinner, with increasing $v_{\infty}/c_{\infty}$. The filled symbols refer to simulations in which the shell is stable, while open symbols refer to simulations with unstable shells.

The solid lines in the figure show $\Delta R_{\mathrm{shell}}/R_{\mathrm{IF}}$ from the analytical model described by Eq. (24) in Sec. 2.2, for several values of $T_{\mathrm{II}}$, together with the arrows indicating the values of $v_{\mathrm{R}}$, at which the thickness becomes zero according to the model. For the moment, we focus on the curve for $T_{\mathrm{II}}=6\times 10^{4}\,\mathrm{K}$ because the temperature inside the HII region has approximately this value in the runs with $n_{\infty}=10^{5}\,\mathrm{cm^{-3}}$ and $M_{\mathrm{BH}}=10^{2}\,M_{\odot}$ (see Figure 3). With the numerical factor set to $\alpha=0.5$, the analytical curve for $T_{\mathrm{II}}=4\times 10^{4}\,\mathrm{K}$ shows good agreement with the simulation results for $n_{\infty}=10^{5}\,\mathrm{cm^{-3}}$ and $M_{\mathrm{BH}}=10^{2}\,M_{\odot}$. The agreement is good also for all the other simulations in the figure, justifying the validity of our analytical model, as well as the choice of $\alpha=0.5$. The analytical model predicts that the thickness approaches zero as $v_{\infty}$ approaches $v_{\mathrm{R}}$. In the runs with $n_{\infty}=10^{5}\,\mathrm{cm^{-3}}$ and $M_{\mathrm{BH}}=10^{2}\,M_{\odot}$, however, the shell becomes unstable before the velocity reaches $v_{\mathrm{R}}$, as indicated by the open symbols.

We also investigate the dependence of shell thickness on $n_{\infty}$ and $M_{\mathrm{BH}}$, in addition to the dependence on $v_{\infty}$. We performed runs with various $v_{\infty}$, assuming $(n_{\infty},\,M_{\mathrm{BH}})=(10^{4}\,\mathrm{cm^{-3}},\,10^{2}\,M_{\odot})$, $(10^{3}\,\mathrm{cm^{-3}},\,10^{2}\,M_{\odot})$, and $(10^{3}\,\mathrm{cm^{-3}},\,10^{3}\,M_{\odot})$, and plot $\Delta R_{\mathrm{shell}}/R_{\mathrm{IF}}$ of the stable D-type flows in Figure 4. We see that $\Delta R_{\mathrm{shell}}/R_{\mathrm{IF}}$ becomes smaller when decreasing $n_{\infty}$ or $M_{\mathrm{BH}}$. It appears that $\Delta R_{\mathrm{shell}}/R_{\mathrm{IF}}$ is proportional to the parameter combination $M_{\mathrm{BH}}\,n_{\infty}$.

According to our model, the shell thickness depends on parameters other than $v_{\infty}$ only because of changes of the sound speed inside the HII region, $c_{\mathrm{II}}$. We will show below that $c_{\mathrm{II}}$ depends on $M_{\mathrm{BH}}\,n_{\infty}$, and that this dependence can be attributed to changes in the temperature profile inside the ionized region. In Fig. 5, we plot the upstream temperature profiles along the axis of BH motion in the runs with $v_{\infty}/c_{\infty}=2$ and different $n_{\infty}$ and $M_{\mathrm{BH}}$. We normalize the radius by the size of the I-front to directly compare the temperature profiles. We see in Fig. 5 that the steepness of the temperature rise inside the HII region has significant differences among the runs. For the run with BH mass $M_{\mathrm{BH}}=10^{2}\,M_{\odot}$ and $n_{\infty}=10^{5}\,\mathrm{cm^{-3}}$, the temperature rapidly reaches the almost constant value $5\,\text{--}\,6\times 10^{4}\,\mathrm{K}$ inside the HII region, while the rise in temperature becomes slower as $n_{\infty}$ decreases. Therefore, in the lower-density case with $n_{\infty}\leq 10^{4}\,\mathrm{cm^{-3}}$, we need to define an effective temperature inside the HII region to be used in the analytical model. We set the effective temperature to $T_{\mathrm{II}}=1.8\times 10^{4}\,\mathrm{K}$ and $3\times 10^{4}\,\mathrm{K}$ for the case with $n_{\infty}=10^{3}$ and $10^{4}\,\mathrm{cm^{-3}}$, respectively. These are values measured near the ionization front. In addition, with these values adopted, the analytical curves for $\Delta R_{\mathrm{shell}}/R_{\mathrm{IF}}$ are in good agreement with the numerical results shown in Fig. 4. For the case with $n_{\infty}=10^{3}$, we also perform runs with temperature upstream of the ionization front $T_{\infty}=T_{\mathrm{I}}=30\,\mathrm{K}$, hence assuming that X-ray preheating is negligible (see Appendix. A.2 about X-ray preheating). In these runs, we observe that the temperature rapidly increases from $T_{\mathrm{I}}=30\,\mathrm{K}$ to $T_{\mathrm{shell}}=10^{4}\,\mathrm{K}$ at the shock. The thickness of the bow-shock is slightly larger than in the corresponding runs with $T_{\infty}=10^{4}\,\mathrm{K}$ and is well reproduced by the analytical model with $T_{\mathrm{I}}=30\,\mathrm{K}$ and $T_{\mathrm{shell}}=10^{4}\,\mathrm{K}$, because the $(1-(v_{\infty}^{2}+c_{\mathrm{I}}^{2})/2c_{\mathrm{II}}v_{\infty})$ term in Eq. (24) is larger when $c_{\mathrm{I}}$ is lower. As for the $M_{\mathrm{BH}}$ dependence, the run with $(n_{\infty},\,M_{\mathrm{BH}})=(10^{4}\,\mathrm{cm^{-3}},\,10^{3}\,M_{\odot})$ has almost the same temperature profile as the run with $(n_{\infty},\,M_{\mathrm{BH}})=(10^{5}\,\mathrm{cm^{-3}},\,10^{2}\,M_{\odot})$.

Next, let us explain how the temperature profile depends on $n_{\infty}$ and $M_{\mathrm{BH}}$. The temperature just inside the I-front is roughly determined by the balance between the photo-ionization heating and the Ly$\alpha$ cooling. The gas inside the I-front is almost fully ionized, i.e., $x_{\mathrm{HI}}\ll 1$ and $x_{\mathrm{HII}}=x_{\mathrm{e}}\approx 1$, where $x_{\mathrm{HI}}$, $x_{\mathrm{HII}}$, and $x_{\mathrm{e}}$ are the hydrogen neutral fraction, ionized fraction, and the electron fraction, respectively. The gas is roughly in photo-ionization equilibrium: $k_{\mathrm{ion}}\,n_{\mathrm{H}}\,x_{\mathrm{HI}}=\alpha_{\mathrm{B}}\,n_{\mathrm{H}}^{2}\,x_{\mathrm{HII}}x_{\mathrm{e}}\approx\alpha_{\mathrm{B}}\,n_{\mathrm{H}}^{2}$, where $k_{\mathrm{ion}}$ is the hydrogen photo-ionization rate. Thus, the photo-ionization heating rate per unit volume is $\Gamma=\Delta E\,k_{\mathrm{ion}}\,n_{\mathrm{H}}\,x_{\mathrm{HI}}\approx\Delta E\,\alpha_{\mathrm{B}}\,n_{\mathrm{H}}^{2}$, where $\Delta E=\langle h\nu\rangle_{\mathrm{ion}}-13.6\,\mathrm{eV}$ is the mean energy deposited into the gas per photo-ionization. The Ly$\alpha$ cooling rate can be written as $\Lambda=\tilde{\Lambda}_{\mathrm{Ly\alpha}}\,n_{\mathrm{H}}^{2}\,x_{\mathrm{e}}\,x_{\mathrm{HI}}\approx\tilde{\Lambda}_{\mathrm{Ly\alpha}}\,n_{\mathrm{H}}^{2}\,x_{\mathrm{HI}}$ with a $T$-dependent coefficient $\tilde{\Lambda}_{\mathrm{Ly\alpha}}$. As $\alpha_{\mathrm{B}}$ and $\tilde{\Lambda}_{\mathrm{Ly\alpha}}$ are decreasing and increasing functions of $T$, respectively, the equilibrium temperature (i.e., $\Gamma=\Lambda$) is higher if the neutral fraction behind the I-front, $x_{\mathrm{HI}}$, is lower.

In ionization equilibrium, $x_{\mathrm{HI}}$ is inversely proportional to the ionization parameter, $\xi\equiv F_{\mathrm{ion}}/n_{\mathrm{II}}$, where $F_{\mathrm{ion}}$ is the flux of ionizing radiation. Since the absorption of the ionizing photons inside the HII region is insignificant except for the extreme vicinity of the I-front, we assume optically thin flux in our estimate. Then, approximating the size of I-front with the Strömgren radius $r_{\mathrm{Strm}}$ (Eq. 13), we obtain $\xi$ slightly inside the I-front as

where we ignore the $T_{\mathrm{II}}$ dependence in the last expression. From Eqs. (8), and (11)-(12), the analytical model predicts that $L$ is proportional to $M_{\mathrm{BH}}^{3}\,n_{\infty}^{2}$ in the RIAF regime. Therefore, from Eq. (27), the dependence of $\xi$ on $M_{\mathrm{BH}}$ and $n_{\infty}$ is

This explains why the temperature rise inside the HII region is steeper in the runs with larger $n_{\infty}$ for fixed $M_{\mathrm{BH}}$, as well as why we obtain the same temperature profiles in runs having the same $n_{\infty}\,M_{\mathrm{BH}}$.

As we increase $v_{\infty}$, the dense shell in the D-type flow becomes unstable. As a reference case for the unstable D-type flow in the Region II of Fig. 1, we perform the run with $n_{\infty}=10^{5}\,\mathrm{cm^{-3}}$, $M_{\mathrm{BH}}=10^{2}\,M_{\odot}$, $T_{\infty}=10^{4}\,\mathrm{K}$, and $v_{\infty}/c_{\infty}=4$. We show a snapshot of the flow structure after the instability is fully developed and saturated in Fig. 6. We also provide a corresponding 3D view in Fig. 7. The shell is destroyed by the eruptive expansions of the I-front launched from various places, with dense clumps forming in the gaps of the I-front.

The instability in this range of $v_{\infty}$ was seen in the previous 2D simulations of PR13. However, the way the instability grows is significantly different in 2D versus 3D simulations (see Fig.8 of PR13), as we have also confirmed by performing 2D comparative runs. Although the structures due to the instability are seen in various directions in 3D, they tend to converge to the axis of BH motion in 2D due to the assumed axisymmetry, making the structures more ordered and stronger in 2D than in 3D. Furthermore, we observe that the instability can be artificially seeded near the axis in 2D due to the singularity of the spherical coordinates.

The instability can be understood as the instability of a thin shell between the D-type I-front and the shock, studied in previous literature. This instability was found analytically using linear analysis by Giuliani (1979) and confirmed in 2D (Garcia-Segura & Franco, 1996) and 3D (Whalen & Norman, 2008a, b) simulations of an expanding D-type I-fronts around massive stars. As mentioned in Garcia-Segura & Franco (1996), this instability is caused by a mechanism similar to that of the thin-shell instability of an expanding hot bubble in a cold ambient gas found by Vishniac (1983). In both cases, the thermal pressure of inner hot gas and the ram pressure of external cold gas balanced each other in the unperturbed state, but a force imbalance is introduced in the presence of shell distortions because the force from the thermal pressure is orthogonal to the surface but the ram-pressure force is parallel to the flow, enhancing the front distortion. Using linear perturbation analysis, under the assumption of plane-parallel and infinitely-thin shell, Giuliani (1979) showed that the modes with wavelength $\lambda$ several times greater than the recombination length, $d_{\mathrm{r}}\equiv c_{\mathrm{II}}/n_{\mathrm{II}}\alpha_{\mathrm{B}}$, are unstable.

This result implies that the shell around a moving BH is unstable if the following two conditions are satisfied. Firstly, the perturbed front displacement $\delta r$ needs to be larger than the shell thickness in order for the mode to grow. Since in the linear regime the front displacement is $\delta r\sim A\lambda$ with $A\ll 1$, we get $\lambda>\Delta R_{\mathrm{shell}}/A$. Since the longest wavelengths growing modes have $\lambda<R_{\mathrm{IF}}$, we find that unstable modes exist only in the range $\Delta R_{\mathrm{shell}}/A<\lambda<R_{\mathrm{IF}}$. Therefore, $\Delta R_{\mathrm{IF}}/A\leq R_{\mathrm{IF}}$, that is $\Delta R_{\mathrm{IF}}/R_{\mathrm{IF}}\leq A$, is a necessary condition for the instability. Secondly, the condition $d_{\mathrm{r}}\ll R_{\mathrm{IF}}$ is required, in order for some of the above modes to be unstable. In the case examined in this paper, the second condition is safely satisfied because $d_{\mathrm{r}}\sim 20\,\mathrm{au}\,(n_{\mathrm{II}}/10^{5}\,\mathrm{cm^{-3}})^{-1}\,(T_{\mathrm{II}}/4\times 10^{4}\mathrm{K})^{3/2}$ is generally much smaller than $R_{\mathrm{IF}}$.

To understand the condition for the instability using the simulations, we focus on the relationship between the instability and the shell thickness. In Fig. 4 we plot $\Delta R_{\mathrm{shell}}/R_{\mathrm{IF}}$ for the unstable simulations (open symbols), along with the results for the stable simulations (solid symbols) discussed in Sec. 4.1. Since the shell thickness is not a well-defined quantity for the unstable simulations, we plot $\Delta R_{\mathrm{shell}}/R_{\mathrm{IF}}$ obtained from lower-resolution runs where the instability was suppressed due to lack of resolution. These values should be viewed with caution because the shell thickness is resolved only with one or two cells. In Fig. 4, except for the case with $(n_{\infty},\,M_{\mathrm{BH}},\,T_{\infty})=(10^{3}\,\mathrm{cm^{-3}},\,10^{2}\,M_{\odot},\,10^{4}\,\mathrm{K})$, we see that the shells are unstable if

as demarcated by the dashed line. Using Eq. (25), we can rewrite the condition $\Delta R_{\mathrm{shell}}/R_{\mathrm{IF}}=0.05$ into the form $(v_{\infty}/c_{\mathrm{I}})^{-2}(1-(v_{\infty}/c_{\mathrm{I}})(c_{\mathrm{I}}/v_{\mathrm{R}}))=0.05$, which is a quadratic equation of $v_{\infty}$ for fixed $c_{\mathrm{I}}=8.1\,\mathrm{km/s}\,(T_{\mathrm{I}}/10^{4}\,\mathrm{K})^{1/2}$ and $v_{\mathrm{R}}\approx 2c_{\mathrm{II}}=46\,\mathrm{km/s}\,(T_{\mathrm{II}}/4\times 10^{4}\,\mathrm{K})^{1/2}$. The positive solution of the above equation is

where we have assumed $c_{\mathrm{I}}=c_{\infty}$ and $\Delta_{T}\equiv T_{\mathrm{II}}/T_{\mathrm{I}}$. Therefore the D-type I-front is unstable in the velocity range $v_{\mathrm{crit}}<v_{\infty}<v_{\mathrm{R}}$. As also evident in Fig. 4, the velocity range in which the shell is unstable is rather large for large temperature jump $\Delta_{T}$ across the I-front, but tends to become very narrow and nearly disappear when $\Delta_{T}\rightarrow 1$. As expected, in the run shown in Fig. 6, we observe the formation of an extremely thin shell before the instability grows and destroys the shell, supporting the hypothesis that the instability seen in this run is the thin-shell instability found by Giuliani (1979).

The instability criteria above, however, is a necessary but not a sufficient condition, as exemplified by the run with $n_{\infty}=10^{3}\,\mathrm{cm^{-3}}$, $M_{\mathrm{BH}}=10^{2}\,M_{\odot}$, $T_{\infty}=10^{4}\,\mathrm{K}$, and $v_{\infty}/c_{\infty}=2.5$. In this run with $\Delta_{T}=1.8$ (see Fig. 4), the I-front is stable, although the shell is thinner than the critical value $0.05$ ($\Delta R_{\mathrm{shell}}/R_{\mathrm{IF}}\sim 0.02$). We attribute the stability of this run to the weak temperature increase across the I-front ($\Delta_{T}\sim 1$), also observed in the run with the same $n_{\infty}$ and $M_{\mathrm{BH}}$ but with $v_{\infty}/c_{\infty}=2$ in Fig. 5. In runs with the same $n_{\infty}$ and $M_{\mathrm{BH}}$ but $\Delta_{T}\approx 600$ due to $T_{\infty}=30\,\mathrm{K}$, the I-fronts are unstable if $\Delta R_{\mathrm{shell}}/R_{\mathrm{IF}}\leq 0.05$ (Fig. 4), supporting our conjecture that $\Delta_{T}$ determines the strength of the imbalance between thermal and ram pressure in the presence of shell distortions, hence is an important parameter for the growth of the instability. From the simulation results, we infer that $\Delta_{T}\geq 3$ is a necessary condition for the existence the instability.

The ionization parameter $\xi$ and $\Delta_{T}$ (see Sec. 4.1) are proportional to $L^{1/3}\,n_{\mathrm{II}}^{1/3}$ (Eq. 27), and hence we can re-formulate the condition $\Delta_{T}\geq 3$ as

To derive this condition, we used $L=2\times 10^{37}\,\mathrm{erg/s}$ for the runs with $n_{\infty}=10^{4}\,\mathrm{cm^{-3}}$ (see Sec. 3) and $n_{\mathrm{II}}\approx 2\,n_{\infty}=2\times 10^{4}\,\mathrm{cm^{-3}}$ with $v_{\infty}$ around $v_{\mathrm{R}}\sim 2\,c_{\mathrm{II}}$ (Eq. 7). In the RIAF regime, $\xi$ is also proportional to $n_{\infty}\,M_{\mathrm{BH}}$ (Eq. 28), and thus the above condition is equivalent to

In addition, using Eqs. (7) and (16), we find $d_{\mathrm{r}}/R_{\mathrm{IF}}\propto(M_{\mathrm{BH}}n_{\infty})^{-1}$. Hence, the condition for the existence of unstable modes, $d_{\mathrm{r}}/R_{\mathrm{IF}}<1$, is more strongly satisfied with increasing $n_{\infty}\,M_{\mathrm{BH}}$. Although we have confirmed that $d_{\mathrm{r}}/R_{\mathrm{IF}}<1$ in our simulations with $n_{\infty}\,M_{\mathrm{BH}}\geq 10^{5}\,M_{\odot}\,\mathrm{cm^{-3}}$, the relatively larger $d_{\mathrm{r}}/R_{\mathrm{IF}}$ in the run with $n_{\infty}=10^{3}\,\mathrm{cm^{-3}}$, $M_{\mathrm{BH}}=10^{2}\,M_{\odot}$ and $v_{\infty}/c_{\infty}=2.5$ may contribute to suppressing the instability. Finally, we also note that the opacity to ionizing radiation of shell is also proportional to $n_{\infty}\,M_{\mathrm{BH}}$, as shown in Eq. (26). Thus, the lower opacity and self-shielding of the shell in the cases with smaller $n_{\infty}\,M_{\mathrm{BH}}$ may also have a stabilizing effect. In summary, the parameter combination $n_{\infty}\,M_{\mathrm{BH}}$ controls the nature of the flow and the stability of I-front around moving BHs. As inferred from our simulation results, we regard the condition $n_{\infty}\,M_{\mathrm{BH}}\geq 10^{6}\,M_{\odot}\,\mathrm{cm^{-3}}$ as the second necessary condition for having unstable and fragmenting shells around moving BHs.

For velocities $v_{\infty}\gtrsim v_{\rm R}$, the dense shell disappears and the flow becomes R-type. For a range of velocities around the critical value, however, the transition between the D- and R-type flows is difficult to define because the I-front is not steady due to the instability. As a reference case for the R-type flow in the Region III in Figure 1, we perform the run with $n_{\infty}=10^{5}\,\mathrm{cm^{-3}}$, $M_{\mathrm{BH}}=10^{2}\,M_{\odot}$, $T_{\infty}=10^{4}\,\mathrm{K}$, and $v_{\infty}/c_{\infty}=7$ and shows the snapshot at a time after the instability is fully developed and saturated in Figure 8. Unexpectedly, with respect to the results in PR13, we find that the R-type I-front is also unstable. Similarly to the case of the unstable D-type flow, in the non-linear regime, the instability strongly distorts the I-front and dense regions form at the receding gaps of the I-front. This type of instability was not observed in PR13, probably because the resolution at the I-front was not sufficiently high due to the logarithmic grid in the radial direction. By performing low-resolution simulations, we find that the instability disappears when the resolution near the I-front is below $\Delta x/R_{\mathrm{IF}}=0.05$. We would also like to emphasize again that the growth of instability is essentially a 3D process and cannot be fully captured by 2D simulations.

The unstable nature of weak R-type I-front has been known since the late 1960s. Newman & Axford (1967) studied the stability of the weak R-type I-front using linear analysis for a plane-parallel and infinitely-thin I-front. Newman & Axford (1967) explained the mechanism of instability as follows: the parts of the perturbed I-front getting displaced upstream, the relative flow velocity with respect to the I-front becomes higher; thus the density behind the I-front becomes smaller according to the jump condition for weak R-type I-front (see Eq. 1 with minus sign). Consequently, the absorption of ionizing photons by recombining neutral hydrogen is suppressed, leading to a further displacement upstream of the I-front. Vice versa for the perturbed I-front displaced in the downstream direction. The linear analysis showed that all modes are unstable irrespective of their wavelength $\lambda$, and the growth timescale becomes shorter at smaller wavelengths, until $\lambda\simeq d_{\mathrm{r}}$, below which the growth time scale remains nearly constant. In less idealized conditions, however, we expect that the instability for small-scale modes may be suppressed by effects not included in the linear analysis, such as the finite thickness of the I-front and possible presence of magnetic fields (Axford, 1964). For a mode with fixed $\lambda$, the growth time scale of the instability is nearly constant with the velocity $v_{\mathrm{I}}$, but decreases at most by a factor of two as $v_{\mathrm{I}}$ approaches $v_{\mathrm{R}}$. The growth time scale of the fastest modes is of the order of the recombination time scale $t_{\mathrm{rec}}=1/n_{\mathrm{II}}\alpha_{\mathrm{B}}$.

While Newman & Axford (1967) showed that the weak R-type I-fronts are always unstable, they concluded that the growth of the instability was insignificant because the duration of the expanding R-type phase of an HII region is only about $t_{\mathrm{rec}}$, and thus the instability does not have enough time to grow substantially. In the current case of moving BH, however, weak R-type I-fronts continue to exist for a much longer time, and thus there is enough time for the instability to grow until it becomes non-linear and saturates.

As shown in Fig. 8, the simulations results confirm the linear analysis expectations: in the runs with $n_{\infty}=10^{5}\,\mathrm{cm^{-3}}$, $M_{\mathrm{BH}}=10^{2}\,M_{\odot}$, and $v_{\infty}/c_{\infty}=7$ and $10$, the I-front is unstable. At lower ambient density, in the runs with $n_{\infty}=10^{4}\,\mathrm{cm^{-3}}$, $M_{\mathrm{BH}}=10^{2}\,M_{\odot}$, and $v_{\infty}/c_{\infty}=4\,\text{--}\,10$, we again observe the instability of the R-type flows. The distortion of I-front and the density contrast of the front fragments, however, becomes less significant than in the higher density runs. With an even weaker temperature jump, in the runs with $n_{\infty}=10^{3}\,\mathrm{cm^{-3}}$, $M_{\mathrm{BH}}=10^{2}\,M_{\odot}$, and $v_{\infty}/c_{\infty}=3\,\text{--}\,10$, the instability of the R-type flows becomes insignificant. We further confirm the requirement of large $\Delta_{T}$ for the growth of I-front instabilities, by performing runs with $\Delta_{T}\approx 600$, assuming $T_{\mathrm{I}}=T_{\infty}=30\,\mathrm{K}$ (X-ray preheating artificially suppressed). We find that the R-type I-front is unstable in the run with $v_{\infty}=4\,c_{\mathrm{s}}(T=10^{4}\,\mathrm{K})=70\,c_{\infty}$. The instability, however, does not grow in the run with $v_{\infty}=7\,c_{\mathrm{s}}(T=10^{4}\,\mathrm{K})=130\,c_{\infty}$. This case suggests that the instability may become weak if the Mach number is very large, although this case is unlikely to be realized in nature because we expect that $T_{\mathrm{I}}\sim 10^{4}\,\mathrm{K}$ due to X-ray preheating, especially for R-type flows in which a shell optically thick to X-rays is not present. In the runs with $n_{\infty}=10^{4}\,\mathrm{cm^{-3}}$, $M_{\mathrm{BH}}=10^{3}\,M_{\odot}$, and $v_{\infty}/c_{\infty}=7$ and $10$, the front is strongly unstable similarly as in the reference run, because the temperature profiles are identical in these two cases that have the same $n_{\infty}\,M_{\mathrm{BH}}$.

In the same way as the thin-shell instability of the D-type flow, the steep temperature jump is a necessary condition in order for the instability to grow. As mentioned above, the instability is driven by the variation of the density behind the I-front as the flow velocity varies behind the perturbed front. The density behind the front depends on the flow velocity more weakly as the ratio $\Delta_{T}$ approaches unity (Eq. 1). Therefore, a steeper temperature jump enhances the strength of the instability. The condition for a steep enough temperature jump is to have a large ionization parameter $\xi$ just behind the I-front. In our constant luminosity simulations of the D-type and R-type flows, we assume a constant luminosity $L$ at the value $v_{\infty}=4\,c_{\infty}$, irrespective of the actual value of $v_{\infty}$ in the runs (see Sec. 3). This is a good approximation for unstable D-type flows in the range of $v_{\mathrm{crit}}<v_{\infty}<v_{\mathrm{R}}$, as well as for R-type flows with $v_{\infty}\gtrsim v_{\mathrm{R}}$, because for this velocity range the luminosity does not change significantly as a function of $v_{\infty}$. With this assumption we have found that the I-front is unstable owing to large enough $\xi$ if $n_{\infty}\,M_{\mathrm{BH}}\geq 10^{6}\,M_{\odot}\,\mathrm{cm^{-3}}$.

However, for the R-type flows with large $v_{\infty}$ (i.e., when $v_{\infty}\gg v_{\mathrm{R}}$), the reduction of $L$ from the value assumed in our constant luminosity simulations can be significant. Therefore, here we derive a more accurate instability criterion allowing a realistic dependence of $L$ on $v_{\infty}$. In the high-velocity limit ($v_{\infty}\gg v_{\mathrm{R}}$), Eqs. (8) and (10)-(12), give $L\propto M_{\mathrm{BH}}^{3}\,n_{\infty}^{2}v_{\infty}^{-6}$, and thus Eq. (27) yields $\xi\propto L^{1/3}\,n_{\mathrm{II}}^{1/3}\propto M_{\mathrm{BH}}\,n_{\infty}\,v_{\infty}^{-2}$. From Eqs. (8) – (12), the actual value of $L$ for $v_{\infty}=100\,\mathrm{km/s}$ is about $10^{3}$ times smaller than the assumed constant value of $L$ for $v_{\infty}=4\,c_{\infty}$. Therefore, by rewriting the condition of $\xi$ larger than the critical value corresponding to $n_{\infty}\,M_{\mathrm{BH}}=10^{6}\,M_{\odot}\,\mathrm{cm^{-3}}$ and $L$ for $v_{\infty}=4\,c_{\infty}$ with Eq. (27), we obtain the instability condition for R-type I-fronts: