© ESO 2019

1. Introduction

The formation of the first supermassive black holes (SMBHs) is one of the main problems in our understanding of the formation of cosmic structures (e.g. Rees 1978, 1984; Volonteri & Begelman 2010; Woods et al. 2019). The most extreme SMBHs have M_•  ∼  10⁹ M_⊙ at redshift z  ∼  7 (Mortlock et al. 2011; Wu et al. 2015; Bañados et al. 2018; Wang et al. 2018). Assuming the seed of these objects formed at z  ∼  30, for example through direct collapse, this value implies the accumulation of mass at average rates of 1–10 M_⊙ yr⁻¹ in compact objects for a billion years. Such rates are found in atomically cooled primordial haloes, where H₂ molecules have been destroyed by an external UV field (e.g. Haiman et al. 1997; Omukai 2001; Bromm & Loeb 2003; Dijkstra et al. 2008; Latif et al. 2013; Regan et al. 2017).

More recently, an alternative route has been found in the merger of massive galaxies (Mayer et al. 2010, 2015; Mayer & Bonoli 2019). This scenario has two advantages compared to atomically cooled haloes: first, star formation does not need to be suppressed and second, much larger inflows are obtained, as high as 10⁴ − 10⁵ M_⊙ yr⁻¹ at sub-parsec scales. The gas accumulates on a rotationally supported disc that reaches ∼10⁹ M_⊙; what happens next however remains unknown. The disc might either collapse directly to a black hole (dark collapse, Mayer & Bonoli 2019), or first form a supermassive star (SMS), which then collapses to form a SMBH once accretion stops or general relativistic (GR) instability sets in (e.g. Woods et al. 2019). Interestingly, the most distant observed quasar seems to be hosted by a galaxy merger (Bañados et al. 2019) and has an inferred mass similar to that of the massive disc obtained in the simulations.

The properties of SMSs accreting at the rates of atomically cooled haloes (0.1–10 M_⊙ yr⁻¹) have been studied both analytically and numerically in the last decade (Begelman 2010; Hosokawa et al. 2012, 2013; Schleicher et al. 2013; Sakurai et al. 2015; Umeda et al. 2016; Woods et al. 2017; Haemmerlé et al. 2018a,b; Haemmerlé & Meynet 2019). They are found to evolve as red supergiant protostars along the Hayashi limit (Hosokawa et al. 2012, 2013), before collapsing due to the GR instability at several 10⁵ M_⊙ (Umeda et al. 2016; Woods et al. 2017; Haemmerlé et al. 2018a). The case of galaxy mergers provides conditions for which rates higher than those of atomically cooled haloes are possible. Here we address the question of whether or not SMSs can form by accretion at rates ≳100 M_⊙ yr⁻¹. We investigate the conditions in which hydrostatic equilibrium can be sustained for such rates by computing the corresponding hydrostatic structures numerically and estimating a posteriori the first-order hydrodynamical corrections. We consider the stability of monolithic SMSs in the context of other formation scenarios in a companion paper (Woods et al., in prep.). For a discussion on the bottlenecks expected in SMS formation, see (Woods et al. 2019).

The paper is organised as follows. In Sect. 2 we proceed to a preliminary analysis based on global timescales in order to illustrate the main ideas. In Sect. 3, we give the numerical method followed to build the stellar models. The stellar structures are described in Sect. 4 and their implications are discussed in Sect. 5. We conclude in Sect. 6.

2. Timescales

In the absence of nuclear reactions, the evolution of an accreting star is governed by the competition between two processes: thermal relaxation and accretion. The first process, appearing mathematically in the energy equation, relies on the variation of entropy ṡ. The second one, appearing in the boundary conditions, relies on the variation of mass Ṁ. The efficiency of these two processes can be approximately estimated from their corresponding global timescales, the Kelvin–Helmholtz (KH) time and the accretion time:

$$\begin{array}{cc}& {\tau }_{\mathrm{KH}}=\frac{G{M}^2}{\mathit{RL}},\hfill \end{array}$$(1)

