Opacity limit for supermassive protostars

TL;DR

This paper models the evolution of supermassive protostars from initial formation to billion-solar-mass scales, identifying the opacity limit and providing a sub-grid recipe for simulations, with implications for feedback processes.

Contribution

It introduces a detailed model of supermassive protostar evolution, including the opacity limit and a new sub-grid prescription for numerical simulations.

Findings

Opacity limit density n_F ≈ 2×10^{17} cm^{-3}

Protostar radius growth from 0.65 AU to 250 AU

Photo-ionization feedback negligible until late stages

Abstract

We present a model for the evolution of supermassive protostars from their formation at $M_\star \simeq 0.1\,\text{M}_\odot$ until their growth to $M_\star \simeq 10^5\,\text{M}_\odot$. To calculate the initial properties of the object in the optically thick regime we follow two approaches: based on idealized thermodynamic considerations, and on a more detailed one-zone model. Both methods derive a similar value of $n_{\rm F} \simeq 2 \times 10^{17} \,\text{cm}^{-3}$ for the density of the object when opacity becomes important, i.e. the opacity limit. The subsequent evolution of the growing protostar is determined by the accretion of gas onto the object and can be described by a mass-radius relation of the form $R_\star \propto M_\star^{1/3}$ during the early stages, and of the form $R_\star \propto M_\star^{1/2}$ when internal luminosity becomes important. For the case of a supermassive protostar, this implies that the radius of the star grows from $R_\star \simeq 0.65 \,{\rm AU}$ to $R_\star \simeq 250 \,{\rm AU}$ during its evolution. Finally, we use this model to construct a sub-grid recipe for accreting sink particles in numerical simulations. A prime ingredient thereof is a physically motivated prescription for the accretion radius and the effective temperature of the growing protostar embedded inside it. From the latter, we can conclude that photo-ionization feedback can be neglected until very late in the assembly process of the supermassive object.

Figures

Equations

$$M_{\mathrm{F}}=\frac{4\pi }{3}\frac{m_{\mathrm{H}}}{X}{n}_{\mathrm{F}}{R}_{\mathrm{F}}^3\text{mbox},$$

$$M_{\mathrm{F}}\simeq M_{\mathrm{J}}={\left(\frac{\pi k_{\mathrm{B}}}{G{m}_{\mathrm{H}}\mu }\right)}^{3/2}{T}^{3/2}{n}_{\mathrm{F}}^{-1/2}\text{mbox},$$

$$\frac{4\pi }{3}{R}_{\mathrm{F}}^3\cdot \frac{n_{\mathrm{F}}{k}_{\mathrm{B}}T}{t_{\mathrm{col}}}=f_{\mathrm{BB}}4\pi R_{\mathrm{F}}^2{\sigma }_{\mathrm{SB}}{T}^4\text{mbox},$$

$$n_{\mathrm{F}}\simeq 4.6\times 10^{16}f{\left(\frac{T}{3000\mathrm{K}}\right)}^{5/2}{\text{cm}}^{-3},$$

$$R_{\mathrm{F}}\simeq 0.33f^{-1}{\left(\frac{T}{3000\mathrm{K}}\right)}^{-3/4}\mathrm{AU}$$

$$M_{\mathrm{F}}\simeq 0.045f^{-1/2}{\left(\frac{T}{3000\mathrm{K}}\right)}^{1/4}{\text{M}}_{\odot }\text{mbox},$$

$$\mathrm{H}+{\mathrm{e}}^{-}\to {\mathrm{H}}^{-}+\gamma ,$$

$$\mathrm{H}+{\mathrm{e}}^{-}\to \mathrm{H}+{\mathrm{e}}^{-}+\gamma ,$$

$${\mathrm{H}}^{-}+\gamma \to \mathrm{H}+{\mathrm{e}}^{-},$$

$$\mathrm{H}+{\mathrm{e}}^{-}+\gamma \to \mathrm{H}+{\mathrm{e}}^{-},$$

$$\mathrm{H}+\gamma \to \mathrm{H}+{\gamma }^{\prime }.$$

