Opacity Limit for Supermassive Protostars

In the theory of star formation, it has been a long-standing quest to understand the limits to fragmentation in a given cloud setting. An influential idea was that fragmentation proceeds hierarchically in a collapsing cloud, as the Jeans mass decreases with increasing density as long as the cloud can collapse almost isothermally (Hoyle 1953). The minimum fragment mass is then set by the scale when opacity prevents the release of the gravitational energy via radiation, such that the Jeans mass would increase again upon further compression (Low & Lynden-Bell 1976; Rees 1976). Fragmentation can be seeded in a number of ways, including from nonspherical perturbations of the Larson–Penston solution (Hanawa & Matsumoto 2000; Lai 2000). We here follow a similar reasoning, applied to the peculiar conditions of isothermally collapsing primordial gas in atomic-cooling halos. Our goal is to robustly derive the characteristic density ${n}_{{\rm{F}}}$, mass ${M}_{{\rm{F}}}$, and radius ${R}_{{\rm{F}}}$ of the emerging protostar, when the gas first becomes optically thick (see third panel of Figure 1). These values will mark the initial stage in the build-up process of the supermassive object.

We start by considering the simple relation between these three quantities:

$$\begin{eqnarray}&&{M}_{{\rm{F}}}=\displaystyle \frac{4\pi }{3}\displaystyle \frac{{m}_{{\rm{H}}}}{X}{n}_{{\rm{F}}}{R}_{{\rm{F}}}^{3},\end{eqnarray} \tag{ 1 }$$

where we have used $\rho ={m}_{{\rm{H}}}n/X$, with X = 0.76 being the primordial hydrogen mass fraction, to translate total mass density to hydrogen number density.

We furthermore assume that the optically thick cloud is gravitationally bound, such that the characteristic mass is of the order of the Jeans mass ${M}_{{\rm{J}}}$,

$$\begin{eqnarray}&&{M}_{{\rm{F}}}\simeq {M}_{{\rm{J}}}={\left(\displaystyle \frac{\pi {k}_{{\rm{B}}}}{{{Gm}}_{{\rm{H}}}\mu }\right)}^{3/2}{T}^{3/2}{n}_{{\rm{F}}}^{-1/2},\end{eqnarray} \tag{ 2 }$$

where $\mu \simeq 1.22$ is the mean molecular weight for a fully neutral primordial mixture of atomic hydrogen and helium.

Finally, we need to account for energy equilibrium. The energy that is to be radiated away originates in the gravitational collapse of the cloud. In this case, the gravitational energy is emitted in a collapse timescale ${t}_{\mathrm{col}}$, as long as the gas remains optically thin to its cooling radiation. Right before the gas cloud becomes opaque, the energy is radiated from the surface as a fraction ${f}_{\mathrm{BB}}$ of the blackbody radiation. Hence, we can equate the gravitational compressional heating rate with the radiation cooling rate

$$\begin{eqnarray}&&\displaystyle \frac{4\pi }{3}{R}_{{\rm{F}}}^{3}\cdot \displaystyle \frac{{n}_{{\rm{F}}}{k}_{{\rm{B}}}T}{{t}_{\mathrm{col}}}={f}_{\mathrm{BB}}4\pi {R}_{{\rm{F}}}^{2}{\sigma }_{\mathrm{SB}}{T}^{4},\end{eqnarray} \tag{ 3 }$$

where ${\sigma }_{\mathrm{SB}}$ is the Stefan–Boltzmann constant. In a one-zone model, where the thermal evolution of the central core of a collapsing cloud is calculated, the collapse timescale is commonly assumed to be ${t}_{\mathrm{col}}=\sqrt{3\pi /32G\rho }$, which is the time for the density of an initially static cloud to reach infinity. On the other hand, the dynamical timescale in the free-fall collapse is ${t}_{\mathrm{col}}=1/\sqrt{24G\rho }$, which is shorter by a factor of $3\pi /2\simeq 4.7$. In reality, the collapse timescale can be between these values. Here, we set the collapse timescale to ${t}_{\mathrm{col}}={f}_{\mathrm{col}}\sqrt{3\pi /32G\rho }$ in order to consider this uncertainty, where $0.2\lesssim {f}_{\mathrm{col}}\lesssim 1$.

We then proceed to solve the system of Equations (1)–(3), and obtain analytic expressions of the characteristic density, radius, and mass

$$\begin{eqnarray}&&{n}_{{\rm{F}}}\simeq 4.6\times {10}^{16}\,f{\left(\displaystyle \frac{T}{3000{\rm{K}}}\right)}^{5/2}\,{\mathrm{cm}}^{-3},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{R}_{{\rm{F}}}\simeq 0.33\,{f}^{-1}{\left(\displaystyle \frac{T}{3000{\rm{K}}}\right)}^{-3/4}\,\mathrm{au}\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{M}_{{\rm{F}}}\simeq 0.045\,{f}^{-1/2}{\left(\displaystyle \frac{T}{3000{\rm{K}}}\right)}^{1/4}\,{M}_{\odot },\end{eqnarray} \tag{ 6 }$$

where we have normalized T to the typical value of isothermally collapsing gas in the high density regime in an atomic-cooling halo, and we have used $f={f}_{\mathrm{col}}{f}_{\mathrm{BB}}\lesssim 1$. We note that this argument for the universal line in the density–temperature plane (Equation (4)), on which the gas cloud becomes optically thick to any continuum opacities (gas and dust grains) has already been discussed by Omukai et al. (2005), and we here reproduce the key result in a simplified way to highlight the basic physics involved.