$$\begin{array}{cc}& {\tau }_{\mathrm{accr}}=\frac{M}{\dot{M}},\hfill \end{array}$$(2)

where M is the stellar mass, R and L are the photospheric radius and luminosity, respectively, and G is the gravitational constant. When τ_accr ≫ τ_KH, the evolution is essentially governed by thermal contraction and the star is expected to be thermally relaxed after a KH time, as nuclear reactions are ignited. When τ_accr ≪ τ_KH, accretion dominates the evolution, leading to the swelling of the envelope, because the entropy cannot be radiated away efficiently enough to restore thermal equilibrium against the perturbative effect of accretion (Hosokawa & Omukai 2009).

In principle, if the accretion time is short enough compared to the KH time (i.e. if accretion is fast enough), thermal processes become negligible (in terms of structure) and the evolution is adiabatic: it is not a loss of entropy that causes the star to contract, but rather the increasing pressure caused by the accreted mass. This implies that above a given threshold in the accretion rate, the stellar structure for a given mass no longer depends on the accretion rate. Indeed, if the evolution is governed by accretion only, the time only appears in the problem through the change in mass. In the range Ṁ ≤ 10 M_⊙ yr⁻¹ already covered in the literature, thermal processes are still at play and this pure accretion regime is not reached.

The existence of a pure accretion regime depends on the ability of the star to maintain hydrostatic equilibrium at the required rates. Mechanical equilibrium in stellar interiors is restored by pressure- and gravity-waves in a sound-crossing time:

$$\begin{array}{c}\hfill {\tau }_{\mathrm{SC}}=\frac{R}{v_{\mathrm{sound}}}=\frac{R}{\sqrt{{\mathrm{\Gamma }}_1\frac{P}{\rho }}},\end{array}$$(3)

where v_sound is the sound-speed, Γ₁ the first adiabatic exponent, P the thermal pressure, and ρ the density of mass. If τ_accr ≲ τ_SC, the mass increases by a factor of two or more during the time needed to restore hydrostatic equilibrium. Thus a hydrostatic core cannot incorporate mass at such rates.

Figure 1 shows the timescales of Eqs. (1)–(3). The KH time is obtained using the fits by Hosokawa et al. (2012) for R and L:

$$\begin{array}{c}\hfill R=260R_{\odot }\times {\left(\frac{M}{M_{\odot }}\right)}^{1/2}L=38000L_{\odot }\times \frac{M}{M_{\odot }}·\end{array}$$(4)

Fig. 1. Global timescales of Eqs. (1)–(3), computed with Eq. (4). The sound-crossing time is shown for ideal gas and radiation pressure. The accretion time is shown for the indicated rates.

For the sound-crossing time, we estimate the sound speed for cases of ideal gas and radiation pressure. We assume constant mean molecular weight (μ = 0.6) and temperature (10⁵ K) relevant for the envelope of accreting SMSs. For the density, we consider only the outer 10% of the stellar mass, which covers most of the radius. For M >  10⁴ M_⊙, the sound-crossing time is longer than a year due to the large radii and low temperatures. We note the effect of radiation pressure, which dominates in this mass range, in shortening the sound-crossing time by an order of magnitude. A comparison with the KH time indicates that a regime of pure accretion might exist in the high mass range for rates around 10³ − 10⁴ M_⊙ yr⁻¹.

The use of global timescales in stellar evolution is not always conclusive however. For instance, Fig. 1 suggests that stars of 1000 M_⊙ accreting at 1 M_⊙ yr⁻¹ are thermally relaxed, which is refuted by all numerical models (e.g. Woods et al. 2019 for a comparison). For the same reason, a precise estimate of the maximal accretion rate that allows for hydrostatic equilibrium requires the consideration of numerical models.

3. Numerical method

