Impact of Ly $\alpha$ radiation force on super-Eddington accretion onto a massive black hole

In order to account for the multiple scattering of Ly$\alpha$ photons in highly optically thick flows, we solve the frequency-dependent $0$-th moment equation of radiation, in which the radiation diffusion approximation and the Fokker–Planck approximation are employed, with terms up to $O(\varv/c)$ being taken into account. Since the photon diffusion time scale is much shorter than the dynamical time in the present situation, the radiation fields can be regarded to become immediately steady. Applying these approximations and assumptions to Eq.(95.18) (95.19) of Mihalas & Mihalas (1984), the basic equation to be solved here is

	$\displaystyle-\frac{1}{r^{2}}\partialderivative{r}\quantity(\frac{c{r^{2}}}{3\kappa_{x_{0}}}\partialderivative{E_{x_{0}}}{r})$
	$\displaystyle+\frac{1}{r^{2}}\partialderivative{r}\quantity(r^{2}\varv E_{x_{0}})+\frac{1}{r^{2}}\partialderivative{r}\quantity(r^{2}\varv)\frac{1}{3}E_{x_{0}}$
	$\displaystyle-\frac{1}{r^{2}}\partialderivative{r}\quantity(r^{2}\varv)\partialderivative{x_{0}}\quantity[\quantity(\frac{c}{\varv_{\rm th}}+x_{0})\frac{1}{3}E_{x_{0}}]$
	$\displaystyle=-\kappa_{x_{0}}^{\rm a}cE_{x_{0}}+4\pi\eta_{x_{0}}^{\rm r}+\frac{1}{2}\partialderivative{x_{0}}\quantity(c\kappa^{\rm s}_{x_{0}}\partialderivative{E_{x_{0}}}{x_{0}}),$		(1)

where $r$, $c,\varv$ and $\varv_{\rm th}$ are the distance from the origin (BH), the speed of light, the fluid velocity and the thermal velocity in the gas, respectively. The normalized frequency $x_{0}$ is defined as $x_{0}=(\nu_{0}-\nu_{\alpha})/\Delta\nu_{D}$, where $\nu_{0},\nu_{\alpha}$ and $\Delta\nu_{D}$ are the frequency of photons in the fluid comoving frame, the frequency at the Ly$\alpha$ line center, $\nu_{\alpha}=2.466\times 10^{15}{\rm Hz}$, and the typical Doppler shift by the thermal motion of gas around the Ly$\alpha$ line, $\Delta\nu_{D}=\varv_{\rm th}\nu_{\alpha}/c$, respectively. $\kappa^{a}_{x_{0}},\kappa^{s}_{x_{0}},\kappa_{x_{0}}$ and $\eta^{r}_{x_{0}}$ are the absorption coefficient related to the two-photon decay and the collisional de-excitation, the scattering coefficient of the resonance scattering, the total opacity defined by $\kappa_{x_{0}}=\kappa_{x_{0}}^{s}+\kappa_{x_{0}}^{a}$ and the emission coefficient, respectively, all of which are measured in the fluid comoving frame and at normalized frequency of $x_{0}$. The radiation energy density $E_{x_{0}}$ is defined in the comoving frame as $E_{x_{0}}\equiv\int\frac{1}{c}I_{x_{0}}d\Omega_{0}$ where $I_{x_{0}}$ is the intensity in the comoving frame at the frequency $x_{0}$, where $\Omega_{0}$ is the solid angle measured in the comoving frame. Eq.1 includes effects related to the fluid velocity, such as trapping effects by the inflowing medium, the work from/to fluid and the Doppler effects by the velocity dispersion, which correspond to the second, third and fourth terms of the left-hand side, respectively.

The scattering coefficient for the pure resonant scattering, $\bar{\kappa}^{s}_{x_{0}}$, is defined as ${\bar{\kappa}_{x_{0}}^{s}}\equiv\sigma_{0}H(a,x_{0})n_{\rm HI}$ where $\sigma_{0},H(a,x_{0})$ and $n_{\rm HI}$ are the cross-section at the Ly$\alpha$ line center, the Voigt function, and the number density of the neutral hydrogen gas, respectively. The cross-section $\sigma_{0}$ and the Voigt function are defined by