Figure 2 shows the density–temperature diagram of a collapsing cloud in an atomic-cooling halo when a protostar, composed of a hydrostatic and adiabatic core, forms at the center of the cloud (red dashed curve). The data is taken from a three-dimensional (3D) hydrodynamical simulation by Inayoshi et al. (2014), including all relevant cooling and chemical reaction networks. In this case, the gas becomes opaque at $n\simeq 5\times {10}^{15}\,{\mathrm{cm}}^{-3}$ (open circle), which agrees with ${n}_{{\rm{F}}}=9.2\times {10}^{15}\,{\mathrm{cm}}^{-3}$ for f = 0.2 within a factor of two (see Equation (4)). Moreover, the radius and mass of the opaque core estimated from Equations (5) and (6) are ${R}_{F}=1.7\,\mathrm{au}$ and ${M}_{F}=0.1\,{M}_{\odot }$ for f = 0.2, respectively. Both values also reasonably agree with the 3D simulation results, where ${R}_{F}\simeq 1\,\mathrm{au}$ and ${M}_{\star }\simeq 0.2\,{M}_{\odot }$ right after protostellar formation.

Zoom In Zoom Out Reset image size

Next, we gain further insight by considering the detailed physics of ${{\rm{H}}}^{-}$ opacity, based on the actual microphysical cross section for this process. To this extent, we rederive the opacity limit by using a one-zone model for the collapse of gas into an atomic-cooling halo.

Following previous works (Omukai 2001; Inayoshi et al. 2014), we implement the cooling function due to H⁻ free–bound (fb) and ff emission

$$\begin{eqnarray}&&{\rm{H}}+{e}^{-}\to {{\rm{H}}}^{-}+\gamma ,\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{\rm{H}}+{e}^{-}\to {\rm{H}}+{e}^{-}+\gamma ,\end{eqnarray} \tag{ 8 }$$

and consider three opacity sources associated with H⁻ bound–free, ff transition and H Rayleigh scattering

$$\begin{eqnarray}&&{{\rm{H}}}^{-}+\gamma \to {\rm{H}}+{e}^{-},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{\rm{H}}+{e}^{-}+\gamma \to {\rm{H}}+{e}^{-},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{\rm{H}}+\gamma \to {\rm{H}}+\gamma ^{\prime} .\end{eqnarray} \tag{ 11 }$$

These processes are treated in a self-consistent way with chemical reaction networks. An updated set of chemical reaction rate coefficients and cross sections is summarized in Inayoshi et al. (2014, 2016).

In what follows, we briefly describe the method introduced by Inayoshi et al. (2014) to calculate the cooling function both in the optically thin and thick regimes. We summarize specific functional forms for the cooling rates and opacities in the Appendix. In the optically thin limit, the ${{\rm{H}}}^{-}$ cooling rate is estimated by integrating emissivities over frequency as

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{\mathrm{thin}}\equiv 4\pi \int ({\eta }_{\nu }^{\mathrm{fb}}+{\eta }_{\nu }^{\mathrm{ff}})d\nu .\end{eqnarray} \tag{ 12 }$$

We here divide the frequency range into two: ${D}_{l}=[0,0.75]$ eV and ${D}_{h}=[0.75,13.6]$ eV, called "low" and "high" frequency, respectively. That is,

$$\begin{eqnarray}{{\rm{\Lambda }}}_{\mathrm{thin}} & = & {\displaystyle \int }_{{D}_{l}}+{\displaystyle \int }_{{D}_{h}}4\pi ({\eta }_{\nu }^{\mathrm{fb}}+{\eta }_{\nu }^{\mathrm{ff}})d\nu ,\\ & \equiv & {{\rm{\Lambda }}}_{\mathrm{thin}}^{(l)}+{{\rm{\Lambda }}}_{\mathrm{thin}}^{(h)}.\end{eqnarray} \tag{ 13 }$$

This distinction between the two ranges is required to calculate the cooling rate in the optically thick case. Figure 3 shows the optical depth at the core of the collapsing cloud due to the H⁻ bound–free/ff transition and H Rayleigh scattering for three different densities. For the lowest density ($n={10}^{12}\,{\mathrm{cm}}^{-3}$, top panel), the gas is optically thin (${\tau }_{\nu }=n{\sigma }_{\nu }{\lambda }_{{\rm{J}}}\lt 1$) to all the continuum opacities at frequencies $\lesssim 2\,\mathrm{eV}$. As the density increases to $n={10}^{15}\,{\mathrm{cm}}^{-3}$ (middle panel), the optical depth at higher frequencies ($\gtrsim 1$ eV) exceeds unity, but the H⁻ ff emission still works as radiation cooling. For the highest density ($n={10}^{16}\,{\mathrm{cm}}^{-3}$, bottom panel), the gas core becomes completely opaque to all the continuum, and hence enters the opacity limit.

Zoom In Zoom Out Reset image size

In the optically thick limit, the cooling function is approximated as

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{\mathrm{thick}}\simeq {\int }_{{D}_{l}}+{\int }_{{D}_{h}}\displaystyle \frac{-4\pi }{3({\kappa }_{\nu }^{{\rm{a}}}+{\kappa }_{\nu }^{{\rm{s}}})}\displaystyle \frac{{\partial }^{2}{B}_{\nu }(T)}{\partial {z}^{2}}d\nu ,\end{eqnarray} \tag{ 14 }$$