The numerical models are computed with the same method as Haemmerlé et al. (2018a), that is, with GENEC, a one-dimensional hydrostatic stellar evolution code that solves the equations of stellar structure numerically. Detailed descriptions of the code are available in the literature (e.g. Eggenberger et al. 2008). In the present case, the central ingredient is accretion (Haemmerlé et al. 2016), included through cold accretion; that is the entropy of the accreted gas matches that of the stellar surface (Palla & Stahler 1992). For rates of atomically cooled haloes, the effect of hot accretion (i.e. advection of entropy; Hosokawa & Omukai 2009) can be neglected since the intrinsic luminosity, which is nearly Eddington, dominates the energy budget at the stellar surface (Hosokawa et al. 2013). For larger rates, this assumption might not be justified, which is discussed in Sect. 5.2. Rotation is not included, since accreting SMSs are found to be slow rotators due to the ΩΓ-limit, which implies negligible impact of rotation on their structure (Haemmerlé et al. 2018b).

For zero metallicity, models with Ṁ ≤ 10 M_⊙ yr⁻¹ are already described in Haemmerlé et al. (2018a). Models of rapidly accreting SMSs at solar metallicity have not yet been published. For the present study, we compute four new models. We first extend the set of Pop III models by considering the rates 100 M_⊙ yr⁻¹ and 1000 M_⊙ yr⁻¹. In order to understand the impact of metals in various regimes, we compute two models at solar metallicity with rates 1 and 1000 M_⊙ yr⁻¹. We highlight the fact that we do not include any mass-loss, even at solar metallicity, due to the short evolutionary timescales. We use initial models of 5 and 10 M_⊙ for numerical stability, with fully convective structures similar to those used in Haemmerlé et al. (2018a). The models are run until the limit of numerical stability, which is suspected to be related to pulsation instability.

4. Results: hydrostatic structures

The internal structures and surface properties of the new models are qualitatively similar to those of the models already published. They evolve along the Hayashi limit and start to burn H in their core at the latest when they reach 10⁵ M_⊙. Hydrogen burning triggers convection in the core, but most of the mass remains radiative. Hydrostatic objects forming at rates ≲1000 M_⊙ yr⁻¹ go through energetically significant nuclear reactions, and are thus stars, sensu stricto.

The main properties of the models are summarised in Tables 1–3. The highest mass we reach is 800 000 M_⊙ for the model of 1000 M_⊙ yr⁻¹ at solar metallicity. The corresponding lifetime is as short as 800 yr. The mass at which H starts to burn increases with the accretion rate until this latter reaches at least 100–1000 M_⊙ yr⁻¹. However, the mass fraction of the convective core at a given mass depends weakly on the rate above 100 M_⊙ yr⁻¹. For Z = 0, this fraction remains always below 1%; at solar metallicity, it grows to at most 8% at the end of therun.

The typical structure of a star accreting at 1000 M_⊙ yr⁻¹ is illustrated in Fig. 2, with the internal profiles of luminosity and entropy of the model at solar metallicity, taken near the endpoint of the run, at a mass of 8 × 10⁵ M_⊙. The convective core contains 8% of the total stellar mass. Most of the mass is in the radiative zone, where the luminosity profile matches the local Eddington luminosity, L_Edd(r) = 4πcGM_r/κ (c is the speed of light and κ the opacity given by the model, which is dominated by electron-scattering), and the entropy profile follows a power law $s\propto M_{\mathrm{r}}^{1/2}$. This power law corresponds to the hylotropic structures of Begelman (2010). The entropy profiles at successive stages match each other in the radiative region, indicating negligible entropy losses, that is, a pure accretion regime in most of the star.

Fig. 2. Internal profiles of the model at Ṁ = 1000 M⊙ yr−1, Z = Z⊙ and M = 8 × 105 M⊙. Upper panel: luminosity profile (solid blue) and Eddington luminosity (dotted white). Lower panel: entropy profile (solid blue) and hylotropic fit $s\propto M_{\mathrm{r}}^{1/2}$ (dotted white).