$$\sigma_{0}=f_{12}\frac{\sqrt{\pi}e^{2}}{m_{e}c\Delta\nu_{D}}\approx 5.88\times 10^{-14}\quantity(\frac{T}{10^{4}{\rm K}})^{-1/2}{\rm~{}cm^{2}},$$

(2)

$$H(a,x)=\frac{a}{\pi}\int_{-\infty}^{\infty}\frac{e^{-y^{2}}}{(x-y)^{2}+a^{2}}dy,$$

(3)

where $e,m_{e},f_{12}$ and $a$ are the elementary charge, the electron mass, the oscillator strength, $f_{12}=0.4162$, and the Voigt parameter, $a=\Delta\nu_{\rm L}/(2\Delta\nu_{D})$, which is the ratio between the natural width $\Delta\nu_{\rm L}$ and the thermal width of the Ly$\alpha$ line, respectively.

The absorption coefficient $\kappa_{x_{0}}^{a}$, due to the two-photon decay and the collisional de-excitation, is written as $\kappa^{a}_{x_{0}}=p_{\rm dest}\bar{\kappa}_{x_{0}}^{s}$ where the photon destruction probability $p_{\rm dest}$ is the probability that a Ly$\alpha$ photon is annihilated at a single resonance-scattering. With photon destruction, the scattering coefficient $\kappa^{s}$ becomes slightly different from that for pure resonant scattering $\bar{\kappa}^{s}$ and can be written as $\kappa^{s}_{x_{0}}=(1-p_{\rm dest})\bar{\kappa}^{s}_{x_{0}}$. Assuming that the hydrogen is in the ground state or in the first excited state and that the higher excited states can be ignored, the photon destruction probability can be calculated as

$$p_{\rm dest}=\frac{A_{2\rm ph}+C_{21}}{A_{\rm Ly\alpha}+A_{2\rm ph}+C_{21}}.$$

(4)

where $A_{{\rm Ly}\alpha}$ is the rate of the ${\rm Ly}\alpha$ production by spontaneous emission, $A_{\rm 2ph}$ is the rate of the two photon decay and $C_{21}$ is the rate of the collisional de-excitation from the first excited state to the ground state, respectively. These rates can be written as

	$\displaystyle A_{\rm Ly\alpha}$	$\displaystyle=\frac{n_{2\rm p}}{n_{2}}A_{\rm 2p1s}=\frac{1}{1+r_{\rm 2s2p}}A_{\rm 2p1s},$		(5)
	$\displaystyle A_{\rm 2ph}$	$\displaystyle=\frac{n_{2\rm s}}{n_{2}}A_{\rm 2s1s}=\frac{r_{\rm 2s2p}}{1+r_{\rm 2s2p}}A_{\rm 2s1s}.$		(6)

where $A_{\rm 2p1s}=6.25\times 10^{8}\>{\rm s^{-1}},A_{\rm 2s1s}=8.25\>{\rm s^{-1}}$ and $r_{\rm 2s2p}$ are the Einstein coefficients for the transition ${\rm 2p}\to{\rm 1s}$, ${\rm 2s}\to{\rm 1s}$ and the ratio between the number densities in the 2s and the 2p states, respectively. Due to the very short timescale to achieve the equilibrium among the number of particles in the 1s, 2s and 2p states, the transition rates from/into the 2s state can be considered to balance and therefore the ratio $r_{\rm 2s2p}$ can be calculated as

$$r_{\rm 2s2p}=\frac{C_{\rm 2p2s}}{C_{\rm 2s2p}+A_{\rm 2s1s}}.$$

(7)

where $C_{\rm 2s2p}$ and $C_{\rm 2p2s}$ are the rates of the transition 2s$\to$ 2p and 2p$\to$2s by collisions, respectively. They are written as