where ${\kappa }_{\nu }^{{\rm{a}}({\rm{s}})}$ is the absorption (scattering) coefficients, ${B}_{\nu }(T)$ is the Planck function, and z is the coordinate along the temperature gradient. Here, we approximate Equation (14) as

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{\mathrm{thick}}\simeq \displaystyle \sum _{i=l,h}\displaystyle \frac{4\pi {\int }_{{D}_{i}}{B}_{\nu }(T)d\nu }{3{\kappa }_{R}^{(i)}{{\ell }}^{2}}=\displaystyle \sum _{i=l,h}\displaystyle \frac{4\pi {\int }_{{D}_{i}}{\eta }_{\nu }d\nu }{3{\kappa }_{R}^{(i)}{\kappa }_{P}^{(i)}{{\ell }}^{2}},\end{eqnarray} \tag{ 15 }$$

where the partial derivative $\partial /\partial z$ is replaced with a characteristic length ℓ, and ${\kappa }_{R}^{l(h)}$ and ${\kappa }_{P}^{l(h)}$ are Rosseland and Planck mean opacities in the low- and high-frequency regimes, respectively. Note that, in this limit, the emissivity is expressed as ${\eta }_{\nu }={\kappa }_{\nu }^{{\rm{a}}}{B}_{\nu }(T)$ because the source function is given by ${B}_{\nu }(T)$.

Finally, in order to connect both the optically thin and thick regimes, we adopt the following functional form

$$\begin{eqnarray}{{\rm{\Lambda }}}_{\mathrm{tot}} & \simeq & \displaystyle \sum _{i=l,h}\displaystyle \frac{4\pi {\displaystyle \int }_{{D}_{i}}{\eta }_{\nu }d\nu }{1+3{\kappa }_{R}^{(i)}{\kappa }_{P}^{(i)}{{\ell }}^{2}},\\ & = & \displaystyle \sum _{i=l,h}\displaystyle \frac{{{\rm{\Lambda }}}_{\mathrm{thin}}^{(i)}}{1+3{\kappa }_{R}^{(i)}{\kappa }_{P}^{(i)}{{\ell }}^{2}}\equiv {{\rm{\Lambda }}}^{(l)}+{{\rm{\Lambda }}}^{(h)},\end{eqnarray} \tag{ 16 }$$

where the first (second) term in the right-hand side in Equation (16) mainly corresponds to H⁻ fb (ff) emission. Figure 4 shows the evolution of the cooling rates for a collapsing cloud in our one-zone calculation (see solid curve in Figure 2). Each solid curve presents the total cooling rate (red), the rate of ${{\rm{\Lambda }}}^{(h)}$ (green) and ${{\rm{\Lambda }}}^{(l)}$ (blue). The H⁻ fb cooling saturates and decreases at $n\gt 6\times {10}^{15}\,{\mathrm{cm}}^{-3}$ because of H Rayleigh scattering and H⁻ bound–free absorption. At $9\times {10}^{15}\,{\mathrm{cm}}^{-3}\lesssim n\lesssim 2\times {10}^{16}\,{\mathrm{cm}}^{-3}$, the H⁻ ff emission acts as the main cooling process instead of H⁻ fb. Since the compressional heating, given by ${{\rm{\Gamma }}}_{\mathrm{comp}}={{nk}}_{{\rm{B}}}T/{t}_{\mathrm{ff}}$ (dashed curve in Figure 4), dominates the total cooling rate during this transition, the temperature begins to increase gradually. Eventually, the gas becomes completely opaque at $n\gt 2\,\times {10}^{16}\,{\mathrm{cm}}^{-3}$, where $T\propto {n}^{2/3}$. Note that this density, here derived by considering the detailed microphysics involved, is very similar to the estimate in Section 2.1. We can thus robustly characterize the conditions at the onset of supermassive protostar formation.

Zoom In Zoom Out Reset image size

After the collapse and formation of the optically thick object, its mass grows through accretion of the surrounding gas and new sources of energy start becoming important (see fourth panel of Figure 1). In that case, its evolution will not be determined by the energy radiated away from the collapse, but by the interplay between internal and accretion radiation from the protostar. During the first stage of protostellar evolution, the energy powering the object will dominantly come from accretion rather than self-gravitating collapse. At some point toward the later evolution of the system, the accretion timescale, ${t}_{\mathrm{acc}}={M}_{\star }/{\dot{M}}_{\star }$, becomes larger than the Kelvin–Helmholtz (KH) timescale, ${t}_{\mathrm{KH}}={{GM}}_{\star }^{2}/{R}_{\star }{L}_{\star }$, and hence the protostellar model needs to be augmented by internal contributions (e.g., Omukai & Palla 2003).

Right after a protostar forms, the accretion timescale is shorter than the KH timescale. In this accretion phase (${t}_{\mathrm{acc}}\lesssim {t}_{\mathrm{KH}}$), we modify the left-hand side of Equation (3) and consider the accretion luminosity released at the stellar surface:

$$\begin{eqnarray}&&\displaystyle \frac{{{GM}}_{\star }{\dot{M}}_{\star }}{{R}_{\star }}\simeq 4\pi {R}_{\mathrm{ph}}^{2}{\sigma }_{\mathrm{SB}}{T}_{\star }^{4},\end{eqnarray} \tag{ 17 }$$