Of course, this is not the case in the convective core, where the entropy released by H-burning is so huge that the local luminosity exceeds the local Eddington luminosity by three orders of magnitude, and even the surface Eddington limit by one order of magnitude (Fig. 2, upper panel). We note that the radiative luminosity is never super-Eddington in the core, meaning that density or temperature inversions are not required. Thermal processes are so inefficient that the released nuclear entropy cannot be transported by radiation and remains captured in the convective layers. As a consequence, the luminosity drops by several orders of magnitude at the boundary of the core, and converges rapidly to the Eddington profile. It follows that the evolution of the envelope is decoupled from that of the core. In particular, the energy radiated at the surface arises from the thermal contraction of the envelope, not from central H-burning. This is in spite of the fact that thermal processes have negligible impact on the structure. Indeed, the dominating process (adiabatic accretion-driven contraction) does not release any entropy, and thus does not contribute to the luminosity.

The impact of metals on stellar structures is twofold: they increase the opacity and allow for the CNO cycle to operate without previous 3α reactions. In the present case, the effect of high opacity is to lock the star on the Hayashi limit for longer, avoiding the drift towards the blue, as obtained for ≳10⁹ L_⊙ in previous studies at zero metallicity (Hosokawa et al. 2013; Haemmerlé et al. 2018a). This is due to the large number of free electrons provided by the metals that are available for H⁻ formation, whose opacity law keeps the photospheric temperature constant. This fact implies a weaker UV feedback in the supermassive range compared to zero metallicity. In solar metallicity models, the possibility for the CNO cycle to start before 3α reactions implies that H-burning is triggered at lower central temperatures, that is, at lower masses (Table 2). The main consequence of this is an increase in the mass fraction of the convective core for a given total mass (Table 3) since entropy starts to be released earlier. Another consequence is a lower central temperature and density for a given mass, since the pre-H-burning contraction is stopped earlier.

The final collapse of SMSs relies on the GR instability, which is a radial pulsation instability and cannot be captured by hydrostatic models. Polytropic models provide a simple criterion (Chandrasekhar 1964) however, which is not conclusive in the present case. This criterion might apply only in the core, which corresponds to a polytrope. For rates of atomically cooled haloes, it fails to predict correct final masses, but still allows the trend of the dependence on the accretion rate to be captured (Woods et al. 2017, 2019; Haemmerlé et al. 2018a). The final masses given by the polytropic criterion in the present case remain in the range of those of previous studies at lower rates, which confirms that above 0.1 M_⊙ yr⁻¹ the final mass depends weakly on the rate. For Z = 0, the instability arises at ∼5 × 10⁵ M_⊙ in the models at 100–1000 M_⊙ yr⁻¹. The models at solar metallicity have larger cores but also lower densities (i.e. lower GR correction), which can delay the instability by a factor of approximately two in mass and age. Overall, if we consider the polytropic criterion in the core only, the most massive stable model we obtain is for 1000 M_⊙ yr⁻¹ and Z = Z_⊙, and has a mass of 6.8 × 10⁵ M_⊙.

5. Discussion

5.1. Hydrostatic equilibrium

As discussed in Sect. 2, hydrostatic equilibrium requires that evolution proceed slowly enough for the pressure- and gravity waves to cross the star before it evolves significantly. Since for the rates considered here the accretion time becomes similar to the sound-crossing time, the consistency of the models must be questioned.

Hydrodynamical effects become important when the acceleration term ṙ in the momentum equation,

$$\begin{array}{c}\hfill \ddot{r}=-\frac{1}{\rho }\mathrm{\nabla }P-g,\end{array}$$(5)

becomes comparable with the other terms, the gravitational acceleration g >  0 and the acceleration due to the pressure gradient ∇P. In order to estimate the departures from hydrostatic equilibrium, we derive hydrodynamical structures from the hydrostatic structures described in Sect. 4, by assuming each layer M_r has, at each age t, the properties (P,  ρ,  g) of the hydrostatic structure at a previous age t′, with a time-delay t − t′ corresponding to the time it takes for a sound wave to cross the layers above M_r. Modifications in the boundary conditions (i.e. accretion) during this delay cannot impact the layers below M_r because this would require supersonic transport of momentum. Since the evolution is driven by accretion, a layer has no reason to evolve except under the effect of accretion. The time-delay is computed by integration of the inverse of the sound speed (denominator of Eq. (3)) over a path comoving with the contracting gas. Assuming the modified profiles of P, ρ, and g mimic the hydrodynamical structures at first order, we can derive the acceleration term with Eq. (5). We note that ∇P is re-derived from the new profiles: adjacent layers are delayed to different hydrostatic models, which induces a correction to ∇P at the origin of the non-zero acceleration term in Eq. (5). The condition for hydrostatic equilibrium is then obtained by comparing this acceleration term with the gravitational term. We assume hydrodynamical effects become significant when $|\ddot{r}|\gtrsim 0.1\mathit{g}$. This method does not allow the actual hydrodynamical evolution to be captured however, but the underlying idea is the following: if, even when accounting for the time it takes for the various layers to come under the effects of accretion, the corrected structures do not induce significant acceleration, and the evolution remains hydrostatic.