$$C_{\rm 2s2p}=5.31\times 10^{-4}\quantity(\frac{n_{\rm e}}{1\mathrm{~{}cm}^{-3}})\mathrm{~{}s}^{-1},$$

(8)

$$C_{\rm 2p2s}=C_{\rm 2s2p}\frac{g_{\rm 2s}}{g_{\rm 2p}}=1.77\times 10^{-4}\quantity(\frac{n_{\rm e}}{1\mathrm{~{}cm}^{-3}})\mathrm{~{}s}^{-1}.$$

(9)

where $g_{\rm 2s}=2$ and $g_{\rm 2p}=6$ are the statistical weight of the 2s state and the 2p state and $n_{e}$ is the electron number density, respectively (see Seaton 1955).

The rate of the collisional de-excitation (Omukai, 2001; Inayoshi et al., 2016), $C_{21}$, is given by

$$C_{21}=\gamma_{21}({\rm e})n_{\rm e}+\gamma_{21}({\rm H})n_{\rm H},$$

(10)

where $n_{\rm H}$ is the hydrogen number density, $\gamma_{21}({\rm e})$ and $\gamma_{21}({\rm H})$ are the de-excitation coefficients for the collisions with electrons and hydrogen atoms, respectively, which are given by substituting $(u,l)=(2,1)$ into the formulae (Sobel’Man et al., 1995; Drawin, 1969) as below:

$$\gamma_{ul}({\rm e})=10^{-8}\left(\frac{l^{2}}{u^{2}-l^{2}}\right)^{3/2}\frac{l^{3}}{u^{2}}\alpha_{lu}\frac{\sqrt{\beta(\beta+1)}}{\beta+\chi_{lu}},\quad\beta=\frac{h\left(v_{l}-v_{u}\right)}{kT},$$

(11)

	$\displaystyle\gamma_{ul}({\rm H})=7.86\times 10^{-15}\left(\frac{l}{u}\right)^{2}\left(\frac{1}{l^{2}}-\frac{1}{u^{2}}\right)^{-2}f_{lu}T^{1/2}$
	$\displaystyle\times\frac{1+1.27\times 10^{-5}\left(1/l^{2}-1/u^{2}\right)^{-1}T}{1+4.76\times 10^{-17}\left(1/l^{2}-1/u^{2}\right)^{-2}T^{2}}.$		(12)

where the coefficients are $\alpha_{12}=24$ and $\chi_{12}=0.28$.

We solve discretized Eq.1 based on the finite volume method by using the alternating-direction implicit (ADI) method. We set the gas temperature, the ionization degree, the mean molecular weight to be $T=10^{4}{\rm K}$, $x_{e}=n_{e}/n_{\rm HI}=10^{-4}$ and $\mu=1$, respectively (see also Inayoshi et al. 2016; Sakurai et al. 2016). The profiles of the gas density and velocity are calculated by solving the fluid equation for the Bondi accretion without radiation effects, using the 11 parameter sets of $(n_{\infty},M_{\rm BH})$ (all parameter sets are shown in Table 1 of the model "Bondi") The Bondi accretion rate, $\dot{M}_{\rm B}$, is given by

$$\dot{M}_{\rm B}\approx 4\times 10^{3}\left(\frac{{L}_{\rm Edd}}{c^{2}}\right)\left(\frac{M_{\rm BH}}{10^{4}M_{\odot}}\right)\left(\frac{n_{\infty}}{10^{5}\>{\rm cm^{3}}}\right)\left(\frac{T}{10^{4}{\rm K}}\right)^{-3/2}.$$

(13)

where $L_{\rm Edd}$ is the Eddington luminosity for the Thomson scattering.

The total optical depth at the Ly$\alpha$ line center, measured over the radial line of the computational domain, $\tau_{0}$, can be approximately written as