where ${R}_{\mathrm{ph}}$ is the photospheric radius. In the early stages of the accretion phase, we assume ${R}_{\mathrm{ph}}\simeq 1.4\,{R}_{\star }$ (Stahler et al. 1986), which is derived for a spherically symmetric, quasi-steady model of an accreting protostar. Specifically, a freely falling envelope is depositing material onto a growing hydrostatic core in an accretion shock. The latter is surrounded by a radiative precursor, transitioning into the optically thin envelope at the photosphere. Here, the factor of 1.4 is determined by ${{\rm{H}}}^{-}$ opacity, but during later evolutionary stages, other opacity effects like electron scattering will play a role. In accretion problems, other radii like the trapping radius, defined as the point where the radiative diffusion and dynamical timescales are of the same order, might be important. During the evolution of the system, the effective opacity can be written as a function of the opacity due to absorption (${\tau }_{\mathrm{abs}}$) and to electron scattering (${\tau }_{\mathrm{es}}$) as ${\tau }_{\mathrm{eff}}=\sqrt{{\tau }_{\mathrm{abs}}({\tau }_{\mathrm{abs}}+{\tau }_{\mathrm{es}})}$. Inside the photosphere, the gas is ionized and hence electron scattering dominates over absorption (${\tau }_{\mathrm{es}}\gg {\tau }_{\mathrm{abs}}$). We can then estimate the trapping radius due to electron scattering as ${R}_{\mathrm{tr}}^{\mathrm{es}}={\dot{M}}_{\star }{\sigma }_{{\rm{T}}}/4\,\pi {m}_{{\rm{H}}}c\simeq 4.45\,\mathrm{au}\left(\tfrac{{\dot{M}}_{\star }}{{M}_{\odot }\,{\mathrm{yr}}^{-1}}\right)$ (Begelman 1978). For our assumed accretion rate of ${\dot{M}}_{\star }\simeq 1\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$, the trapping radius due to electron scattering is of the order of a few astronomical units, which is consistent with the approximation ${R}_{\mathrm{tr}}^{\mathrm{es}}\lesssim {R}_{\mathrm{ph}}\sim {R}_{\star }$. On the other hand, we can also estimate the trapping radius due to ${{\rm{H}}}^{-}$ absorption, given by equating the diffusion timescale and the free-fall time as ${\tau }_{{{\rm{H}}}^{-}}R/c\simeq \sqrt{3\pi /32G\rho }$. To evaluate this expression, we take the results from the one-zone model and solve for the radius. This gives us a value of ${R}_{\mathrm{tr}}^{{{\rm{H}}}^{-}}\simeq 4\times {10}^{-2}\,\mathrm{au}$, which lies inside the photospheric radius at all times. As a result, we can adequately assume that the trapping radius might not influence the evolution of the object throughout the protostellar assembly. Future, self-consistent radiation-hydrodynamical simulations will provide a more complete understanding. With these assumptions, we then evaluate the stellar radius as a function of stellar mass ${M}_{\star }$, accretion rate ${\dot{M}}_{\star }$, and surface temperature ${T}_{\star }$:

$$\begin{eqnarray}&&{R}_{\star }\simeq 1.1{\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right)}^{1/3}{\left(\displaystyle \frac{{\dot{M}}_{\star }}{{M}_{\odot }{\mathrm{yr}}^{-1}}\right)}^{1/3}{\left(\displaystyle \frac{{T}_{\star }}{6000{\rm{K}}}\right)}^{-4/3}\,\mathrm{au}.\end{eqnarray} \tag{ 18 }$$

Later in the evolution of the object, the internal luminosity from the star begins to dominate because the opacity decreases as the temperature increases inside the star ($\kappa \propto {T}^{-7/2}$), resulting in ${t}_{\mathrm{acc}}\gt {t}_{\mathrm{KH}}$. In a typical case for Population III star formation with a moderate accretion rate of ${\dot{M}}_{\star }\,\simeq {10}^{-3}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$, the star contracts losing the thermal energy via radiative diffusion and forms a main-sequence star. However, when the accretion rate is sufficiently high, the total luminosity (i.e., the accretion luminosity and the internal luminosity) tends to exceed the Eddington luminosity during the KH contraction. Then, the stellar surface expands in order to regulate the increase of the total luminosity, and the protostar evolves into a red-giant-like structure with a contracting core and an expanding envelope. The critical accretion rate to bifurcate the protostellar evolution is estimated as ${\dot{M}}_{\star }\gtrsim 4\,\times {10}^{-3}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ (Omukai & Palla 2003). According to stellar evolution calculations, even if ${t}_{\mathrm{acc}}\gt {t}_{\mathrm{KH}}$, the stellar surface continues to expand without contraction phases when the accretion rate is higher than ${\dot{M}}_{\star }\gtrsim {10}^{-2}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ (Hosokawa et al. 2012). In this case, the stellar luminosity approaches the Eddington value for the corresponding mass. Hence, the energy equilibrium equation can be written as

$$\begin{eqnarray}&&{L}_{\mathrm{Edd},\star }=\displaystyle \frac{4\pi {{GM}}_{\star }{m}_{{\rm{H}}}c}{{\sigma }_{{\rm{T}}}}\simeq 4\pi {R}_{\star }^{2}{\sigma }_{\mathrm{SB}}{T}_{\star }^{4},\end{eqnarray} \tag{ 19 }$$

from where we can derive an expression for the stellar radius during the expansion phase, as a function of the mass of the star and the surface temperature (Hosokawa et al. 2012, 2013):