$${\Lambda }_{\mathrm{thin}}\equiv 4\pi \int ({\eta }_{\nu }^{\mathrm{fb}}+{\eta }_{\nu }^{\mathrm{ff}})d\nu .$$

$${\Lambda }_{\mathrm{thin}}$$

$$\equiv {\Lambda }_{\mathrm{thin}}^{(l)}+{\Lambda }_{\mathrm{thin}}^{(h)}.$$

$${\Lambda }_{\mathrm{thick}}\simeq {\int }_{D_l}+{\int }_{D_h}\frac{-4\pi }{3({\kappa }_{\nu }^{\mathrm{a}}+{\kappa }_{\nu }^{\mathrm{s}})}\frac{{\partial }^2{B}_{\nu }(T)}{\partial z^2}d\nu ,$$

$${\Lambda }_{\mathrm{thick}}\simeq \sum _{i=l,h}\frac{4\pi {\int }_{D_i}{B}_{\nu }(T)d\nu }{3{\kappa }_R^{(i)}{\ell }^2}=\sum _{i=l,h}\frac{4\pi {\int }_{D_i}{\eta }_{\nu }d\nu }{3{\kappa }_R^{(i)}{\kappa }_P^{(i)}{\ell }^2},$$

$${\Lambda }_{\mathrm{tot}}$$

$$=\sum _{i=l,h}\frac{{\Lambda }_{\mathrm{thin}}^{(i)}}{1+3{\kappa }_R^{(i)}{\kappa }_P^{(i)}{\ell }^2}\equiv {\Lambda }^{(l)}+{\Lambda }^{(h)},$$

$$\frac{G{M}_{\star }{\dot{M}}_{\star }}{R_{\star }}\simeq 4\pi R_{\mathrm{ph}}^2{\sigma }_{\mathrm{SB}}{T}_{\star }^4,$$

$$R_{\star }\simeq$$

$$L_{\mathrm{Edd},\star }=\frac{4\pi G{M}_{\star }{m}_{\mathrm{H}}c}{{\sigma }_{\mathrm{T}}}\simeq 4\pi R_{\star }^2{\sigma }_{\mathrm{SB}}{T}_{\star }^4,$$

$$R_{\star }\simeq 0.78{\left(\frac{M_{\star }}{{\text{M}}_{\odot }}\right)}^{1/2}{\left(\frac{T_{\star }}{6000\mathrm{K}}\right)}^{-2}\mathrm{AU}.$$

$$T_{\star }\simeq 4000\mathrm{K}{\left(\frac{M_{\star }}{{\mathrm{M}}_{\odot }}\right)}^{0.044}{\left(\frac{{\dot{M}}_{\star }}{{\mathrm{M}}_{\odot }{\mathrm{yr}}^{-1}}\right)}^{-0.055}.$$

$$R_{\mathrm{B}}\simeq 17\left(\frac{M_{\star }}{{\text{M}}_{\odot }}\right){\left(\frac{T}{6000\mathrm{K}}\right)}^{-1}\mathrm{AU}.$$

$$n_{\mathrm{B}}$$

$$\simeq 4.66\times 10^{13}{f}^{-2}{\left(\frac{T}{6000\mathrm{K}}\right)}^3{\left(\frac{M_{\star }}{{\text{M}}_{\odot }}\right)}^{-2}{\text{cm}}^{-3}\text{mbox},$$

$$n_{\star }$$

$$\simeq 5.13\times 10^{15}{f}^{-1}{\left(\frac{T}{6000\mathrm{K}}\right)}^{9/2}{\left(\frac{M_{\star }}{{\text{M}}_{\odot }}\right)}^{5/4}{\text{cm}}^{-3}\text{mbox},$$

$${\dot{M}}_{\mathrm{B}}\simeq 0.78f^{-2}{\left(\frac{T}{6000\mathrm{K}}\right)}^{3/2}{\text{M}}_{\odot }{\text{yr}}^{-1}\text{mbox},$$