$\displaystyle\tau_{0}\approx 4.3\times 10^{11}\left(\frac{n_{\infty}}{10^{5}\>{\rm cm^{3}}}\right)\left(\frac{M_{\rm BH}}{10^{4}M_{\odot}}\right)\left(\frac{T}{10^{4}{\rm K}}\right)^{-3/2}\quantity(\frac{r_{\min}}{10^{-3}r_{\rm B}})^{-1/2},$

(14)

assuming that the number density distirbution $n(r)\propto r^{-3/2}$ and the outer boundary of the radial integral of $\sigma_{0}n(r)$ is set to infinity, where $r_{\rm B}$ is the Bondi radius which depends on $M_{\rm BH}$, $T$, and $r_{\rm min}$ which is the innermost radius of the computational domain defined below.

The radii of the inner and outer boundaries of the computational domain are $r_{\rm min}=10^{-3}r_{\rm B}$ and $r_{\rm max}=10^{-1}r_{\rm B}$, respectively. The innermost radius is the same as in Inayoshi et al. (2016); Sakurai et al. (2016). The number of radial grid cells is set to be $N=800$. The radial grid is set so that the sequence of the radial length of cells becomes the geometrical sequence with the initial value $\Delta r_{0}=10^{-7}r_{\rm min}$ and the common ratio $\epsilon$. Here, $\Delta r_{0}$ corresponds to the size of the finest (innermost) grid cell and the common ratio $\epsilon\approx 1.02131$ satisfies the equation below:

$\displaystyle r_{\rm max}-r_{\rm min}=\Delta r_{0}\frac{\epsilon^{N}-1}{\epsilon-1}.$

(15)

Such a grid design is convenient since we can control the size of the finest grid cell without changing $r_{\rm min},r_{\rm max}$ and $N$. Note that it is irrelevant that the simulation domain does not include the Bondi radius because both the gravitational and radiation forces at such a distant region are small and their impacts on the flow are negligible (see also Fig.1).

The grid for the normalized frequency spans over the range of $[x_{\rm min},x_{\rm max}]$, where we set $x_{\min}=-2(a\tau_{0})^{1/3}$, and $x_{\max}=5(a\tau_{0})^{1/3}$. The grid interval $\Delta x_{0}$ is set to unity within the range $[-10,10]$, and $\Delta x_{0}=(x_{\rm max}-x_{\rm min})/100$ otherwise, considering the significant and steep change in the line profile near the line center.

The situation where Ly$\alpha$ photons are entering the computational domain, passing through the inner boundary at the rate of $L_{{\rm Ly}\alpha}$, is reproduced by setting the emissivity at the innermost cell as follows:

$$\eta_{x_{0}}(r_{0})=\frac{L_{\rm Ly\alpha}}{4\pi}\frac{1}{\Delta V_{0}}\frac{H(x_{0})}{\sqrt{\pi}}.$$

(16)

where $\Delta V_{0}$ is the volume of the innermost cell. The inner boundary condition of the spatial grid is set as the mirror boundary.

Assuming the radiation is isotropic and can leak freely from outer boundary, the outer boundary condition can be described as

$$F_{0}(r_{\max},x_{0})=\frac{1}{2}cE_{0}(r_{\max},x_{0}).$$

(17)

The radiation energy density at the frequency grid boundaries $x_{0}=x_{\min}$ and $x_{\max}$ is set to be zero at all radii.

We investigate the two different cases of accretion flow named as the "slim" and "Eddington" cases. In the "slim" case, the Ly$\alpha$ luminosity is given as

$$L_{\rm Ly\alpha}=0.2\quantity[1+\ln\quantity(\dfrac{{\dot{M}_{\rm B}c^{2}}}{20L_{\rm Edd}})]L_{\rm Edd},$$

(18)

by assuming $10\%$ of the disk luminosity is converted into Ly$\alpha$ luminosity (Sakurai et al., 2016) and using the disk luminosity of the slim disk model from Watarai et al. (2000). It should be noted that the accretion rate of our models is typically $10^{3}$ times larger than the Eddington limit as seen in Eq.13. In the "Eddington" case, we employ $L_{\rm Ly\alpha}=0.1L_{\rm Edd}$ (Inayoshi et al., 2016). We will primarily focus on discussing the result of the "slim" case while the "Eddington" case will also be used to compare with the previous studies.