$$\begin{eqnarray}&&{R}_{\star }\simeq 0.78{\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right)}^{1/2}{\left(\displaystyle \frac{{T}_{\star }}{6000{\rm{K}}}\right)}^{-2}\,\mathrm{au}.\end{eqnarray} \tag{ 20 }$$

In Figure 5, we compare the stellar radius (blue solid) for ${T}_{\star }\simeq 6000\,{\rm{K}}$ to other important length scales in our problem. Among them, we consider the Schwarzschild radius, ${R}_{\mathrm{Sch}}=2{{GM}}_{\star }/{c}^{2}$ (red solid), the innermost stable circular orbit (ISCO) radius, ${R}_{\mathrm{ISCO}}=6{{GM}}_{\star }/{c}^{2}$ (green solid), and the Bondi radius, ${R}_{{\rm{B}}}={{GM}}_{\star }/{c}_{{\rm{s}}}^{2}$ (purple solid). The protostar's radius grows from ${R}_{\star }\simeq 0.65\,\mathrm{au}$ at ${M}_{\star }\simeq 0.1{M}_{\odot }$ to ${R}_{\star }\simeq 250\,\mathrm{au}$ at ${M}_{\star }\simeq {10}^{5}{M}_{\odot }$, well below the Bondi radius, but above both the Schwarzschild and ISCO radii for the whole range of masses. Additionally, the massive protostar does not enter the region where it becomes GR unstable (pink shaded region) during its evolution. This might indicate that the whole star does not collapse due to GR instability, but only its core, in agreement with Hosokawa et al. (2013).

Zoom In Zoom Out Reset image size

For the Bondi radius, we here assume for simplicity a nearly constant sound speed, corresponding to $T(r)\simeq {T}_{\star }\simeq 6000\,{\rm{K}}$ (black dashed line in Figure 6). This argument is based on the Larson–Penston collapse (Larson 1969; Penston 1969), which approximately describes the evolution of atomic-cooling halos. In such a case, the density follows a profile $\rho \propto {r}^{-2}$ and hence the evolution of the temperature is nearly isothermal, up to radii of a few parsecs. Furthermore, we have verified that our simulations of the initial stages of collapse (see Section 3.2) exhibit such near-isothermality out to $\simeq 100$ pc, as shown by the red solid line in Figure 6. It is clear that we oversimplify the situation here. In reality, the infalling matter will heat up when transitioning to optically thick conditions inside the Bondi radius, and eventually follow an adiabatic evolution. In such a case, $\gamma =5/3$ and the temperature has a radial dependency $T(r)\propto {r}^{-1}$ (Shapiro & Teukolsky 1983), as represented for a stellar mass ${M}_{\star }\simeq {10}^{4}\,{M}_{\odot }$ by the black dotted line in the same figure. In general, the transonic radius for Bondi accretion is given by ${R}_{{\rm{s}}}=\left(\tfrac{5-3\gamma }{4}\right)\tfrac{{GM}}{{c}_{{\rm{s}},\infty }^{2}}$ (Shapiro & Teukolsky 1983). If the gas evolution is characterized by a different γ value (e.g., $\gamma =1.1$ for classical Population III star formation, Omukai & Nishi 1998), the Bondi radius might change by a factor of a few. Additionally, the transition to the adiabatic stage depends on how the diffusion and free-fall timescales compare to each other. Since the values for the trapping and photospheric radii are similar, this implies that gas outside the photospheric radius is not affected by this increase in temperature, further validating our assumption.

Zoom In Zoom Out Reset image size

One factor that might dramatically influence the evolution described in this section is the radiation emitted from the central accreting protostar. Here, we have assumed a constant surface temperature because of strong temperature dependence of H⁻ opacity. Throughout the evolution of the central object, however, its temperature will vary and eventually reach a point where radiative feedback becomes important. Previous studies have analyzed the formation of primordial supermassive stars in the rapid mass accretion regime and have found that the effective temperature of the object remains well below ${10}^{4}\,{\rm{K}}$, for protostellar masses up to ${10}^{4}{M}_{\odot }$, or so (Hosokawa et al. 2012, 2013), suggesting that radiative feedback might not become important up to those mass scales.

It is useful to explore the likely temperature evolution of the growing supermassive protostar, in response to a realistic mass accretion history provided by a cosmological simulation. To this extent, we consider the photospheric temperature, given by the general stellar evolution calculations of Stahler et al. (1986):

$$\begin{eqnarray}&&{T}_{\star }\simeq 4000\,{\rm{K}}{\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right)}^{0.044}{\left(\displaystyle \frac{{\dot{M}}_{\star }}{{M}_{\odot }{\mathrm{yr}}^{-1}}\right)}^{-0.055}.\end{eqnarray} \tag{ 21 }$$

We have plotted this relation in Figure 7 for accretion rates of 10⁻² (red), 10⁻¹ (blue), and $1\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ (green). Because this relation is only valid when ${t}_{\mathrm{acc}}\lesssim {t}_{\mathrm{KH}}$, we have used solid lines up to the mass where this inequality inverts. For higher masses, we make a rough extrapolation based on the same expression (dotted lines), although the evolution of the temperature for this stage is unclear. As can be seen, the temperature at which the photosphere begins to emit H-ionizing radiation, ${T}_{\mathrm{ion}}\simeq$ 10,000 K (black dotted), is not reached in the range of accretion rates explored here. For lower values of ${\dot{M}}_{\star }$ the stellar radius is smaller than the photospheric one, and follows the zero-age main-sequence evolution. In such a case, the ionizing temperature can be reached well before ${M}_{\star }\simeq {10}^{4}\,{M}_{\odot }$.