The larger the stellar mass, the weaker the impact of accretion of a given amount of mass on the structure. Therefore, for a constant accretion rate the evolution slows down as the mass grows, which is reflected by the increase of the accretion time with mass (Eq. (2), Fig. 1). As a consequence, the star reaches stability more easily as its mass grows, and the condition for hydrostatic equilibrium corresponds to a minimum mass for any given rate. Moreover, due to their low temperatures, large radii, and rapid contraction, the outer layers are the least stable, meaning that hydrostatic equilibrium is broken at higher masses in the envelope than in the core. Thus for each accretion rate, we derive the minimum mass at which the whole star becomes hydrostatic as well as that at which the core only can stabilise.

Figure 3 summarises the limits to hydrostatic equilibrium for the case at zero metallicity. The models at solar metallicity do not exhibit significant differences. The limit to equilibrium arising from accretion is shown by the two black curves, between which only a fraction of the star can stabilise. Departures from equilibrium remain negligible for rates ≤1 M_⊙ yr⁻¹, and the model at 10 M_⊙ yr⁻¹ becomes stable as soon as it exceeds 100 M_⊙. Hydrodynamical effects become significant for rates ≥100 M_⊙ yr⁻¹. The core can only stabilise above a stellar mass of 100–1000 yr⁻¹ M_⊙, and the star as a whole above ≳10⁴ M_⊙. Therefore, accretion at rates ≳100 M_⊙ yr⁻¹ is only compatible with mechanical equilibrium for stars more massive than 10⁴ M_⊙. The red line on Fig. 3 shows the limits to equilibrium resulting from GR instability obtained by Woods et al. (2017). An extrapolation of the curve (see also Sect. 4) suggests that a hylotropic, pure-accretion regime can be stable over one or two orders of magnitude in mass. Additional effects, such as rotation, could extend the mass range further by stabilising the star against the GR instability (Fowler 1966; Bisnovatyi-Kogan et al. 1967).

Fig. 3. Limits to hydrostatic equilibrium for accretion above 10−3 M⊙ yr−1 (without rotation). Sound waves can maintain full hydrostatic equilibrium only below the solid black line. Above the dashed black line, the whole mass evolves hydrodynamically. Stars become GR unstable at the right-hand side of the red dashed line. The maximum rate of atomically cooled haloes is shown by a grey dashed line. The location of the hylotropic, pure-accretion regime is indicated.

5.2. Hylotropic structures

Hylotropic models have been developed by Begelman (2010) to analytically study the structures of accreting SMSs. These were applied to the quasi-star model (Begelman et al. 2008), where a hydrostatic envelope is supported by the energy of accretion onto a central BH. The stability of such systems requires the BH to be ≲1% of the envelope mass. The question remains open as to whether or not accreting SMSs collapsing through the GR instability can form such objects (Hosokawa et al. 2013; Woods et al. 2019). Since the GR instability is thought to appear in the core first, a core of ≲1% of the total mass is probably required. The models of SMSs accreting at the rates of atomically cooled haloes have larger cores (≳10% M), but interestingly in the models of the present study at 100–1000 M_⊙ yr⁻¹ and Z = 0 the core mass has the required fraction, while the rest of the star corresponds indeed to a hylotrope.

We note however that for the rates required for hylotropic structures, the gravitational energy liberated by accretion becomes similar or exceeds the energy radiated at the photosphere. Indeed, with the fits of Eq. (4), the ratio of the accretion luminosity L_accr = GMṀ/R to the intrinsic luminosity is