In Fig.1 we present radial profiles of the Ly$\alpha$ radiation force per unit volume,

$\displaystyle f_{\rm rad}=-\int\frac{1}{3}\partialderivative{E_{x_{0}}}{r}dx_{0},$

(19)

along with gravity for the case with $M_{\rm BH}=10^{5.5}M_{\odot}$ and $n_{\infty}=10^{5}{\rm~{}cm^{-3}}$ both for the "slim" and "Eddington" cases.

It is found from the figure that the Ly$\alpha$ radiation force is higher in the "slim" case than in the "Eddington" case since the Ly$\alpha$ photons are emitted more efficiently in the "slim" case.

In both cases, within the region $r-r_{\min}\lesssim 10^{12}{\rm~{}cm}$, where $r_{\rm min}\simeq 6.2\times 10^{16}\>{\rm cm}$, the radiation force is significantly enhanced due to numerous photon scatterings, exceeding gravity by approximately $10^{6}$ times. However, as the radius increases, the radiation force decreases monotonically and drastically in both cases, falling below gravity for $r-r_{\min}\gtrsim 5\times 10^{16}{\rm~{}cm}$. This decline is attributed to the two-photon decay process. Ly$\alpha$ photons are destroyed in the vicinity of the inner boundary since the photon destruction length, which is the distance a photon travels before it is destroyed via the two photon decay, is very small $\sim 10^{12}{\rm\>cm}$ (see also Eq. 33 in Appendix B). Thus, the number of photons decreases dramatically with increasing $r-r_{\rm min}$, leading to a reduction in the radiation force.

The observed trend, where the radiation force surpasses gravity in the innermost region and decreases with increasing radius, is consistent across all parameter sets. For computational domains with smaller optical depths, the radius (normalized by the Bondi radius) of the region where the radiation force exceeds gravity is larger. In the case with the smallest optical depth in this study $(M_{\rm BH},n_{\infty})=(10^{4}{M_{\odot}},10^{4}{\rm cm^{-3}})$, the radiation force in the "slim" case overwhelms gravity throughout the entire computational domain.

Fig.2 depicts the relationship between the total optical depth at the Ly$\alpha$ line center across the entire computational domain and the force multiplier $M_{\rm F}$¹¹1It is worth noting that the force multiplier is sometimes defined as the ratio between the actual radiation force and the radiation force considering only Thomson opacity (Castor et al., 1975). However, the definition used here is different from it., which is defined as:

$$M_{\rm F}=\frac{c}{L_{\rm Ly\alpha}}\int f_{\rm rad}dV.$$

(20)

We can observe that the force multiplier remains almost independent of the optical depth, maintaining a constant value of $M_{\rm F}\simeq 130$. The reason for this is primarily attributed to the effect of the two-photon decay, rather than the Doppler shift caused by the velocity and its gradient.

To illustrate the reason, we examine the following two cases. The first is the case of a static, homogeneous hydrogen medium where Ly$\alpha$ photons are assumed not to be destroyed. In this case, $M_{\rm F}$ can be analytically described²²2We use the definition ${\rm arctanh}(x)=\frac{1}{2}\ln{\left(\frac{1+x}{1-x}\right)}.$ as (Tomaselli & Ferrara, 2021)

$\displaystyle M_{\rm F}^{\rm ana}(\tau_{0})=\frac{4}{\pi}\sqrt{\frac{2}{3}}\int_{-\infty}^{\infty}{\rm arctanh}(e^{-\pi|\sigma(x)|/\tau_{0}}){dx},$

(21)

where $\sigma(x)=\int_{0}^{x}\sqrt{2/3}/H(y,a)dy$ In the limit of a large optical depth, this equation can be approximated as:

$$M_{\rm F}^{\rm ana}(\tau_{0})\approx 3.51(a\tau_{0})^{1/3}=610\quantity(\frac{\tau_{0}}{10^{10}})^{1/3}\left(\frac{T}{10^{4}{\rm K}}\right)^{-1/6},$$

(22)

which is also represented in Fig.2 by the black dot-dashed line. The second is the case where destruction of Ly$\alpha$ photons are taken into consideration in a static homogeneous hydrogen gas. The results of numerical calculation are shown by the green points. It is observed that the resulting $M_{\rm F}$ is considerably smaller than that of the first case, closely agreeing with the outcomes of our present study. This indicates that the radiation force diminishes and $M_{\rm F}$ becomes constant with $\tau_{0}$, not due to the velocity gradient, but owing to photon destruction. This fact can be understood as follows: Ly$\alpha$ photons undergo numerous scattering and vanish when traveling the distance equivalent to the mean free path for the photon destruction process. Consequently, the Ly$\alpha$ radiation force effectively work only within the extent that Ly$\alpha$ photons can penetrate, and it becomes approximately constant regardless of $\tau_{0}$, even when the optical depth of the system is considerably large. Since the region penetrated by Ly$\alpha$ photons is extremely narrow, the velocity change within the region is small. Therefore, the effect of the velocity field on $M_{\rm F}$ is negligibly small (see Appendix B for more details). Note that the constancy of the force multiplier holds true only when the gas number density $n_{\infty}$ is between $10^{8}{\>\rm cm^{-3}}$ and $10^{12}{\>\rm cm^{-3}}$ as mentioned in Appendix B. It is also noteworthy that the properties of $M_{\rm F}$ discussed in this subsection are independent of $L_{\rm Ly\alpha}$. This is because the radiation force, $f_{\rm rad}$, is proportional to $L_{\rm Ly\alpha}$ via $E_{x_{0}}$ (see Eq. 19).

Fig.3 summarizes the possibility of the supercritical accretion on the $n_{\infty}$-$M_{\rm BH}$ plane for the "slim" models. The filled symbols show the models where the super-Eddington accretion is possible as $F_{\rm grav}>F_{\rm rad}$. Here, $F_{\rm grav}$ and $F_{{\rm rad}}$ are gravitational and radiation forces integrated across the entire simulation domain, defined as follows:

	$\displaystyle F_{\rm grav}$	$\displaystyle=4\pi\int_{r_{\rm min}}^{r_{\rm max}}{G\rho M_{\rm BH}}dr,$		(23)
	$\displaystyle F_{\rm rad}$	$\displaystyle=\int_{r_{\rm min}}^{r_{\rm max}}f_{\rm rad}4\pi r^{2}dr=\frac{M_{F}}{c}L_{{\rm Ly}\alpha}.$		(24)

The radiation force $F_{\rm rad}$ is proportional to the force multiplier $M_{F}$ and the Ly$\alpha$ luminosity. The latter is given by the product of the Eddington luminosity and a factor that is either fixed to $0.1$ for the "Eddington" case or depends on $n_{\infty}M_{\rm BH}$ from Eqs.13 and 18 for the "slim" case.

This figure shows that the Ly$\alpha$ radiation force elevates the lower limit of $n_{\infty}M_{\rm BH}$ for realizing super-Eddington accretion. The grey area represents the parameter range where the super-Eddington accretion is unattainable even in the absence of the Ly$\alpha$ radiation force (Inayoshi et al., 2016; Sakurai et al., 2016)³³3Note that Toyouchi et al. (2019) derives the relation where the required value for $n_{\infty}M_{\rm BH}$ is one order of magnitude larger than the previous studies by considering the pressure balance on the boundary of HII region, but without the Ly$\alpha$ radiation.. By incorporating the Ly$\alpha$ radiation force, the region of non-supercritical accretion expands into the shaded area. Our calculations confirm that the seven models within the shadow region are incapable of supporting super-Eddington accretion.