Zoom In Zoom Out Reset image size

In addition to our constant accretion rate assumption, we have included a time-dependent toy model of the form ${\dot{M}}_{\star }(t)=1\,{M}_{\odot }\,{\mathrm{yr}}^{-1}{e}^{-t/{t}_{\mathrm{ff},0}}$ (purple dashed–dotted line), with ${t}_{\mathrm{ff},0}={10}^{5}\,\mathrm{year}$ being the free-fall time in the core of the atomic-cooling halo before the collapse (e.g., Safranek-Shrader et al. 2012) and ${M}_{\star }(t=0)={M}_{{\rm{F}}}\simeq 0.1\,{M}_{\odot }$. We also compare with the Hosokawa et al. (2013) results for the effective temperature ${T}_{\mathrm{eff}}$ at ${\dot{M}}_{\star }\simeq 1\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ (orange dashed line). Both models agree well, but in the former case the increase in the effective temperature is due to a drop in the accretion rate, while in the latter case this is a result of the decrease in the opacity because of the expansion of the stellar radius. Similar to the case of constant ${\dot{M}}_{\star }$, none of these models reach ${T}_{\mathrm{ion}}$ during the evolution of the protostar up to masses of ${10}^{5}\,{M}_{\odot }$. Since the photospheric temperature at the characteristic accretion rate ${\dot{M}}_{\star }\simeq 1\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ never surpasses 10⁴ K, we can safely neglect photoionization feedback from the central protostar during most of its evolution. This radiation only becomes important for accretion rates a few orders of magnitude lower than $1\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$, or for masses ${M}_{\star }\simeq {10}^{5}\,{M}_{\odot }$ in the case of our time-dependent model.

Lastly, the final mass of the protostar can be affected by continuum radiation-driven mass loss once its total luminosity exceeds the Eddington luminosity (e.g., Fiacconi & Rossi 2016), or by mass loss due to pulsations (Inayoshi et al. 2013). As previously discussed, in the early stages of the evolution, the total luminosity remains below the Eddington value and hence it is not affected by radiation-driven mass loss. However, as the mass of the protostar grows its luminosity increases and this scenario changes. In such a case, the final mass of the protostar can vary significantly, as studied by Fiacconi & Rossi (2016). Furthermore, the accreting supermassive protostar might become pulsationally unstable, but the estimated mass-loss rates are too low to effectively prevent protostellar growth (Inayoshi et al. 2013). In summary, mass loss, either due to continuum radiation or pulsations, should not affect the early evolutionary stages, but continuum opacity might become important later on, when the protostellar mass approaches $\simeq {10}^{5}-{10}^{6}\,{M}_{\odot }$ (see the fifth panel of Figure 1).

In the case of numerical simulations, the evolution of the central protostar requires either the implementation of sink particles (e.g., Latif et al. 2013b) or an artificially stiffened equation of state (e.g., Hirano & Bromm 2016). For the former, we can use the treatment in the previous section of the protostellar evolution to construct a physically motivated subgrid model.

The formation of the central object in the optically thick regime is characterized by a central fragment of density ${n}_{{\rm{F}}}$ and radius ${R}_{{\rm{F}}}$, and a isothermal profile of the form $n\propto {r}^{-2}$ outside that scale, as represented by the dashed line in the density–radius diagram in Figure 8. The Bondi radius of the object is given by

$$\begin{eqnarray}&&{R}_{{\rm{B}}}\simeq 17\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right){\left(\displaystyle \frac{T}{6000{\rm{K}}}\right)}^{-1}\,\mathrm{au}.\end{eqnarray} \tag{ 22 }$$

As seen from Figure 8, ${R}_{{\rm{B}}}$ is larger than ${R}_{{\rm{F}}}$, and hence the relation between both quantities follow the isothermal profile, from which we can derive:

$$\begin{eqnarray}{n}_{{\rm{B}}} & = & {n}_{{\rm{F}}}{\left(\displaystyle \frac{{R}_{{\rm{F}}}}{{R}_{{\rm{B}}}}\right)}^{2}\\ & \simeq & 4.66\times {10}^{13}{f}^{-2}{\left(\displaystyle \frac{T}{6000{\rm{K}}}\right)}^{3}{\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right)}^{-2}\,{\mathrm{cm}}^{-3},\end{eqnarray} \tag{ 23 }$$

where ${n}_{{\rm{B}}}=n(r={R}_{{\rm{B}}})$ is defined as the density at the Bondi radius.

Zoom In Zoom Out Reset image size

The solid line in Figure 8 corresponds to the subsequent evolution, which is characterized by the growth of the stellar radius following Equation (20) for the stage where ${t}_{\mathrm{acc}}\gtrsim {t}_{\mathrm{KH}}$ and accretion rates ${\dot{M}}_{\star }\gtrsim {10}^{-2}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$. Outside the Bondi radius the evolution is still described by the isothermal profile $n\propto {r}^{-2}$, but inside ${R}_{{\rm{B}}}$ material falls within a free-fall time, and hence it is represented by the relation $n\propto {r}^{-3/2}$ in the density–radius space. This allows us to relate the stellar and the Bondi radius as