$$R_{\mathrm{acc}}\simeq \{\begin{array}{l}\text{vspace}1em17\left(\frac{M_{\star }}{{\text{M}}_{\odot }}\right){\left(\frac{T}{6000\mathrm{K}}\right)}^{-1}\mathrm{AU},\text{if}{n}_{\mathrm{th}}\gtrsim n_{\mathrm{B}}\\ {\left(\frac{1.49\times 10^{17}{\text{cm}}^{-3}}{n_{\mathrm{th}}}\right)}^{1/2}{f}^{-1/2}{\left(\frac{T}{6000\mathrm{K}}\right)}^{1/2}\mathrm{AU},\\ \text{if}{n}_{\mathrm{th}}<n_{\mathrm{B}}\end{array}$$

$$k^{(l)}$$

$$k^{(h)}$$

$$\times \left(1.4924+0.07815T_3+0.0063T_3^2\right)$$

Full text

hydrodynamics – stars: formation – galaxies: formation – galaxies: high-redshift– cosmology: theory – early universe

1 Introduction

Recent observations at redshifts $z\gtrsim 6$ suggest that quasars were already powered by supermassive black holes (SMBHs) with masses $\gtrsim 10^{9}\,\text{M}_{\odot}$ when the Universe was less than one billion years old (Fan et al., 2003, 2006; Mortlock et al., 2011; Wu et al., 2015). These SMBHs most likely grew from smaller seed BHs that formed earlier, but the origin of these seeds remains unclear (Haiman, 2006, 2009; Bromm & Yoshida, 2011; Greene, 2012; Volonteri, 2012; Volonteri & Bellovary, 2012). Furthermore, feedback and self-regulation of the seeds make the study of their formation and growth even more complex (Milosavljević et al., 2009). The two most promising theories concerning the formation of seed BHs at high redshift are the remnants of massive Population III stars (Madau & Rees, 2001; Li et al., 2007; Johnson et al., 2012), and the direct collapse of primordial gas in haloes with virial temperatures $T_{\rm vir}\gtrsim 10^{4}\,$K, the so-called atomic cooling haloes (Bromm & Loeb, 2003; Begelman et al., 2006; Spaans & Silk, 2006).

In the direct collapse scenario, high temperatures are reached in haloes where cooling by molecular hydrogen and metal lines to below $\simeq 10^{4}\,$K has been suppressed, which implies that the only coolant acting on the gas is atomic hydrogen (Omukai, 2001; Oh & Haiman, 2002). In the case of molecular hydrogen, which naturally forms at the center of the halo, its photo-dissociation can be achieved by an external soft ultraviolet (UV) background in the Lyman-Werner (LW) bands. Previous studies have found that this leads to a nearly isothermal collapse at $T_{\rm vir}\simeq 10^{4}\,{\rm K}$ due to initially Lyman-$\alpha$ cooling, and subsequently ${\rm H}^{-}$ bound-free and free-free emission, when higher densities are reached (Regan & Haehnelt, 2009; Latif et al., 2013a; Inayoshi et al., 2014; Becerra et al., 2015; Chon et al., 2016). High-resolution simulations have shown that, as the gas collapses and reaches densities of $\simeq 10^{17}\text{cm}^{-3}$, it becomes optically thick to ${\rm H}^{-}$ radiation, and a massive protostar with accretion rate $\simeq 1\,\text{M}_{\odot}\,\text{yr}^{-1}$ forms at the center of the halo (Inayoshi et al., 2014; Van Borm et al., 2014; Becerra et al., 2015; Latif et al., 2016). Due to this high accretion rate, the central object can easily become a supermassive star of $\simeq 10^{5}-10^{6}\,\text{M}_{\odot}$ within a million years (Regan & Haehnelt, 2009; Latif et al., 2013a), which later might collapse into a SMBH due to relativistic instabilities (Baumgarte & Shapiro, 1999; Umeda et al., 2016; Woods et al., 2017, see also Figure 1).

In this work, we study the physics of the central object when it approaches the optically thick regime. In particular, we investigate the properties of the emerging protostar when the optical depth due to ${\rm H}^{-}$ emission becomes unity, thus extending the classical theory of opacity-limited fragmentation developed for present-day star formation (Rees, 1976; Low & Lynden-Bell, 1976). Previous studies have explored this scenario using detailed one-zone models (e.g. Omukai, 2001). Here, we present an alternative approach based on both simplified dimensional arguments and a fitting formula for the cooling and heating processes within the non-equilibrium chemistry of H and ${\rm H}^{-}$ ions. In addition, we develop an idealized model for the subsequent evolution of the accreting protostar, until the formation of a supermassive object. Based on this modeling of the growing protostar, we deduce parameters for a physically-motivated sink particle algorithm, to be used as a sub-grid recipe in large-scale, hydrodynamic simulations of the formation of SMBH seeds in a fully cosmological context. Such simulations are needed to derive detailed diagnostics for the SMBH assembly process at high redshifts, to be probed with next-generation observational facilities (Pacucci et al., 2015), such as the James Webb Space Telescope (JWST), the ATHENA X-ray mission, and the Laser Interferometer Space Antenna (LISA) gravitational-wave observatory.

2 Physics of the opacity limit

2.1 Classical picture

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 (Rees, 1976; Low & Lynden-Bell, 1976). Fragmentation can be seeded in a number of ways, including from non-spherical 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 haloes. 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:

[TABLE]

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 order the Jeans mass $M_{\rm J}$,

[TABLE]

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_{\rm 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_{\rm BB}$ of the black-body radiation. Hence, we can equate the gravitational compressional heating rate with the radiation cooling rate

[TABLE]

where $\sigma_{\rm 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_{\rm 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_{\rm 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_{\rm col}=f_{\rm col}\sqrt{3\pi/32G\rho}$ in order to consider this uncertainty, where $0.2\lesssim f_{\rm col}\lesssim 1$.

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

[TABLE]

[TABLE]

[TABLE]

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_{\rm col}f_{\rm 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}~{}\text{cm}^{-3}$ (open circle), which agrees with $n_{\rm F}=9.2\times 10^{15}~{}\text{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$ AU and $M_{F}=0.1~{}\text{M}_{\odot}$ for $f=0.2$, respectively. Both values also reasonably agree with the 3D simulation results, where $R_{F}\simeq 1$ AU and $M_{\star}\simeq 0.2~{}\text{M}_{\odot}$ right after protostellar formation.

2.2 Detailed modeling: One-zone model

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 re-derive 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 free-free (ff) emission

[TABLE]

[TABLE]

and consider three opacity sources associated with H- bound-free, free-free transition and H Rayleigh scattering

[TABLE]

[TABLE]

[TABLE]

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 regime. 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

[TABLE]

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,

[TABLE]

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 H- bound-free/free-free transition and H Rayleigh scattering for three different densities. For the lowest density ($n=10^{12}~{}\text{cm}^{-3}$, top panel), the gas is optically thin ($\tau_{\nu}=n\sigma_{\nu}\lambda_{\rm J}<1$) to all the continuum opacities at frequencies $\lesssim 2$ eV. As the density increases to $n=10^{15}~{}\text{cm}^{-3}$ (middle panel), the optical depth at higher frequencies ($\gtrsim 1$ eV) exceeds unity, but the H- free-free emission still works as radiation cooling. For the highest density ($n=10^{16}~{}\text{cm}^{-3}$, bottom panel), the gas core becomes completely opaque to all the continuum, and hence enters the opacity limit.

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

[TABLE]

where $\kappa_{\nu}^{\rm a(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

[TABLE]

where the partial derivative $\partial/\partial z$ is replaced with a characteristic length $\ell$, and $\kappa_{R}^{l(h)}$ and $\kappa_{P}^{l(h)}$ are Rosseland and Planck mean opacity in the low and high frequency regime, 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 regime, we adopt the following functional form

[TABLE]

where the first (second) term in the right-hand-side in Equation (16) mainly corresponds to H- free-bound (free-free) 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 $\Lambda^{(h)}$ (green) and $\Lambda^{(l)}$ (blue). The H- free-bound cooling saturates and decreases at $n>6\times 10^{15}~{}\text{cm}^{-3}$ because of H Rayleigh scattering and H- bound-free absorption. At $9\times 10^{15}~{}\text{cm}^{-3}\lesssim n\lesssim 2\times 10^{16}~{}\text{cm}^{-3}$, the H- free-free emission acts as the main cooling process instead of H- free-bound. Since the compressional heating, given by $\Gamma_{\rm comp}=nk_{\rm B}T/t_{\rm 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>2\times 10^{16}~{}\text{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.

2.3 Protostellar evolution

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_{\rm acc}=M_{\star}/\dot{M}_{\star}$, becomes larger than the Kelvin-Helmholtz (KH) timescale, $t_{\rm 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_{\rm acc}\lesssim t_{\rm KH}$), we modify the left-hand-side of Equation (3) and consider the accretion luminosity released at the stellar surface:

[TABLE]

where $R_{\rm ph}$ is the photospheric radius. In the early stages of the accretion phase, we assume $R_{\rm 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_{\rm abs}$) and to electron scattering ($\tau_{\rm es}$) as $\tau_{\rm eff}=\sqrt{\tau_{\rm abs}(\tau_{\rm abs}+\tau_{\rm es})}$. Inside the photosphere the gas is ionized and hence electron scattering dominates over absorption ($\tau_{\rm es}\gg\tau_{\rm abs}$). We can then estimate the trapping radius due to electron scattering as $R_{\rm tr}^{\rm es}=\dot{M}_{\star}\sigma_{\rm T}/4\pi m_{\rm H}c\simeq 4.45\,{\rm AU}\left(\dot{M}_{\star}\over\text{M}_{\odot}/{\rm yr}\right)$ (Begelman, 1978). For our assumed accretion rate of $\dot{M}_{\star}\simeq 1\,\text{M}_{\odot}\,\text{yr}^{-1}$, the trapping radius due to electron scattering is of the order of a few AU, which is consistent with the approximation $R_{\rm tr}^{\rm es}\lesssim R_{\rm 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_{\rm tr}^{{\rm H}^{-}}\simeq 4\times 10^{-2}\,{\rm 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}$:

[TABLE]

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_{\rm acc}>t_{\rm KH}$. In a typical case for Population III star formation with a moderate accretion rate of $\dot{M}_{\star}\simeq 10^{-3}~{}\text{M}_{\odot}\,\text{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}~{}\text{M}_{\odot}\,\text{yr}^{-1}$ (Omukai & Palla, 2003). According to stellar evolution calculations, even if $t_{\rm acc}>t_{\rm KH}$, the stellar surface continues to expand without contraction phases when the accretion rate is higher than $\dot{M}_{\star}\gtrsim 10^{-2}~{}\text{M}_{\odot}\,\text{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

[TABLE]

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 cthe surface temperature (Hosokawa et al., 2012, 2013):

[TABLE]

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_{\rm Sch}=2GM_{\star}/c^{2}$ (red solid), the innermost stable circular orbit (ISCO) radius, $R_{\rm ISCO}=6GM_{\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\,{\rm AU}$ at $M_{\star}\simeq 0.1\text{M}_{\odot}$ to $R_{\star}\simeq 250\,{\rm AU}$ at $M_{\star}\simeq 10^{5}\text{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).

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 pc. 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}\,\text{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(\frac{5-3\gamma}{4}\right)\frac{GM}{c_{{\rm s,}\infty}^{2}}$ (Shapiro & Teukolsky, 1983). If the gas evolution is characterized by a different $\gamma$ 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.

2.4 Onset of radiative feedback

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}\text{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):

[TABLE]

We have plotted this relation in Figure 7 for accretion rates of $10^{-2}$ (red), $10^{-1}$ (blue), and $1\,\text{M}_{\odot}\,\text{yr}^{-1}$ (green). Because this relation is only valid when $t_{\rm acc}\lesssim t_{\rm 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 it can be seen, the temperature at which the photosphere begins to emit H-ionizing radiation, $T_{\rm ion}\simeq 10000\,{\rm 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}\,\text{M}_{\odot}$.

In addition to our constant accretion rate assumption, we have included a time-dependent toy model of the form $\dot{M}_{\star}(t)=1\,\text{M}_{\odot}\,\text{yr}^{-1}e^{-t/t_{\rm ff,0}}$ (purple dash-dotted line), with $t_{\rm ff,0}=10^{5}\,{\rm yr}$ 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\,\text{M}_{\odot}$. We also compare with the Hosokawa et al. (2013) results for the effective temperature $T_{\rm eff}$ at $\dot{M}_{\star}\simeq 1\,\text{M}_{\odot}\,\text{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_{\rm ion}$ during the evolution of the protostar up to masses of $10^{5}\,\text{M}_{\odot}$. Since the photospheric temperature at the characteristic accretion rate $\dot{M}_{\star}\simeq 1\,\text{M}_{\odot}\,\text{yr}^{-1}$ never surpasses $10^{4}$ K, we can safely neglect photo-ionization 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\,\text{M}_{\odot}\,\text{yr}^{-1}$, or for masses $M_{\star}\simeq 10^{5}\,\text{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 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}\,\text{M}_{\odot}$ (see fifth panel of Figure 1).

3 Lessons for sink algorithm

3.1 Accretion radius

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 sub-grid 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

[TABLE]

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:

[TABLE]

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

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_{\rm acc}\gtrsim t_{\rm KH}$ and accretion rates $\dot{M}_{\star}\gtrsim 10^{-2}\,\text{M}_{\odot}\,\text{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

[TABLE]

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:

[TABLE]

The framework here depicted can be used to derive a prescription for the accretion radius of the sink particle, $R_{\rm 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_{\rm acc}\lesssim R_{\rm B}$ if the simulations resolves 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_{\rm th}<n_{\rm B}$. In such 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_{\rm acc}=R_{\rm F}\left(n_{\rm F}/n_{\rm th}\right)^{1/2}$.

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

[TABLE]

For threshold densities of $n_{\rm th}=10^{8}\,\text{cm}^{-3}$ and $n_{\rm th}=10^{10}\,\text{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_{\rm acc}\simeq 3.8\times 10^{4}\,{\rm AU}$ and $R_{\rm acc}\simeq 3.8\times 10^{3}\,{\rm AU}$, respectively. These values are kept until $n_{\rm B}=n_{\rm th}$, which occurs when the mass of the protostar is $\simeq 10^{3}\,\text{M}_{\odot}$ for the former and $\simeq 10^{2}\,\text{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_{\rm th}=10^{8}\,\text{cm}^{-3}$ and $n_{\rm th}=10^{10}\,\text{cm}^{-3}$ are shown in Figure 5 as black dashed line and 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-10^{5}\,\text{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_{\rm 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); Shlosman et al. (2016) assumed $R_{\rm 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_{\rm acc}\sqrt{1.2544c_{\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\text{M}_{\odot}\,\text{yr}^{-1}$, which corresponds to our estimation from Equation (26).

3.2 Cosmological boundary conditions

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 $\Lambda$ cold dark matter ($\Lambda$CDM) cosmology. We then follow the evolution of the halo until the highest density cell reaches the threshold density $n_{\rm th}=10^{8}\,\text{cm}^{-3}$. For that we have used a primordial chemistry network that includes five species (H, $\text{H}_{2}$, ${\rm H}^{-}$, ${\rm H}^{+}$, and ${\rm e}^{-}$) and cooling processes such as ${\rm H}^{-}$ cooling, $\text{H}_{2}$ line cooling, $\text{H}_{2}$ collision-induced emission, Ly-$\alpha$ 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\,{\rm pc}$. This implies that the object reaches spherical symmetry at scales larger than the accretion radius at that point ($R_{\rm acc}\simeq 0.2\,{\rm pc}$, as calculated in Section 3.1), which is plotted in black dashed lines.

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 red solid line, which is calculated as $\dot{M}=-4\pi r^{2}\rho v_{\rm rad}$, with $r$ the distance to the highest density cell, $\rho$ the mass density of hydrogen, and $v_{\rm rad}$ the radial component of the velocity. For reference, we have also included the Shu accretion rate for spherical collapse, $\dot{M}_{\rm Shu}\simeq 0.975c_{\rm s}^{3}/G$ (Shu, 1977), and the Larson-Penston accretion rate for dynamical collapse, $\dot{M}_{\rm LP}\simeq 46.9c_{\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\,\text{M}_{\odot}\,\text{yr}^{-1}$ at $r\simeq 0.5\,{\rm pc}$ and then it decreases to values $\dot{M}\simeq 0.1\,\text{M}_{\odot}\,\text{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}_{\rm Shu}\simeq 0.18\,\text{M}_{\odot}\,\text{yr}^{-1}$ and $\dot{M}_{\rm LP}\simeq 8.5\,\text{M}_{\odot}\,\text{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\,\text{M}_{\odot}\,\text{yr}^{-1}$, which is consistent with the values assumed throughout this study.

3.3 Disk accretion

Throughout the paper, we discuss a sub-grid model within a sink assuming spherical symmetric accretion flows. However, gas material with angular momentum form 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~{}\text{M}_{\odot}\,\text{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 supermassive black holes 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}>10^{2}-10^{3}~{}\text{M}_{\odot}$ with an almost constant effective temperature of $T_{\rm 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~{}\text{M}_{\odot}\,\text{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).

4 Summary and Conclusions

In this paper we have developed a model for the early evolution of supermassive protostars. After the formation of the initial protostar the surrounding gas becomes optically thick to ${\rm H}^{-}$ radiation, at which point we can robustly calculate the properties of the object using the equations of hydrostatic and thermal equilibrium. From that we obtain a characteristic density, radius, and mass of $n_{\rm F}\simeq 4.6\times 10^{16}\,\text{cm}^{-3}$, $R_{\rm F}\simeq 0.33\,{\rm AU}$, and $M_{\rm F}\simeq 0.045\,\text{M}_{\odot}$, respectively, for a temperature $T=3000\,{\rm K}$. An alternative approach to model the same situation is to use one-zone models. For that, we describe in detail the methods introduced in Inayoshi et al. (2014) and provide explicit numerical fits for the ${\rm H}^{-}$ cooling rate and opacity. Combined with the adiabatic heating rate, we can then calculate the critical density at which the gas becomes optically thick, which results in $n_{\rm crit}\simeq 2\times 10^{16}\,{\rm cm}^{-3}$, consistent with the previous estimate. Hence we can robustly characterize the properties of the protostar in the initial optically thick regime.

The early stages of protostellar evolution, where $t_{\rm acc}\lesssim t_{\rm KH}$, are described by the accretion of material onto the central object. For this case, we derive an expression for the protostellar radius as a function of the mass and accretion rate. Using a characteristic value of $\dot{M}_{\star}\simeq 1\,\text{M}_{\odot}\,\text{yr}^{-1}$ for the accretion rate, we find that the protostellar radius grows as $R_{\star}\propto M_{\star}^{1/4}$ during this phase. Once internal sources of radiation start dominating, $t_{\rm KH}\lesssim t_{\rm acc}$ and hence the radius-mass relation changes to $R_{\star}\propto M_{\star}^{1/2}$. For the case of a supermassive protostar, the radius varies from $R_{\star}\simeq 0.65\,{\rm AU}$ at $M_{\star}\simeq 0.1\text{M}_{\odot}$ to $R_{\star}\simeq 250\,{\rm AU}$ at $M_{\star}\simeq 10^{5}\text{M}_{\odot}$. For the surface temperature of the object, we base our analysis on the prescription of Stahler et al. (1986), deducing that it remains well below the ionizing temperature of $T_{\rm ion}\simeq 10^{4}\,{\rm K}$ during most of its evolution. We can thus safely neglect UV ionizing radiation until the late stages of the assembly process.

In numerical simulations, supermassive protostars are commonly represented by sink particles. Our model allows us to derive the properties of such particles and implement a physically-motivated sub-grid model for their evolution in hydrodynamical codes. In particular, we derive an expression for the accretion radius ($R_{\rm acc}$) as a function of the threshold density at which the sink particle is inserted ($n_{\rm th}$), relating it to the physical conditions on the surface of the protostar. For high threshold densities our model proposes a numerical value based on the isothermal profile of the atomic cooling halo, but this value will eventually transition to $R_{\rm acc}=R_{\rm B}$ once the Bondi radius is resolved further in the evolution of the protostar. We can thus verify throughout the simulation that the accretion radius is well adjusted, in the sense that it is larger than the stellar radius at every moment during its evolution. Our new prescription for sink particles implies changes in the early stages of the evolution up to a protostar mass of $\simeq 10^{2-3}\,\text{M}_{\odot}$. After that, it follows the Bondi accretion scenario, and hence it has an accretion rate of $\simeq 1\,\text{M}_{\odot}\,\text{yr}^{-1}$. The final mass of the object corresponds to $\simeq 10^{5-6}\,\text{M}_{\odot}$, similar to previous estimates in the literature. We will track the accretion and KH timescales during the actual simulation to determine when the accretion rate becomes low enough and the star enters the KH phase. At that point, the radiation hydrodynamic effects from an ionizing central source would have to be taken into account.

The ultimate goal of this line of work is to simulate the assembly process of the first supermassive objects in the Universe in an ab-initio fashion. One key question then is: When will this build-up enter a radiation-hydrodynamical phase, where the strong radiative feedback from the growing protostar will eventually turn the object into hyper-luminous beacons from the end of the cosmic dark ages? Those will be probed with next-generation observational facilities, such as the James Webb Space Telescope (JWST), to be launched in 2018. An ideally complementary window into the formation of the first supermassive objects is provided by the gravitational wave signal accompanying the possible merger of binary black holes, which is a prime target for the planned Laser Interferometer Space Antenna (LISA). In light of this suite of next-generation facilities, simulations will have a key role to play in providing physically robust predictions, based on well-motivated sub-grid prescriptions.

Acknowledgements

We thank Kazuyuki Omukai for kindly providing permission to use the numerical models for ${\rm H}^{-}$ cooling described in Section 2.2. KI acknowledges support by the Simons Foundation through the Simons Society of Fellows. VB was supported by NSF grant AST-1413501. We also thank the anonymous referee for the constructive comments that helped to improve our paper.

Following our discussion in Section 2.2, we can derive numerical fits to the ${\rm H}^{-}$ cooling rate based on the description of Equations (12) to (16) as introduced by Inayoshi et al. (2014). The terms in the right-hand side of Equation (13) can be approximated as $\Lambda^{(h)}_{\rm thin}=k^{(h)}n_{\rm HI}n_{\rm e}$ and $\Lambda^{(l)}_{\rm thin}=k^{(l)}n_{\rm HI}n_{\rm e}$, respectively. The cooling rate coefficients are given by

[TABLE]

where $T_{3}=T/10^{3}\,{\rm K}.$

Additionally, we estimate the Planck mean opacity for the ${\rm H}^{-}$ free-free emission in both the low ($\kappa^{(l)}_{\rm ff,P}$) and high ($\kappa^{(h)}_{\rm ff,P}$) frequency regime, and for the ${\rm H}^{-}$ bound-free emission for the high ($\kappa^{(h)}_{\rm bf,P}$) frequency regime as

[TABLE]

Another opacity source is the Rosseland mean opacity for ${\rm H}^{-}$ free-free emission in the low frequency regime, which can be modeled as

[TABLE]

Finally, we also estimate the opacity due to Rayleigh scattering for the high frequency regime as

[TABLE]

The total opacities for both regimes can then be written as $\kappa^{(h)}_{\rm R}=\kappa^{(h)}_{\rm Ray}$, $\kappa^{(l)}_{\rm R}=\kappa^{(l)}_{\rm ff,R}$, $\kappa^{(h)}_{\rm P}=\kappa^{(h)}_{\rm ff,P}+\kappa^{(h)}_{\rm bf,P}$, and $\kappa^{(l)}_{\rm P}=\kappa^{(l)}_{\rm ff,P}$. With these approximations and Equation (16) we can rewrite the total ${\rm H}^{-}$ cooling rate as a function of the optically thin cooling rates and the opacities:

[TABLE]

Here the characteristic length $\ell$ has been set to the Jeans length $\ell=\lambda_{\rm J}\simeq c_{\rm s}t_{\rm ff}$, with $c_{\rm s}=\sqrt{\gamma k_{\rm B}T/\mu m_{\rm H}}$ being the sound speed.