The boundary of the shaded area and the unshaded region (super-Eddington are) can be derived through analytical considerations. Assuming the density profile of the Bondi accretion at the inner limit, $n(r)=n_{\infty}\left(r/r_{\rm B}\right)^{-3/2}$, the graviational force $F_{\rm grav}$ (Eq.23) can be written as

$\displaystyle{F_{\rm grav}\simeq 4.2\frac{L_{\rm Edd}}{c}\left(\frac{n_{\infty}}{10^{5}{\rm cm^{-3}}}\right)\left(\frac{M_{{\rm BH}}}{10^{4}M_{\odot}}\right)\left(\frac{T}{10^{4}{\rm\>K}}\right)^{-1}\left(\frac{r_{\rm min}}{10^{-3}r_{\rm B}}\right)^{-1/2}}$

(25)

and thus the ratio of $F_{\rm rad}$ and $F_{\rm grav}$ can be expressed as

	$\displaystyle\frac{F_{\rm rad}}{F_{\rm grav}}\sim 2$	$\displaystyle\left(\frac{n_{\infty}}{10^{5}{\rm cm^{-3}}}\right)^{-1}\left(\frac{M_{\rm BH}}{10^{4}M_{\odot}}\right)^{-1}$
	$\displaystyle\times$	$\displaystyle\quantity(\frac{M_{\rm F}}{100})\quantity(\frac{L_{\rm Ly\alpha}}{0.1L_{\rm Edd}})\quantity(\frac{r_{\min}}{10^{-3}r_{\rm B}})^{1/2}\left(\frac{T}{10^{4}{\rm K}}\right).$		(26)

The force multiplier, $M_{\rm F}$, is approximately $100$ in the parameter sets employed in our study (see also fig.2), allowing us to deduce $(n_{\infty}/10^{5}{\rm cm^{-3}})(M_{{\rm BH}}/10^{4}M_{\odot})={40}$ from $F_{\rm rad}=F_{\rm grav}$ for the "slim" case. Consequently, super-Eddington accretion occurs when the gas density and the BH mass satisfy the condition $(n_{\infty}/10^{5}{\rm cm^{-3}})(M_{{\rm BH}}/10^{4}M_{\odot})\gtrsim{40}$ for the "slim" case. A similar argument yields the condition $(n_{\infty}/10^{5}{\rm cm^{-3}})(M_{{\rm BH}}/10^{4}M_{\odot})\gtrsim 2$ for the "Eddington" case. Note that, in the "slim" case, the condition is more stringent because of the luminosity exceeding the Eddington limit, resulting in a stronger radiation force.

The reason why the relation depends only on $M_{BH}n_{\infty}$ even in the "slim" case is attributed to the fact that the ratio $L_{{\rm Ly}\alpha}/L_{\rm Edd}$ depends only on $\dot{M}_{\rm B}/(L_{\rm Edd}/c^{2})$, which is proportional to $n_{\infty}M_{\rm BH}$ (see Eqs.13 and 18). The boundary drawn under this assumption agrees with our results. It should be noted that the analytical description of $F_{\rm rad}/F_{\rm grav}$ includes the inner boundary radius, $r_{\rm min}$, which is fixed to $r_{\rm min}=10^{-3}r_{\rm B}$ in our model but nearly unknown in realistic cases. If the inner boundary radius $r_{\rm min}$, corresponding to the HII region redius, increases, the gravitational force becomes weaker in the innermost region, facilitating the dominance of the radiation force over gravity. According to the simulation of the super-Eddington flow in Inayoshi et al. (2016), the size of the HII region in the early phase is similar to the Bondi radius. At this stage, ${\rm Ly}\alpha$ radiation becomes more important. Under the condition of $r_{\rm min}\sim r_{\rm B}$, the required value of $n_{\infty,5}M_{{\rm BH},4}$ to realize the super-Eddington accretion becomes about thirty times larger than that with $r_{\rm min}=10^{-3}r_{\rm B}$ and reaches $\sim 1500$ for the "slim" case.