$$\begin{eqnarray}&&{n}_{\star }={n}_{{\rm{B}}}{\left(\displaystyle \frac{{R}_{{\rm{B}}}}{{R}_{\star }}\right)}^{3/2}\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&\simeq \,5.13\times {10}^{15}{f}^{-1}{\left(\displaystyle \frac{T}{6000{\rm{K}}}\right)}^{9/2}{\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right)}^{5/4}\,{\mathrm{cm}}^{-3},\end{eqnarray} \tag{ 25 }$$

where ${n}_{\star }=n(r={R}_{\star })$ is defined as the density at the stellar radius.

Furthermore, inside the Bondi radius, we can estimate the accretion rate of infalling gas as ${\dot{M}}_{{\rm{B}}}=4\pi {n}_{{\rm{B}}}{m}_{{\rm{H}}}{c}_{{\rm{s}}}{R}_{{\rm{B}}}^{2}$. Using Equation (23) for ${n}_{{\rm{B}}}$ and Equation (22) for ${R}_{{\rm{B}}}$, we derive

$$\begin{eqnarray}&&{\dot{M}}_{{\rm{B}}}\simeq 0.78{f}^{-2}{\left(\displaystyle \frac{T}{6000{\rm{K}}}\right)}^{3/2}{M}_{\odot }\,{\mathrm{yr}}^{-1},\end{eqnarray} \tag{ 26 }$$

The framework depicted here can be used to derive a prescription for the accretion radius of the sink particle, ${R}_{\mathrm{acc}}$. In the ideal case of a simulation with high enough resolution, we would set the accretion radius to ${R}_{{\rm{F}}}$ or ${R}_{\star }$. Unfortunately, this case is not always achievable and hence we need to choose ${R}_{\mathrm{acc}}\lesssim {R}_{{\rm{B}}}$ if the simulations resolve densities $n\gtrsim {n}_{{\rm{B}}}$, with an accretion rate estimated by Equation (26). On the other hand, the scenario becomes more complex when the simulation is only able to resolve a certain threshold density ${n}_{\mathrm{th}}\lt {n}_{{\rm{B}}}$. In such a case, the choice of accretion radius should enclose the Bondi radius and, as suggested in Figure 8, follow the isothermal profile. Hence, the value of the accretion rate for a given threshold density can be expressed as ${R}_{\mathrm{acc}}\,={R}_{{\rm{F}}}{({n}_{{\rm{F}}}/{n}_{\mathrm{th}})}^{1/2}$.

In summary, we can write a formula for the accretion radius based on the threshold density as