$$\begin{array}{c}\hfill \frac{L_{\mathrm{accr}}}{L}\sim {\left(\frac{M}{{10}^6{M}_{\odot }}\right)}^{-1/2}\frac{\dot{M}}{100M_{\odot }{\mathrm{yr}}^{-1}},\end{array}$$(6)

which shows that hot accretion (Sect. 3) becomes dominant if the rate exceeds 100 M_⊙ yr⁻¹ before the mass reaches 10⁶ M_⊙. Nevertheless, assuming that the heat advected per mass unit, TΔs_hot, scales with the specific energy of free-fall, GM/R, and that the temperature at accretion is constant (Hayashi limit), the fit of Eq. (4) for R implies that $\mathrm{\Delta }{s}_{\mathrm{hot}}\propto M_{\mathrm{r}}^{1/2}$, that is, the structure remains hylotropic.

6. Summary and conclusions

The models of the present study extend the previous works on rapidly accreting SMSs by considering larger rates (100–1000 M_⊙ yr⁻¹) and solar chemical composition. The new stellar structures display most of the main features of the previous models: the star evolves as a red supergiant protostar with most of the mass being radiative, and burns hydrogen in a convective core before GR instability leads to collapse, forming a SMBH. The polytropic criterion indicates an onset of GR instability at masses of 10⁵ − 10⁶ M_⊙.

The main quantitative differences are in the size of the convective core, which represents only ≲1% of the total stellar mass for the models at 100–1000 M_⊙ yr⁻¹ and zero metallicity. Adding metals at solar abundances changes the core mass fraction by a factor of a few. For rates above the atomic cooling limit (≳10 − 100 M_⊙ yr⁻¹), thermal processes become inefficient and the evolution is governed by accretion only, leading to hylotropic structures.

However, the evolution at ≳100 M_⊙ yr⁻¹ is so fast that hydrostatic equilibrium can only be sustained if the star is already supermassive. At lower masses, pressure- and gravity-waves cannot restore hydrostatic equilibrium fast enough. This implies that stars can become supermassive by accretion only at the rates of atomically cooled haloes. Once they are supermassive, larger rates are possible. This result implies that forming hylotropic stars by accretion requires rates that increase with time. The evolution in the dynamical accretion regime, in particular the possibility for dark collapse, shall be addressed in a future study.

Acknowledgments

This work was sponsored by the Swiss National Science Foundation (project number 200020-172505). RSK acknowledges financial support from the German Science Foundation (DFG) via the collaborative research centre (SFB 881) “The Milky Way System” (subprojects B1, B2, and B8) and from the Heidelberg cluster of excellence EXC 2181 “STRUCTURES: A unifying approach to emergent phenomena in the physical world, mathematics, and complex data” funded by the German Excellence Strategy. TEW acknowledges support from the NRC-Canada Plaskett Fellowship. We thank Kazuyuki Omukai and Takashi Hosokawa for fruitful discussions.

References

All Tables

All Figures

	Fig. 1. Global timescales of Eqs. (1)–(3), computed with Eq. (4). The sound-crossing time is shown for ideal gas and radiation pressure. The accretion time is shown for the indicated rates.
In the text

	Fig. 2. Internal profiles of the model at Ṁ = 1000 M⊙ yr−1, Z = Z⊙ and M = 8 × 105 M⊙. Upper panel: luminosity profile (solid blue) and Eddington luminosity (dotted white). Lower panel: entropy profile (solid blue) and hylotropic fit $s\propto M_{\mathrm{r}}^{1/2}$ (dotted white).
In the text

Fig. 3. Limits to hydrostatic equilibrium for accretion above 10−3 M⊙ yr−1 (without rotation). Sound waves can maintain full hydrostatic equilibrium only below the solid black line. Above the dashed black line, the whole mass evolves hydrodynamically. Stars become GR unstable at the right-hand side of the red dashed line. The maximum rate of atomically cooled haloes is shown by a grey dashed line. The location of the hylotropic, pure-accretion regime is indicated.

In the text