$$\begin{eqnarray}{R}_{\mathrm{acc}}\simeq \left\{\begin{array}{ll}17\left(\tfrac{{M}_{\star }}{{M}_{\odot }}\right){\left(\tfrac{T}{6000{\rm{K}}}\right)}^{-1}\,\mathrm{au}, & \mathrm{if}\,{n}_{\mathrm{th}}\gtrsim {n}_{{\rm{B}}}\\ {\left(\tfrac{1.49\times {10}^{17}{\mathrm{cm}}^{-3}}{{n}_{\mathrm{th}}}\right)}^{1/2}{f}^{-1/2} & {\left(\tfrac{T}{6000{\rm{K}}}\right)}^{1/2}\,\mathrm{au},\\ & \mathrm{if}\,{n}_{\mathrm{th}}\lt {n}_{{\rm{B}}}.\end{array}\right.\end{eqnarray} \tag{ 27 }$$

For threshold densities of ${n}_{\mathrm{th}}={10}^{8}\,{\mathrm{cm}}^{-3}$ and ${n}_{\mathrm{th}}\,={10}^{10}\,{\mathrm{cm}}^{-3}$, the initial values of the accretion radius in the low-mass regime (where the Bondi radius is not resolved) are given by ${R}_{\mathrm{acc}}\simeq 3.8\times {10}^{4}\,\mathrm{au}$ and ${R}_{\mathrm{acc}}\simeq 3.8\times {10}^{3}\,\mathrm{au}$, respectively. These values are kept until ${n}_{{\rm{B}}}={n}_{\mathrm{th}}$, which occurs when the mass of the protostar is $\simeq {10}^{3}\,{M}_{\odot }$ for the former and $\simeq {10}^{2}\,{M}_{\odot }$ for the latter. From then on, the accretion radius is given by the Bondi radius and the final mass is determined by the Bondi accretion rate, which is independent of the initial resolution of the simulation. Hence, the final mass does not depend on the choice for the threshold density.

The two initial accretion radii for ${n}_{\mathrm{th}}={10}^{8}\,{\mathrm{cm}}^{-3}$ and ${n}_{\mathrm{th}}={10}^{10}\,{\mathrm{cm}}^{-3}$ are shown in Figure 5 as a black dashed line and a black dotted line respectively. In both cases, the accretion radius is larger than the stellar radius (calculated at $T\,=6000\,{\rm{K}}$) in the mass range $0.1\mbox{--}{10}^{5}\,{M}_{\odot }$ once the sink is formed. As the star evolves, the Bondi radius increases with the mass of the star, eventually reaching the point when it is resolved, in which case, the accretion radius should transition to ${R}_{\mathrm{acc}}={R}_{{\rm{B}}}$. This method ensures that the central star, modeled as a sink particle, will always be enclosed by the accretion radius during its evolution, from its formation until it becomes a supermassive star.

Previous works have used different strategies to implement sink particles. For example, Latif et al. (2013b) and Shlosman et al. (2016) assumed ${R}_{\mathrm{acc}}={{GM}}_{\star }/({c}_{{\rm{s}}}^{2}+{v}_{\infty }^{2})$ and then estimated the accretion rate as ${\dot{M}}_{\star }=4\pi {\rho }_{\infty }{R}_{\mathrm{acc}}\sqrt{1.2544{c}_{{\rm{s}}}^{2}+{v}_{\infty }^{2}}$, while Regan & Downes (2017) used a fixed accretion radius of four cells in the maximum refinement level. We expect all of these recipes to give a similar accretion rate of $\simeq 1{M}_{\odot }\,{\mathrm{yr}}^{-1}$, which corresponds to our estimation from Equation (26).

We perform a low-resolution simulation of the collapse of an atomic-cooling halo following a similar approach to the one described in Becerra et al. (2015). We start from cosmological initial conditions at redshift z = 99 and box size of 2 Mpc (comoving) in a Λ cold dark matter (ΛCDM) cosmology. We then follow the evolution of the halo until the highest density cell reaches the threshold density ${n}_{\mathrm{th}}={10}^{8}\,{\mathrm{cm}}^{-3}$. For that, we have used a primordial chemistry network that includes five species (H, ${{\rm{H}}}_{2}$, ${{\rm{H}}}^{-}$, ${{\rm{H}}}^{+}$, and ${e}^{-}$) and cooling processes such as ${{\rm{H}}}^{-}$ cooling, ${{\rm{H}}}_{2}$ line cooling, ${{\rm{H}}}_{2}$ collision-induced emission, Ly-α cooling, and inverse Compton cooling. The refinement criteria ensures that the Jeans length is resolved by at least 64 cells at every stage of the evolution.

We show the number density (top) and temperature (bottom) projections of the central object at scales of 200 (left) and 5 pc (right) at that instant in time in Figure 9. At large scales, the cloud shows an irregular morphology, but it becomes nearly spherical at scales of $\simeq 10\,\mathrm{pc}$. This implies that the object reaches spherical symmetry at scales larger than the accretion radius at that point (${R}_{\mathrm{acc}}\simeq 0.2\,\mathrm{pc}$, as calculated in Section 3.1), which is plotted in black dashed lines.

Zoom In Zoom Out Reset image size

Finally, we analyze the accretion rate onto the object at the moment when the simulation reaches the threshold density in Figure 10. The radial profile of the accretion rate is shown as a red solid line, which is calculated as $\dot{M}=-4\pi {r}^{2}\rho {v}_{\mathrm{rad}}$, where r is the distance to the highest density cell, ρ is the mass density of hydrogen, and ${v}_{\mathrm{rad}}$ is the radial component of the velocity. For reference, we have also included the Shu accretion rate for spherical collapse, ${\dot{M}}_{\mathrm{Shu}}\simeq 0.975{c}_{{\rm{s}}}^{3}/G$ (Shu 1977), and the Larson–Penston accretion rate for dynamical collapse, ${\dot{M}}_{\mathrm{LP}}\simeq 46.9{c}_{{\rm{s}}}^{3}/G$ (Larson 1969; Penston 1969), as blue and green dotted lines, respectively. The mass accretion reaches a maximum of $\dot{M}\simeq 1.4\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ at $r\simeq 0.5\,\mathrm{pc}$ and then it decreases to values $\dot{M}\simeq 0.1\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ at larger scales, consistent with our estimation from Equation (26). Up to scales of $\simeq 15$ pc, its value lies in between the Shu and the Larson–Penston accretion rates, which remain roughly constant at ${\dot{M}}_{\mathrm{Shu}}\simeq 0.18\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ and ${\dot{M}}_{\mathrm{LP}}\simeq 8.5\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ for the whole radial range. At the accretion radius (shown as a vertical black dashed line), the value of the mass infall rate is $\dot{M}\simeq 0.6\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$, which is consistent with the values assumed throughout this study.

Zoom In Zoom Out Reset image size

Throughout the paper, we discuss a subgrid model within a sink assuming spherical symmetric accretion flows. However, gas material with angular momentum forms an accretion disk, through which the central protostar is fed. In supermassive star formation, the disk would be unstable against its self-gravity because of high accretion rates from the parent cloud (${\dot{M}}_{\star }\sim 1\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$). In such an unstable disk, the disk is likely to fragment into multiple clumps, which could migrate inward losing their orbital angular momentum due to gravitational interaction with the disk and other clumps. Angular momentum redistribution induced by the clumps can drive the evolution of the disk and predict the formation of SMBHs in its nuclei as described by Lodato & Natarajan (2006). Eventually, most of the clumps can feed the gas into the central protostar episodically before the clumps evolve to main-sequence stars, which could suppress the gas accretion through the disk due to ionizing radiation (Inayoshi & Haiman 2014; Latif & Schleicher 2015). Moreover, the radius of the central protostar monotonically increases at ${M}_{\star }\gt {10}^{2}\mbox{--}{10}^{3}\,{M}_{\odot }$ with an almost constant effective temperature of ${T}_{\mathrm{eff}}\simeq 5000\,{\rm{K}}$, resulting in weak radiation feedback. Since the average accretion rate through the disk is as high as $\sim 0.1\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ and the duration of clump accretion episodes is shorter than the KH timescale at the stellar surface (Sakurai et al. 2016), the evolution of the stellar structure is not affected by details of episodic accretion (Sakurai et al. 2015).