© ESO 2021

1. Introduction

The possibility for stars with masses ≳10⁵ M_⊙ to exist in a stable state of equilibrium has been considered for half a century, in connection with the physics of quasars (e.g., Hoyle & Fowler 1963b,a; Fowler 1966; Bisnovatyi-Kogan et al. 1967; Appenzeller & Kippenhahn 1971; Appenzeller & Fricke 1971, 1972a,b; Shapiro & Teukolsky 1979; Fuller et al. 1986; Baumgarte & Shapiro 1999b,a; Begelman 2010; Hosokawa et al. 2013; Umeda et al. 2016; Woods et al. 2017; Haemmerlé et al. 2018b,a, 2019; see Woods et al. 2019 and Haemmerlé et al. 2020 for recent reviews). One of the particularities of these supermassive stars (SMSs) is to end their life through the general-relativistic (GR) instability (Chandrasekhar 1964). The pressure support of SMSs is dominated by radiation (∼99%), so that these stars are always close to the Eddington limit. Since the Eddington limit corresponds to a state of marginal stability in Newtonian gravity, small GR corrections, of the order of a percent, are sufficient to make the star unstable with respect to radial pulsations. For Population III (Pop III) SMSs, the subsequent collapse is expected to lead to the direct formation of a supermassive black hole (Fricke 1973; Shapiro & Teukolsky 1979; Shibata & Shapiro 2002; Shapiro & Shibata 2002; Liu et al. 2007; Uchida et al. 2017; Sun et al. 2017, 2018).

Such direct collapse black holes are invoked to explain the existence of the most massive quasars observed around redshift 7 (e.g., Rees 1978, 1984; Volonteri & Begelman 2010; Volonteri 2010; Valiante et al. 2017; Zhu et al. 2020). The inferred black hole masses in these objects are in excess of 10⁹ M_⊙, so that accretion must have proceeded at average rates ≳1−10 M_⊙ yr⁻¹ because of the finite age of the Universe. The most studied channel of direct collapse is the case of primordial, atomically cooled haloes, in which mass inflows of 0.1–10 M_⊙ yr⁻¹ are expected below a parsec (e.g., Bromm & Loeb 2003; Dijkstra et al. 2008; Latif et al. 2013; Regan et al. 2016, 2017; Chon et al. 2018; Patrick et al. 2020). Larger inflows, up to ∼10⁵ M_⊙ yr⁻¹, have been found in the simulations of the mergers of massive galaxies at redshifts 8–10 (Mayer et al. 2010, 2015; Mayer & Bonoli 2019). Such galaxies are thought to have experienced significant star formation, so that the progenitors of direct collapse black holes in this channel might be Population I (Pop I) SMSs. The hydrodynamical simulations of unstable SMSs show that for metallicities as high as solar, the collapse generally triggers a thermonuclear explosion that prevents black hole formation (Appenzeller & Fricke 1972a,b; Fuller et al. 1986; Montero et al. 2012). Only for masses ≳10⁶ M_⊙ does the deep potential well allow for gravity to bind the collapsing star. In both cases, direct collapse or supernova explosion, the high energies involved in the death of a SMS might allow for observational signatures of the existence of these objects. The simulations of collapsing Pop III SMSs suggest the possibility for gravitational wave emission and ultra-long gamma-ray bursts (Liu et al. 2007; Shibata et al. 2016a; Uchida et al. 2017; Sun et al. 2017, 2018; Li et al. 2018). But the detectability of such signatures sensitively depends on the mass and rotational properties of the SMS at collapse. In particular, both gravitational waves and gamma-ray bursts require spherical symmetry to be broken, so that rotation is crucial in that respect.

The properties of SMSs at the onset of GR instability have been studied in the case of polytropic structures (e.g., Fowler 1966; Bisnovatyi-Kogan et al. 1967; Baumgarte & Shapiro 1999a; Shibata et al. 2016b; Butler et al. 2018). Stars with masses ∼10⁸ − 10⁹ M_⊙ are found to be stable, provided rapid enough rotation. Such structures correspond to fully convective, thermally relaxed stars, which are assumed to have formed at once (‘monolithic’ SMSs, Woods et al. 2020). They are found to feature universal rotational properties, in particular a universal value for the spin parameter, which is key regarding black hole formation and gravitational wave emission. But polytropic models become irrelevant when we account for the formation process of SMSs. The minimum accretion rates ∼1 M_⊙ yr⁻¹ that are required to explain the existence of the most massive and distant quasars are also necessary for the formation of a SMS, because the H-burning timescale near the Eddington limit is expected to be of the order of a million years, so that gathering 10⁵ M_⊙ before fuel exhaustion requires rates ≳0.1 M_⊙ yr⁻¹. Accretion at such rates critically impacts the structure of SMSs by preventing the thermal relaxation of their envelope (Begelman 2010; Hosokawa et al. 2012, 2013; Schleicher et al. 2013; Sakurai et al. 2015; Haemmerlé et al. 2018b, 2019). The high entropy in the outer regions stabilises most of the star against convection, which develops only in the inner ∼10% of the mass, triggered by H-burning. The stability of such structures against GR has been addressed in several works (Umeda et al. 2016; Woods et al. 2017; Haemmerlé 2020, 2021). SMSs accreting at rates from 0.1 to 10 M_⊙ yr⁻¹ are found to collapse at masses 1−4 × 10⁵ M_⊙, while masses in excess of 10⁶ M_⊙ would require rates ≳1000 M_⊙ yr⁻¹. However, these results do not account for the effect of rotation.

Accretion onto a hydrostatic structure requires low enough angular momentum for the centrifugal force to be cancelled by an excess of gravity over the pressure forces. The rotation velocity of SMSs is constrained by the ΩΓ-limit (Maeder & Meynet 2000), which imposes velocities of less than ∼10% of the Keplerian values for masses ≳10⁵ M_⊙, due to the prominent role of radiation pressure in this mass range (Haemmerlé et al. 2018a; Haemmerlé & Meynet 2019). For such slow velocities, rotation hardly impacts the stellar structure and, in particular, rotational flattening is negligible. As we see later, the dynamical effects of rotation are of the same order as those of gas pressure and GR corrections, so that, similar to these effects, rotation plays a key role in the stability of the star, but can be neglected regarding its equilibrium properties.

In the present work, we use this advantage to define rotation profiles on the non-rotating structures of rapidly accreting SMSs on the basis of the hylotropic models of Begelman (2010). In the absence of convection, large differential rotation is expected to develop in the interior of SMSs (Haemmerlé et al. 2018a; Haemmerlé & Meynet 2019). Thus, while solid rotation must be imposed in the convective core, we assume that each radiative layer contracts with local angular momentum conservation. In Haemmerlé (2020), we applied the relativistic pulsation equation of Chandrasekhar (1964) to non-rotating hylotropes. Here, we add a term accounting for rotation to this equation, following the method of Fowler (1966), and apply it to the rotating hylotropes. The hydrostatic structures and rotation profiles are defined in Sects. 2.1 and 2.2, respectively. The rotation term for the GR instability is introduced in Sect. 2.3 (see also Appendix A). The numerical results are presented in Sect. 3 and their implications are discussed in Sect. 4. We summarise the main conclusions in Sect. 5.

2. Method

2.1. Hylotropic structures

Classical spherical stellar structures are determined by the equations of hydrostatic equilibrium and continuity of mass in the following form:

$$\begin{array}{cc}& \frac{\mathrm{d}P}{\mathrm{d}r}=-\rho \frac{G{M}_r}{r^2}\hfill \end{array}$$(1)

$$\begin{array}{cc}& \frac{\mathrm{d}{M}_r}{\mathrm{d}r}=4\pi r^2\rho \hfill \end{array}$$(2)

where r is the radial distance, M_r is the mass enclosed inside r, P is the pressure, ρ is the density of mass, and G is the gravitational constant. The structure is fully determined, provided an additional constraint of the form f(r, M_r, P, ρ) = 0, for example a constraint on the entropy profile: s(P, ρ) = s(M_r). In convective regions, where entropy is uniformly distributed, a constraint is given by s(M_r) = cst. But the entropy gradients that characterise radiative regions generally prevent from closing the equation system without the use of additional differential equations, energy transfer, and conservation.

For rapid enough accretion, the entropy profile of SMSs is well described by a hylotropic law (Begelman 2010; Haemmerlé et al. 2019; Haemmerlé 2020):

$$\begin{array}{c}\hfill s=\{\begin{array}{cc}{s}_c& \mathrm{if}{M}_r<M_{\mathrm{core}}\\ s_c{\left(\frac{M_r}{M_{\mathrm{core}}}\right)}^{1/2}& \mathrm{if}{M}_r>M_{\mathrm{core}}\end{array}\end{array}$$(3)

where M_core is the mass of the isentropic, convective core and s_c is the specific entropy in the core. Equation (3) provides the required constraint s(P, ρ) = s(M_r) that closes (1) and (2), with a function s(P, ρ) given by the equation of state.

The entropy of SMSs is dominated by radiation,

$$\begin{array}{c}\hfill s\simeq s{}_{\mathrm{rad}}=\frac{4}{3}\frac{a_{\mathrm{SB}}{T}^3}{\rho },\end{array}$$(4)

where T is the temperature and a_SB is the Stefan-Boltzmann constant. The contribution from gas reduces to its contribution to pressure:

$$\begin{array}{c}\hfill P=P{}_{\mathrm{gas}}+P{}_{\mathrm{rad}}=\frac{k_{\mathrm{B}}}{\mathrm{\mu }{m}_{\mathrm{H}}}\rho T+\frac{1}{3}{a}_{\mathrm{SB}}{T}^4,\end{array}$$(5)

where μ is the mean molecular weight, m_H is the mass of a proton, and k_B is the Boltzmann constant. This contribution is expressed by the following ratio (∼1%):

$$\begin{array}{c}\hfill \beta :=\frac{P{}_{\mathrm{gas}}}{P}\simeq \frac{P{}_{\mathrm{gas}}}{P{}_{\mathrm{rad}}}=\frac{3}{a_{\mathrm{SB}}}\frac{k_{\mathrm{B}}}{\mathrm{\mu }{m}_{\mathrm{H}}}\frac{\rho }{T^3}\simeq \frac{k_{\mathrm{B}}}{\mathrm{\mu }{m}_{\mathrm{H}}}\frac{4}{s}.\end{array}$$(6)

Inserting Eq. (3) into Eq. (6), we obtain

$$\begin{array}{c}\hfill \beta =\{\begin{array}{cc}{\beta }_c& \mathrm{if}{M}_r<M_{\mathrm{core}}\\ {\beta }_c{\left(\frac{M_r}{M_{\mathrm{core}}}\right)}^{-1/2}& \mathrm{if}{M}_r>M_{\mathrm{core}},\end{array}\end{array}$$(7)

where ${\beta }_c=\frac{k_{\text{B}}}{\mathrm{\mu }{m}_{\text{H}}}\frac{4}{s_c}$. Pressure can be expressed as a function of ρ and β instead of ρ and T, using the definition (6) of β to eliminate T in Eq. (5):

$$\begin{array}{c}\hfill P\simeq {\left(\frac{3}{a_{\mathrm{SB}}}\right)}^{1/3}{\left(\frac{k_{\mathrm{B}}}{\mathrm{\mu }{m}_{\mathrm{H}}}\right)}^{4/3}\frac{{\rho }^{4/3}}{{\beta }^{4/3}}.\end{array}$$(8)

Inserting Eq. (7) into Eq. (8), we obtain

$$\begin{array}{c}\hfill P=\{\begin{array}{cc}K{\rho }^{4/3}& \mathrm{if}{M}_r<M_{\mathrm{core}}\\ K{\left(\frac{M_r}{M_{\mathrm{core}}}\right)}^{2/3}{\rho }^{4/3}& \mathrm{if}{M}_r>M_{\mathrm{core}}\end{array},\end{array}$$(9)

where

$$\begin{array}{c}\hfill K={\left(\frac{3}{a_{\mathrm{SB}}}\right)}^{1/3}{\left(\frac{k_{\mathrm{B}}}{\mathrm{\mu }{m}_{\mathrm{H}}}\right)}^{4/3}{\beta }_c^{-4/3}.\end{array}$$(10)

The power-law (9) allows one to close the system (1) and (2), and thus to solve the stellar structure numerically for any length or mass scales, in a way similar to polytropes (Begelman 2010). To that aim, we define dimensionless functions ξ, θ, φ, and ψ as follows:

$$\begin{array}{c}\hfill \xi =\alpha r,\rho ={\rho }_c{\theta }^3,M_r=\frac{4}{\sqrt{\pi }}{\left(\frac{K}{G}\right)}^{3/2}\phi ,P=P_c\psi .\end{array}$$(11)

The length scale is defined as a function of the central density and pressure

$$\begin{array}{c}\hfill {\alpha }^2=\frac{\pi G{\rho }_c^2}{P_c}=\frac{\pi G{\rho }_c^{2/3}}{K},\end{array}$$(12)

where we used the fact that $P_c=K{\rho }_c{}^{4/3}$ (Eq. (9)). We notice that the mass scale is directly given by the constant K, that is to say by β_c and μ (Eq. (10)). It can be written alternatively as

$$\begin{array}{c}\hfill \frac{4}{\sqrt{\pi }}{\left(\frac{K}{G}\right)}^{3/2}=4\pi {\left(\frac{K}{\pi G}\right)}^{3/2}=\frac{4\pi {\rho }_c}{{\alpha }^3}.\end{array}$$(13)

In these dimensionless quantities, Eqs. (1) and (2) and constraint (9) read as follows:

$$\begin{array}{cc}& \frac{\mathrm{d}\psi }{\mathrm{d}\xi }=-4\frac{\phi {\theta }^3}{{\xi }^2}\hfill \end{array}$$(14)

$$\begin{array}{cc}& \frac{\mathrm{d}\phi }{\mathrm{d}\xi }={\xi }^2{\theta }^3\hfill \end{array}$$(15)

$$\begin{array}{c}\hfill \psi =\{\begin{array}{cc}{\theta }^4& \mathrm{if}\phi <{\phi }_{\mathrm{core}}\\ {\left(\frac{\phi }{{\phi }_{\mathrm{core}}}\right)}^{2/3}{\theta }^4& \mathrm{if}\phi >{\phi }_{\mathrm{core}}.\end{array}\end{array}$$(16)

We see that, in contrast to polytropes, a free parameter remains in the dimensionless structure itself, φ_core, the dimensionless mass of the convective core. Thus, exploring the full parameter space requires one to solve a series of structures, labelled by the various values of φ_core in the interval 0.46–2.02. These two limits correspond to the limit of gravitational binding and to polytropic structures, respectively, where φ_surf = φ_core (Begelman 2010). As shown in Haemmerlé (2020), these Newtonian structures allow for a precise determination of the onset point of the GR instability in non-rotating, rapidly accreting SMSs, thanks to the weakness of GR corrections (≲1%).

2.2. Rotation profiles

Because SMSs are close to the Eddington limit, the rotation of their surface is constrained by the ΩΓ-limit (Maeder & Meynet 2000), which accounts for radiation pressure, instead of the Keplerian limit, based on the simple balance between gravity and the centrifugal force. The ΩΓ-limit is extremely restrictive for SMSs due to their high Eddington factor (0.9–0.99), and it implies surface velocities of ≲10% of the Keplerian velocity (Haemmerlé et al. 2018a):

$$\begin{array}{c}\hfill \mathrm{\Omega }\lesssim 0.1{\mathrm{\Omega }}_K,\end{array}$$(17)

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_K=\sqrt{\frac{G{M}_r}{r^3}}.\end{array}$$(18)

For such velocities, the centrifugal acceleration represents less than a percent of the gravitational one:

$$\begin{array}{c}\hfill \frac{r{\mathrm{\Omega }}^2}{\left(\frac{G{M}_r}{r^2}\right)}=\frac{{\mathrm{\Omega }}^2}{{\mathrm{\Omega }}_K^2}\lesssim 0.01.\end{array}$$(19)

Interestingly, this dynamical contribution is of the same order as GR corrections and departures from the Eddington limit. Thus, similar to these effects, rotation shall play a significant role on the stability of SMSs, but it can be neglected in the definition of the hydrostatic structures.

Here, we use this advantage to define rotation profiles on the non-rotating hylotropes of Sect. 2.1. As a result of the short evolutionary timescales, angular momentum is expected to be locally conserved in the radiative envelope of SMSs (Haemmerlé et al. 2018a; Haemmerlé & Meynet 2019). Thus, we assume that the specific angular momentum of each radiative layer is given by that at accretion. Then, as the successive layers join the convective core, they advect this angular momentum, which is redistributed inside the core according to the constraint of solid rotation. We neglect the effect of the outer convective zone found in models of rapidly accreting SMSs (Hosokawa et al. 2013; Haemmerlé et al. 2018b). This envelope plays a negligible role in the GR instability as a result of its small mass and low density (Haemmerlé 2020). It plays a critical role in the rotation of the surface because its small mass covers large radii, over which angular momentum is instantaneously redistributed (Haemmerlé et al. 2018a; Haemmerlé & Meynet 2019). However, since it represents a percent of the total stellar mass, this redistribution impacts the angular momentum profile only locally once the layers contract in the radiative regions. Indeed, at each time, the specific angular momentum j₋ of the mass that leaves the envelope at its bottom r₋ is given by the angular velocity Ω_env of the envelope in solid rotation:

$$\begin{array}{c}\hfill j_{-}=\frac{2}{3}{r}_{-}^2{\mathrm{\Omega }}_{\mathrm{env}}=\frac{2}{3}{r}_{-}^2\frac{J_{\mathrm{env}}}{I_{\mathrm{env}}}=\frac{\frac{2}{3}{r}_{-}^2{M}_{\mathrm{env}}}{I_{\mathrm{env}}}⟨j_{\mathrm{accr}}⟩\end{array}$$(20)

where M_env, J_env, and I_env are the mass, the angular momentum, and the moment of inertia of the envelope, respectively, and ⟨j_accr⟩=J_env/M_env is the average specific angular momentum advected at accretion by the layers of the envelope, which provides its angular momentum content. The large density contrasts in the envelope imply that most of its mass is located near r₋, so that $I_{\mathrm{env}}\simeq \frac{2}{3}{r}_{-}^2{M}_{\mathrm{env}}$ and j₋ ≃ ⟨j_accr⟩. Since the layers of the envelope, over which the average ⟨ ⟩ is done, represent only a percent of the total mass, we can assume j_accr ≃ cst = ⟨j_accr⟩ for all these layers, so that the angular momentum j₋ advected by a layer in the radiative region is always close to what it advected at accretion: j₋ ≃ j_accr. In other words, regarding the inner rotation profiles, the convective envelope can be considered as a unique layer located in r₋. Only when we consider the surface velocity does the depth of the envelope become important.

The accreted angular momentum is chosen as a fraction f of the Keplerian momentum at the equator

$$\begin{array}{c}\hfill j_{\mathrm{accr}}=f\sqrt{GM{R}_{\mathrm{accr}}(M)},\end{array}$$(21)

where R_accr(M) is the accretion radius, that is to say the radius of the accretion shock when the star has a mass M. We use the mass-radius relation of rapidly accreting SMSs (Hosokawa et al. 2012, 2013; Haemmerlé et al. 2018b, 2019):

$$\begin{array}{c}\hfill R_{\mathrm{accr}}(M)=260R_{\odot }{\left(\frac{M}{M_{\odot }}\right)}^{1/2}.\end{array}$$(22)

This power-law, which reflects the evolution along the Eddington and Hayashi limits, reproduces the models well for masses ≲10⁵ M_⊙. Possible departures from this relation at larger masses are discussed in Sect. 4.2. The dependence R ∝ M^1/2 implies a constant surface gravity as follows:

$$\begin{array}{c}\hfill \overline{g}=\frac{\mathit{GM}}{R_{\mathrm{accr}}^2}=\frac{1}{{260}^2}\frac{G{M}_{\odot }}{R_{\odot }^2}\simeq 0.4cm{s}^{-2}.\end{array}$$(23)

The initial angular momentum of each layer M_r, which is conserved until it joins the core, is given by the angular momentum advected at the accretion of the layer (Eq. (21)), when the mass of the star was M = M_r:

$$\begin{array}{c}\hfill j(M_r)=f\sqrt{G{M}_r{R}_{\mathrm{accr}}(M_r)}.\end{array}$$(24)

With the mass-radius relation (22), it is given uniquely by the mass-coordinate M_r of the layer and the fraction f:

$$\begin{array}{c}\hfill j(M_r)=f\frac{{(G{M}_r)}^{3/4}}{{\overline{g}}^{1/4}}.\end{array}$$(25)

The angular velocity in the radiative envelope is then given by

$$\begin{array}{c}\hfill \mathrm{\Omega }(M_r)=\frac{3}{2}\frac{j(M_r)}{r^2}=\frac{3}{2}f{\mathrm{\Omega }}_K{\left(\frac{g}{\overline{g}}\right)}^{1/4},\end{array}$$(26)

where we used the Keplerian velocity (18) and the local gravitational acceleration g = GM_r/r². We see that the ratio of Ω to its Keplerian value is given uniquely by f and the local gravity. The total angular momentum of the core is the integral of the angular momentum advected by the layers M_r < M_core:

$$\begin{array}{c}\hfill J_{\mathrm{core}}={\displaystyle {\int }_0^{M_{\mathrm{core}}}j(M_r)\mathrm{d}{M}_r=\frac{4}{7}f\frac{G^{3/4}{M}_{\mathrm{core}}^{7/4}}{{\overline{g}}^{1/4}}.}\end{array}$$(27)

From the assumption of solid rotation, the angular velocity of the core is

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\mathrm{core}}=\frac{J_{\mathrm{core}}}{I_{\mathrm{core}}},\end{array}$$(28)

where I_core, the moment of inertia of the core, is given by the non-rotating structure.

In order to apply these rotation profiles to the dimensionless hylotropic structures of Sect. 2.1, we rewrote them with the use of the dimensionless functions (11). The gravitational acceleration reads as follows:

$$\begin{array}{c}\hfill g=\frac{G{M}_r}{r^2}=\frac{4\pi G{\rho }_c}{\alpha }\frac{\phi }{{\xi }^2}=\sqrt{16\pi G{P}_c}\frac{\phi }{{\xi }^2}.\end{array}$$(29)

It is fully determined by the dimensionless structure (ξ, φ) and the central pressure. At first order, the central pressure is given by radiation, and thus it depends on the central temperature only. The rotation profile in the radiative envelope is given by (Eq. (26)) the following:

$$\begin{array}{c}\hfill \frac{\mathrm{\Omega }}{{\mathrm{\Omega }}_K}=\frac{3}{2}f{\left(\frac{\sqrt{16\pi G{P}_c}}{\overline{g}}\right)}^{1/4}\frac{{\phi }^{1/4}}{{\xi }^{1/2}}.\end{array}$$(30)

From Eq. (27), the angular momentum of the core is

$$\begin{array}{c}\hfill J_{\mathrm{core}}=\frac{4}{7}f{\left(\frac{4\pi {\rho }_c}{{\alpha }^3}\right)}^{7/4}\frac{G^{3/4}{\phi }_{\mathrm{core}}^{7/4}}{{\overline{g}}^{1/4}}.\end{array}$$(31)

By definition, the moment of inertia is given by

$$\begin{array}{c}\hfill I_{\mathrm{core}}=\frac{2}{3}{\displaystyle {\int }_0^{M_{\mathrm{core}}}{r}^2\mathrm{d}{M}_r=\frac{2}{3}\frac{4\pi {\rho }_c}{{\alpha }^5}{\int }_0^{{\phi }_{\mathrm{core}}}{\xi }^2\mathrm{d}\phi .}\end{array}$$(32)

Equation (28) gives

$$\begin{array}{c}\hfill \frac{{\mathrm{\Omega }}_{\mathrm{core}}}{\sqrt{4\pi G{\rho }_c}}=\frac{6}{7}f{\left(\frac{\sqrt{16\pi G{P}_c}}{\overline{g}}\right)}^{1/4}\frac{{\phi }_{\mathrm{core}}^{7/4}}{{\int }_0^{{\phi }_{\mathrm{core}}}{\xi }^2\mathrm{d}\phi }.\end{array}$$(33)

The ratio to the Keplerian velocity (Eq. (18))

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_K=\sqrt{\frac{G{M}_r}{r^3}}=\sqrt{4\pi G{\rho }_c}\frac{{\phi }^{1/2}}{{\xi }^{3/2}}\end{array}$$(34)

is

$$\begin{array}{c}\hfill \frac{{\mathrm{\Omega }}_{\mathrm{core}}}{{\mathrm{\Omega }}_K}=\frac{6}{7}f{\left(\frac{\sqrt{16\pi G{P}_c}}{\overline{g}}\right)}^{1/4}\frac{{\phi }_{\mathrm{core}}^{7/4}}{{\int }_0^{{\phi }_{\mathrm{core}}}{\xi }^2\mathrm{d}\phi }\frac{{\xi }^{3/2}}{{\phi }^{1/2}}.\end{array}$$(35)

Thus, if we define

$$\begin{array}{c}\hfill \stackrel{\sim }{g}:=\frac{\sqrt{16\pi G{P}_c}}{\overline{g}}\end{array}$$(36)

$$\begin{array}{c}\hfill \stackrel{\sim }{\omega }:=\{\begin{array}{cc}\frac{{\xi }^{3/2}}{{\phi }^{1/2}}·\frac{\frac{4}{7}{\phi }_{\mathrm{core}}^{7/4}}{{\int }_0^{{\phi }_{\mathrm{core}}}{\xi }^{\prime 2}\mathrm{d}{\phi }^{\prime }}& \mathrm{if}\phi <{\phi }_{\mathrm{core}}\\ \frac{{\phi }^{1/4}}{{\xi }^{1/2}}& \mathrm{if}\phi >{\phi }_{\mathrm{core}}\end{array},\end{array}$$(37)

we can write the full rotation profile, that is to say the core plus envelope, in the following way:

$$\begin{array}{c}\hfill \frac{\mathrm{\Omega }}{{\mathrm{\Omega }}_K}=\frac{3}{2}f{\stackrel{\sim }{g}}^{1/4}\stackrel{\sim }{\omega }.\end{array}$$(38)

The gravity-ratio $\stackrel{\sim }{g}$ is given uniquely by the central pressure, that is by the central temperature, while the dimensionless profile $\stackrel{\sim }{\omega }$ is given uniquely by the dimensionless hylotropic structure. Thus, for a given hylotrope, the profile of Ω/Ω_K is given uniquely by the fraction f and the central temperature. We notice that for the central temperatures of SMSs (10⁷ − 10⁸ K), the gravity ratio $\stackrel{\sim }{g}$ takes values of ∼10⁶, which reflects the strong contraction of the central layers of SMSs from their large accretion radii (22). Equation (38) suggests that typical values of f ∼ 0.1−1% are required for the typical velocities of ∼0.1Ω_K imposed by the ΩΓ-limit. The Ω-profile additionally requires the scaling factor for Ω_K, given by central density (Eq. (34))

$$\begin{array}{c}\hfill \frac{\mathrm{\Omega }}{\sqrt{4\pi G{\rho }_c}}=\frac{3}{2}f{\stackrel{\sim }{g}}^{1/4}\stackrel{\sim }{\mathrm{\Omega }},\end{array}$$(39)

where

$$\begin{array}{c}\hfill \stackrel{\sim }{\mathrm{\Omega }}:=\frac{{\phi }^{1/2}}{{\xi }^{3/2}}\stackrel{\sim }{\omega }=\{\begin{array}{cc}\frac{\frac{4}{7}{\phi }_{\mathrm{core}}^{7/4}}{{\int }_0^{{\phi }_{\mathrm{core}}}{\xi }^{\prime 2}\mathrm{d}{\phi }^{\prime }}& \mathrm{if}\phi <{\phi }_{\mathrm{core}}\\ \frac{{\phi }^{3/4}}{{\xi }^2}& \mathrm{if}\phi >{\phi }_{\mathrm{core}}.\end{array}\end{array}$$(40)

From these rotation profiles, we can derive the profiles of the spin parameter:

$$\begin{array}{cc}& a_r:=\frac{c{J}_r}{G{M}_r^2}=\frac{c}{G{M}_r^2}{\displaystyle {\int }_0^{M_r}\frac{2}{3}{\mathrm{\Omega }}^{\prime }{r}^{\prime 2}\mathrm{d}{M}_r^{\prime }}\hfill \end{array}$$(41)

$$\begin{array}{cc}& =f\frac{c}{{\overline{g}}^{1/4}}{\left(4\pi G{\left(\frac{K}{\pi G}\right)}^{3/2}\right)}^{-1/4}\frac{{\int }_0^{\phi }{\stackrel{\sim }{\mathrm{\Omega }}}^{\prime }{\xi }^{\prime 2}\mathrm{d}{\phi }^{\prime }}{{\phi }^2}.\hfill \end{array}$$(42)

We can see that the scaling factor of the spin parameter is given by the mass scale and the fraction f. In the radiative envelope (M_r > M_core), the spin parameter can be directly obtained by integration of Eq. (25):

$$\begin{array}{c}\hfill a_r=f\frac{c}{{\overline{g}}^{1/4}}\frac{4}{7}{\left(G{M}_r\right)}^{-1/4}(\text{envelope}).\end{array}$$(43)

Since convection only redistributes the angular momentum in the core, the amount J_r included in the mass M_r remains unchanged, and it is set directly by the accretion law (21), independently of the hylotropic structure. The choice of the structure impacts the spin parameter only in the core. For M_r < M_core, Eq. (42) gives

$$\begin{array}{c}\hfill a_r=f\frac{c}{{\overline{g}}^{1/4}}\frac{4}{7}{\left(G{M}_r\right)}^{-1/4}·{\left(\frac{M_r}{M_{\mathrm{core}}}\right)}^{-7/4}\frac{I_r}{I_{\mathrm{core}}}(\text{core}),\end{array}$$(44)

where $I_r:=\frac{2}{3}{\int }_0^{M_r}{r}^{\prime 2}\text{d}{M}_r^{\prime }$ is the moment of inertia of the mass M_r. We see that the spin parameter is changed by a factor that depends on the two dimensionless ratios M_r/M_core and I_r/I_core, uniquely given by the dimensionless structure. Thus, the full profile of the spin parameter, in the core and the envelope, is given by the dimensionless structure, the fraction f, and the mass scale.

Finally, the rotational kinetic energy is given by (Eq. (A.25))

$$\begin{array}{c}\hfill E_{\mathrm{\Omega }}=\frac{1}{3}{\displaystyle \int }\frac{{\mathrm{\Omega }}^2}{{\mathrm{\Omega }}_K^2}\frac{G{M}_r}{r}\mathrm{d}{M}_r.\end{array}$$(45)

Its ratio to the absolute value of the gravitational energy

$$\begin{array}{c}\hfill |W|={\displaystyle \int }\frac{G{M}_r}{r}\mathrm{d}{M}_r\end{array}$$(46)

corresponds to the average of $\frac{1}{3}{\mathrm{\Omega }}^2/{\mathrm{\Omega }}_K^2$ weighted on the gravitational energy per mass unit. With the rotation profile (38), we obtain

$$\begin{array}{c}\hfill \frac{E_{\mathrm{\Omega }}}{|W|}=\frac{3}{4}{f}^2{\stackrel{\sim }{g}}^{1/2}\frac{\int {\stackrel{\sim }{\omega }}^2\frac{\phi }{\xi }\mathrm{d}\phi }{\int \frac{\phi }{\xi }\mathrm{d}\phi }=\frac{3}{4}{f}^2{\stackrel{\sim }{g}}^{1/2}\frac{\int {\stackrel{\sim }{\omega }}^2\phi {\theta }^3\xi \mathrm{d}\xi }{\int \phi {\theta }^3\xi \mathrm{d}\xi },\end{array}$$(47)

where we used Eq. (15) to substitute dξ to dφ in the integrals. Since E_Ω/|W| is an average of ${\mathrm{\Omega }}^2/{\mathrm{\Omega }}_K^2$, these two quantities share the same scaling factor. As a consequence, the profile of Ω/Ω_K (Eq. (38)) is fully determined by E_Ω/|W| and the dimensionless structure:

$$\begin{array}{c}\hfill \frac{{\mathrm{\Omega }}^2}{{\mathrm{\Omega }}_K^2}=\frac{3E_{\mathrm{\Omega }}}{|W|}\frac{{\stackrel{\sim }{\omega }}^2\int \phi {\theta }^3\xi \mathrm{d}\xi }{\int {\stackrel{\sim }{\omega }}^2\phi {\theta }^3\xi \mathrm{d}\xi }.\end{array}$$(48)

Equation (48) suggests E_Ω/|W|∼0.01 for the velocities of ∼0.1Ω_K of the ΩΓ-limit.

2.3. Impact of rotation on the GR instability

The relativistic equation for the frequency ω of adiabatic pulsations was derived by Chandrasekhar (1964) from Einstein’s equations in full generality, with the only assumption of spherical symmetry. In Haemmerlé (2020, 2021), we have shown that this equation allows one to determine, with a high precision, the stability of any spherical (i.e. non-rotating) hydrostatic structure with respect to GR, if used in the following form:

$$\begin{array}{c}\hfill \frac{{\omega }^2}{c^2}{I}_0={\displaystyle \sum _{i=1}^4{I}_i}\end{array}$$(49)

with

$$\begin{array}{cc}& I_0={\displaystyle {\int }_0^R{e}^{a+3b}(P+\rho c^2)r^4\mathrm{d}r}\hfill \end{array}$$(50)

$$\begin{array}{cc}& I_1=9{\displaystyle {\int }_0^R}{e}^{3a+b}({\mathrm{\Gamma }}_1-\frac{4}{3})P{r}^2\mathrm{d}r\hfill \end{array}$$(51)

$$\begin{array}{cc}& I_2=-12{\displaystyle {\int }_0^R}{e}^{3a+b}(a^{\prime }+\frac{b^{\prime }}{3})P{r}^3\mathrm{d}r\hfill \end{array}$$(52)

$$\begin{array}{cc}& I_3=\frac{8\pi G}{c^4}{\displaystyle {\int }_0^R}{e}^{3(a+b)}P(P+\rho c^2)r^4\mathrm{d}r\hfill \end{array}$$(53)

$$\begin{array}{cc}& I_4=-{\displaystyle {\int }_0^R}{e}^{3a+b}\frac{P^{\prime 2}}{P+\rho c^2}{r}^4\mathrm{d}r,\hfill \end{array}$$(54)

where ′ indicates the derivatives with respect to r (r = R at the surface), Γ₁ is the first adiabatic exponent, and a and b are the coefficients of the metric ds² = −e^2a(cdt)² + e^2bdr² + r²dΩ², given by the following:

$$\begin{array}{cc}& a^{\prime }=\frac{G{M}_r}{r^2{c}^2}\frac{1+\frac{4\pi r^3}{M_r{c}^2}P}{1-\frac{2G{M}_r}{r{c}^2}},e^{2a(R)}=1-\frac{2G{M}_R}{R{c}^2}\hfill \end{array}$$(55)

$$\begin{array}{cc}& e^{-2b}=1-\frac{2G{M}_r}{r{c}^2}.\hfill \end{array}$$(56)

Thanks to the weakness of GR corrections in SMSs, we can restrict our analysis to the first order post-Newtonian component of Eq. (49). Moreover, due to the proximity with the Eddington limit, Eq. (49) can be linearised in β, using the following expression for Γ₁ (Haemmerlé 2020):

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_1=\frac{4}{3}+\frac{\beta }{6}+\mathcal{O}({\beta }^2).\end{array}$$(57)

In Haemmerlé (2020), we show that the first order post-Newtonian and sub-Eddington components of integrals (51)–(54) can be written in the following form with the help of the dimensionless functions defined in Eq. (11):

$$\begin{array}{cc}& I_1=\frac{3}{2}{\beta }_c\frac{P_c}{{\alpha }^3}{\displaystyle \int }\frac{\beta }{{\beta }_c}\psi {\xi }^2\mathrm{d}\xi \hfill \end{array}$$(58)

$$\begin{array}{cc}& I_2=-16\sigma \frac{P_c}{{\alpha }^3}{\displaystyle \int }(\xi {\theta }^3+\frac{2\phi }{{\xi }^2})\psi {\xi }^3\mathrm{d}\xi \hfill \end{array}$$(59)

$$\begin{array}{cc}& I_3=8\sigma \frac{P_c}{{\alpha }^3}{\displaystyle \int }\psi {\theta }^3{\xi }^4\mathrm{d}\xi \hfill \end{array}$$(60)

$$\begin{array}{cc}& I_4=-16\sigma \frac{P_c}{{\alpha }^3}{\displaystyle \int }{\phi }^2{\theta }^3\mathrm{d}\xi .\hfill \end{array}$$(61)

The dimensionless parameter σ, defined by (Tooper 1964)

$$\begin{array}{c}\hfill \sigma :=\frac{P_c}{{\rho }_c{c}^2}=\frac{K{\rho }_c^{1/3}}{c^2}=\frac{\pi G{\rho }_c}{{\alpha }^2{c}^2},\end{array}$$(62)

represents the departures from the Newtonian limit, while β_c represents the departures from the Eddington limit. Integrals (58)–(61) have been obtained by linearisation with respect to these two parameters, that is to say neglecting terms in 𝒪(σ²), $\mathcal{O}({\beta }_c^2)$, or 𝒪(σβ_c). The compactness is related to σ by

$$\begin{array}{c}\hfill \frac{2G{M}_r}{r{c}^2}=8\sigma \frac{\phi }{\xi }.\end{array}$$(63)

Since all integrals (58)–(61) have a global σ or β_c factor in front, their content can be evaluated with Newtonian equations. In particular, Eqs. (14)–(15) and an integration by parts allow one to write

$$\begin{array}{cc}& {\displaystyle \int \psi {\theta }^3{\xi }^4\mathrm{d}\xi =\int \psi {\xi }^2\frac{\mathrm{d}\phi }{\mathrm{d}\xi }\mathrm{d}\xi =-\int \phi \frac{\mathrm{d}}{\mathrm{d}\xi }\left({\xi }^2\psi \right)\mathrm{d}\xi }\hfill \end{array}$$(64)

$$\begin{array}{cc}& =-2{\displaystyle \int \psi \phi \xi \mathrm{d}\xi +4\int {\phi }^2{\theta }^3\mathrm{d}\xi ,}\hfill \end{array}$$(65)

so that

$$\begin{array}{c}\hfill I_2=-64\sigma \frac{P_c}{{\alpha }^3}{\displaystyle \int }{\phi }^2{\theta }^3\mathrm{d}\xi =4I_4.\end{array}$$(66)

Thus, in the post-Newtonian and sub-Eddington limit, the sum (49) of the integrals reduces to

$$\begin{array}{c}\hfill \frac{2}{3}\frac{{\alpha }^3}{P_c}\frac{{\omega }^2}{c^2}{I}_0={\beta }_c{\displaystyle \int \frac{\beta }{{\beta }_c}\psi {\xi }^2\mathrm{d}\xi -\frac{32}{3}\sigma (\int \psi \phi \xi \mathrm{d}\xi +3\int {\phi }^2{\theta }^3\mathrm{d}\xi ).}\end{array}$$(67)

This sum can be re-expressed in terms of physical quantities, using Eqs. (11)–(12) and (62):

$$\begin{array}{cc}& \frac{2}{3}\frac{{\omega }^2}{c^2}{I}_0={\displaystyle \int \beta P{r}^2\mathrm{d}r-\frac{8}{3}\frac{G}{c^2}\int P{M}_{r}r\mathrm{d}r-\frac{2G^2}{c^2}\int M_r^2\rho \mathrm{d}r}\hfill \end{array}$$(68)

$$\begin{array}{cc}& ={\displaystyle \int \beta P{r}^2\mathrm{d}r-\int (\frac{2G{M}_r}{r{c}^2}+\frac{8}{3}\frac{P}{\rho c^2})\frac{G{M}_r}{r}\rho r^2\mathrm{d}r.}\hfill \end{array}$$(69)

Using the volume dV = 4πr²dr and dM_r from Eq. (2), the equation becomes

$$\begin{array}{c}\hfill {\omega }^{2}I={\displaystyle \int \beta P\mathrm{d}V-\int (\frac{2G{M}_r}{r{c}^2}+\frac{8}{3}\frac{P}{\rho c^2})\frac{G{M}_r}{r}\mathrm{d}{M}_r,}\end{array}$$(70)

where $I=\frac{8}{3}\frac{\pi }{c^2}{I}_0$ reduces to the moment of inertia in the Newtonian limit.

Equation (70) contains all the terms that are required to capture the GR instability on any hydrostatic structure in the post-Newtonian and sub-Eddington limit, in the case of pure spherical symmetry, that is without rotation. The condition for the instability is an imaginary frequency of ω² < 0. Since I is always positive, it happens when the negative term on the right-hand side, which contains the GR corrections in c⁻², exceeds, in absolute value, the positive term, which contains β, the departures from the Eddington limit. This positive term is purely Newtonian since crossed GR and β corrections would give only second-order terms. Since β is the ratio of gas to the total pressure, this integral actually corresponds to the integrated gas pressure, which is two-thirds of the total internal energy contained in the gas. On the other hand, the negative integral, which scales with the GR corrections, also scales with the total gravitational energy of the star since the integrand scales with the gravitational energy per mass unit −GM_r/r.

Since both β and the GR corrections have typical values of 0.01 in SMSs, the positive and negative integrals in Eq. (70) represent about a percent of the total internal and gravitational energies, respectively. As we have seen in Sect. 2.2, the dynamical contribution by rotation is of the same order. Thus, similar to the effect of gas pressure, the effect of rotation on the stability of the star can be reduced to its Newtonian component since relativistic rotation terms would be of second order. This method has been used by Fowler (1966), Bisnovatyi-Kogan et al. (1967), and Baumgarte & Shapiro (1999a) for polytropic models. For slow rotation near the Eddington limit, a Newtonian pulsation analysis leads to (see Appendix A, Eq. (A.21)) the following:

$$\begin{array}{c}\hfill {\omega }^{2}I={\displaystyle \int \beta P\mathrm{d}V+\frac{4}{9}\int \frac{{\mathrm{\Omega }}^2}{{\mathrm{\Omega }}_K^2}\frac{G{M}_r}{r}\mathrm{d}{M}_r.}\end{array}$$(71)

The Ω²-term in Eq. (71) is the rotational energy (45) multiplied by a factor of 4/3. Since this term is positive, we see that rotation acts as a stabilising agent, similar to gas pressure, against the destabilising GR effects. It is expressed by adding this term in the post-Newtonian Eq. (70):

$$\begin{array}{c}\hfill {\omega }^{2}I={\displaystyle \int \beta P\mathrm{d}V-\int (\frac{2G{M}_r}{r{c}^2}+\frac{8}{3}\frac{P}{\rho c^2}-\frac{4}{9}\frac{{\mathrm{\Omega }}^2}{{\mathrm{\Omega }}_K^2})\frac{G{M}_r}{r}\mathrm{d}{M}_r.}\end{array}$$(72)

In the sum of the I_i that appears in Eq. (49), this term translates into a fifth integral:

$$\begin{array}{c}\hfill I_5:=\frac{2}{3}{\displaystyle \int {\mathrm{\Omega }}^2\rho r^4\mathrm{d}r=\frac{8}{3}\frac{P_c}{{\alpha }^3}\int \frac{{\mathrm{\Omega }}^2}{{\mathrm{\Omega }}_K^2}\phi {\theta }^3\xi \mathrm{d}\xi .}\end{array}$$(73)

With the rotation profiles defined in Sect. 2.2 (Eq. (38)), it gives

$$\begin{array}{c}\hfill I_5=6f^2{\stackrel{\sim }{g}}^{1/2}\frac{P_c}{{\alpha }^3}{\displaystyle \int {\stackrel{\sim }{\omega }}^2\phi {\theta }^3\xi \mathrm{d}\xi .}\end{array}$$(74)

Adding this term to Eq. (67), we obtain

$$\begin{array}{cc}& \frac{2}{3}\frac{{\alpha }^3}{P_c}\frac{{\omega }^2}{c^2}{I}_0={\beta }_c{\displaystyle \int \frac{\beta }{{\beta }_c}\psi {\xi }^2\mathrm{d}\xi +4f^2{\stackrel{\sim }{g}}^{1/2}\int {\stackrel{\sim }{\omega }}^2\phi {\theta }^3\xi \mathrm{d}\xi }\hfill \\ \hfill & -\frac{32}{3}\sigma {\displaystyle (\int \psi \phi \xi \mathrm{d}\xi +3\int {\phi }^2{\theta }^3\mathrm{d}\xi ).}\hfill \end{array}$$(75)

We see that the rotation term scales with the independent fraction f², and with the square root of the gravity ratio $\stackrel{\sim }{g}$. From its definition (36), ${\stackrel{\sim }{g}}^{1/2}$ scales with $P_c^{1/4}$, that is with the central temperature:

$$\begin{array}{c}\hfill {\stackrel{\sim }{g}}^{1/2}(T_c)={\left(\frac{\sqrt{16\pi G}}{\overline{g}}\right)}^{1/2}{\left(\frac{a_{\mathrm{SB}}}{3}\right)}^{1/4}{T}_c.\end{array}$$(76)

The product of the other two scaling factors, β_c and σ, is also given by the central temperature and the chemical composition:

$$\begin{array}{c}\hfill {\beta }_c\sigma =\frac{k_{\mathrm{B}}{T}_c}{\mathrm{\mu }{m}_{\mathrm{H}}{c}^2}\end{array}$$(77)

by definitions (6)–(62). Thus, the limit of stability, at which the right-hand side of Eq. (75) vanishes, is given by a second-order polynomial equation in β_c:

$$\begin{array}{c}\hfill {\beta }_{c,\mathrm{crit}}^2+4f^2{\stackrel{\sim }{g}}^{1/2}(T_c){\stackrel{\sim }{I}}_1·{\beta }_{c,\mathrm{crit}}-\frac{k_{\mathrm{B}}{T}_c}{\mathrm{\mu }{m}_{\mathrm{H}}{c}^2}{\stackrel{\sim }{I}}_0=0,\end{array}$$(78)

where the integral ratios

$$\begin{array}{cc}& {\stackrel{\sim }{I}}_0:=\frac{32}{3}\frac{\int \psi \phi \xi \mathrm{d}\xi +3\int {\phi }^2{\theta }^3\mathrm{d}\xi }{\int \frac{\beta }{{\beta }_c}\psi {\xi }^2\mathrm{d}\xi }\hfill \end{array}$$(79)

$$\begin{array}{cc}& {\stackrel{\sim }{I}}_1:=\frac{\int {\stackrel{\sim }{\omega }}^2\phi {\theta }^3\xi \mathrm{d}\xi }{\int \frac{\beta }{{\beta }_c}\psi {\xi }^2\mathrm{d}\xi }\hfill \end{array}$$(80)

are uniquely given by the dimensionless structure $(\xi ,\theta ,\phi ,\psi ,\stackrel{\sim }{\omega })$. The ratio ${\stackrel{\sim }{I}}_0$ corresponds to the critical ratio β_c/σ in the non-rotating case (Haemmerlé 2020), while ${\stackrel{\sim }{I}}_1$ gives the correction from rotation. Equation (78) has a unique positive root of

$$\begin{array}{c}\hfill {\beta }_{c,\mathrm{crit}}=\sqrt{\frac{k_{\mathrm{B}}{T}_c}{\mathrm{\mu }{m}_{\mathrm{H}}{c}^2}{\stackrel{\sim }{I}}_0+{\left(2f^2{\stackrel{\sim }{g}}^{1/2}(T_c){\stackrel{\sim }{I}}_1\right)}^2}-2f^2{\stackrel{\sim }{g}}^{1/2}(T_c){\stackrel{\sim }{I}}_1.\end{array}$$(81)

From Eq. (10), the corresponding mass scale (13) is directly given by ${\beta }_{c,\mathrm{crit}}$ and the chemical composition μ:

$$\begin{array}{c}\hfill 4\pi {\left(\frac{K_{\mathrm{crit}}}{\pi G}\right)}^{3/2}=\frac{4}{{\pi }^{1/2}{G}^{3/2}}{\left(\frac{3}{a_{\mathrm{SB}}}\right)}^{1/2}{\left(\frac{k_{\mathrm{B}}}{\mathrm{\mu }{m}_{\mathrm{H}}}\right)}^2{\beta }_{c,\mathrm{crit}}^{-2}.\end{array}$$(82)

Thus, the critical masses for the star and its core

$$\begin{array}{cc}& M_{\mathrm{crit}}=4\pi {\left(\frac{K_{\mathrm{crit}}}{\pi G}\right)}^{3/2}{\phi }_{\mathrm{surf}}\hfill \end{array}$$(83)

$$\begin{array}{cc}& M_{\mathrm{core},\mathrm{crit}}=4\pi {\left(\frac{K_{\mathrm{crit}}}{\pi G}\right)}^{3/2}{\phi }_{\mathrm{core}}\hfill \end{array}$$(84)

are uniquely determined by the dimensionless structure $(\xi ,\theta ,\phi ,\psi ,\stackrel{\sim }{\omega })$ and the choice of (f, T_c, μ).

In the case where f = 0, when the stability limit is set by gas only, Eq. (81) leads to ${\beta }_c/\sigma ={\stackrel{\sim }{I}}_0$. In the opposite limit, when the stability is dominated by the effect of rotation,

$$\begin{array}{c}\hfill \sqrt{\frac{k_{\mathrm{B}}{T}_c}{\mathrm{\mu }{m}_{\mathrm{H}}{c}^2}{\stackrel{\sim }{I}}_0}\ll 2f^2{\stackrel{\sim }{g}}^{1/2}(T_c){\stackrel{\sim }{I}}_1,\end{array}$$(85)

a linear development of Eq. (81) gives

$$\begin{array}{c}\hfill {\beta }_{c,\mathrm{crit}}\simeq \frac{1}{4f^2}{\left(\frac{\overline{g}}{\sqrt{16\pi G}}\right)}^{1/2}{\left(\frac{a_{\mathrm{SB}}}{3}\right)}^{-1/4}\frac{k_{\mathrm{B}}}{\mathrm{\mu }{m}_{\mathrm{H}}{c}^2}\frac{{\stackrel{\sim }{I}}_0}{{\stackrel{\sim }{I}}_1}.\end{array}$$(86)

Expression (86) of β_c,crit is independent of T_c, which implies that for large f, the critical mass scale (82) is independent of T_c and μ:

$$\begin{array}{c}\hfill 4\pi {\left(\frac{K_{\mathrm{crit}}}{\pi G}\right)}^{3/2}\simeq f^4\frac{256}{\overline{g}}\frac{c^4}{G}\frac{{\stackrel{\sim }{I}}_1^2}{{\stackrel{\sim }{I}}_0^2}.\end{array}$$(87)

In other words, in the limit of large f, the critical mass of a given hylotrope is set uniquely by f. This fraction appears to the power 4 in Eq. (87), which shows the strong dependence of the critical masses on the accreted angular momentum. Moreover, it implies that for large f, the profile of the spin parameter (Eq. (42)) at the limit of stability is independent of f, and is given uniquely by the following dimensionless structure:

$$\begin{array}{c}\hfill a_r=\frac{1}{4}\sqrt{\frac{{\stackrel{\sim }{I}}_0}{{\stackrel{\sim }{I}}_1}}\frac{{\int }_0^{\phi }{\stackrel{\sim }{\mathrm{\Omega }}}^{\prime }{\xi }^{\prime 2}\mathrm{d}{\phi }^{\prime }}{{\phi }^2}.\end{array}$$(88)

In particular, at the surface, there is the following:

$$\begin{array}{c}\hfill a_R=\frac{1}{7}\sqrt{\frac{{\stackrel{\sim }{I}}_0}{{\stackrel{\sim }{I}}_1}}{\phi }_{\mathrm{surf}}^{-1/4}.\end{array}$$(89)

We emphasise, however, that this universality is broken when gas plays a significant role in the determination of the stability limit.

Finally, we notice that the critical condition (81) can be expressed as a function of E_Ω/|W| (Eq. (47)) instead of f:

$$\begin{array}{c}\hfill {\beta }_{c,\mathrm{crit}}=\sqrt{\frac{k_{\mathrm{B}}{T}_c}{\mathrm{\mu }{m}_{\mathrm{H}}{c}^2}{\stackrel{\sim }{I}}_0+{\left(2{\stackrel{\sim }{I}}_2\frac{E_{\mathrm{\Omega }}}{|W|}\right)}^2}-2{\stackrel{\sim }{I}}_2\frac{E_{\mathrm{\Omega }}}{|W|}\end{array}$$(90)

with

$$\begin{array}{c}\hfill {\stackrel{\sim }{I}}_2:=\frac{4}{3}\frac{\int \phi {\theta }^3\xi \mathrm{d}\xi }{\int \frac{\beta }{{\beta }_c}\psi {\xi }^2\mathrm{d}\xi }=\frac{\int \psi {\xi }^2\mathrm{d}\xi }{\int \frac{\beta }{{\beta }_c}\psi {\xi }^2\mathrm{d}\xi }.\end{array}$$(91)

The integral at the numerator has been transformed by an integration by parts, using Eq. (14).

3. Results

3.1. Rotation profiles

The dimensionless rotation profiles of Eq. (40) are shown in Fig. 1 for a choice of hylotropes with the indicated values of φ_core. The convective core, where solid rotation is imposed, clearly appears in each profile. In particular, the case φ_core = 2.02, which corresonds to the polytropic limit, is fully rotating as a solid body. In contrast, strong differential rotation appears in the envelope of the other models; and for φ_core ≲ 1, the angular velocity in the core exceeds, by several orders of magnitude, that of the surface, consistently with the models accounting for full stellar evolution (Haemmerlé et al. 2018a; Haemmerlé & Meynet 2019). We also see that the instantaneous angular momentum transport in the core introduces a discontinuity in the $\stackrel{\sim }{\mathrm{\Omega }}$-profiles at the interface with the envelope. Physically, such infinite gradients shall be smoothed out by shears, but here we can assume that such shears only lead to a local redistribution of the angular momentum, without a significant impact on the rotation profiles.

Fig. 1. Dimensionless rotation profiles of Eq. (40) for a set of hylotropes with the indicated values of φcore.

From these dimensionless rotation profiles, the profiles of Ω/Ω_K (Eq. (38)) are given by the ratio E_Ω/|W| (Eq. (48)). They are shown in Fig. 2 for E_Ω/|W| = 0.01. We see that the maximum values of Ω/Ω_K are always reached at the surface of the convective core. The rotation velocity typically remains at 0.1–0.2 Ω_K, except in the very outer layers of the polytropic model (φ_core = 2.02). In this case, the ΩΓ-limit shall be exceeded at the surface, which shows that accretion onto fully convective SMSs requires E_Ω/|W|≪0.01. On the other hand, the models with φ_core ≲ 2 keep a surface velocity of ∼0.1Ω_K, which is consistent with the ΩΓ-limit. We see that for a given constraint on the surface velocity, hylotropic structures allow for larger E_Ω/|W|, thanks to differential rotation in the envelope. However, as already noticed in Sect. 2.2, the surface rotation of rapidly accreting SMSs depends on the depth of their convective envelope, which is not included in the hylotropic structures. The constraint on f that arises from this effect is discussed in Sect. 4.1.

Fig. 2. Profiles of Ω/ΩK (Eq. (48)) for EΩ/|W| = 0.01, for a set of hylotropes with the indicated values of φcore.

3.2. Critical masses

The critical masses obtained with Eqs. (83)–(84) for the full series of hylotropes are shown in Fig. 3 for T_c = 1.8 × 10⁸ K, μ = 0.6, and a series of values of the fraction f. These values of T_c and μ are relevant for Pop III SMSs (Hosokawa et al. 2013; Umeda et al. 2016; Woods et al. 2017; Haemmerlé et al. 2018b, 2019). In Fig. 4, the critical masses are shown for the lower T_c = 8.5 × 10⁷ K relevant for Pop I SMSs (Haemmerlé et al. 2019). For these central temperatures, the central densities are typically a few g cm⁻³ for M ∼ 10⁵ − 10⁶ M_⊙, down to ∼0.01−0.1 g cm⁻³ for M ∼ 10⁸ − 10⁹ M_⊙.

Fig. 3. Limits of stability for Pop III SMSs accreting angular momentum at a fraction f of the Keplerian momentum. The solid red lines show the hylotropic limits for the indicated f (Eq. (81)), with Tc = 1.8 × 108 K and μ = 0.6, relevant for Pop III SMSs. The dotted red lines indicate the limit of stability for constant ratios 0.001–0.01–0.1 of rotational to gravitational energy (Eq. (90)). The coloured areas reflect the relative contribution by gas and rotation in the determination of the stability limit (81): From yellow to red, the rotation term $2f^2{\stackrel{\sim }{g}}^{1/2}(T_c){\stackrel{\sim }{I}}_1$ represents < 0.1, 0.1−1, 1−10, and > 10 times the gas term ${(\frac{k_{\text{B}}{T}_c}{\mathrm{\mu }{m}_{\text{H}}{c}^2}{\stackrel{\sim }{I}}_0)}^{1/2}$, respectively. The Pop III GENEC tracks of Haemmerlé et al. (2018b, 2019) are shown for the indicated accretion rates. The mass fraction of the convective core is indicated by the grey diagonals.

As in the non-rotating case (Haemmerlé 2020), for every given f, the maximum total mass consistent with GR stability is a decreasing function of the core mass fraction. We see that the effect of rotation on the stability of the star becomes significant as soon as f ≳ 0.1%. For f = 0.2−0.3%, the mass limit is shifted by a significant factor, and for f = 0.3−0.5% by an order of magnitude. In this last case, the GR instability requires masses ≳10⁶ M_⊙. For f = 1%, the GR instability cannot be reached if the mass does not exceed several 10⁷ M_⊙, and the star can be stable with masses ∼10⁹ M_⊙ if the core does not exceed 10⁶ M_⊙. In the range of f = 1.5−2%, the star remains stable up to masses of 10⁸ − 10⁹ M_⊙ for any core mass. We notice the small knee of the mass limits near the polytropic structures (M_core ≳ 75%M). The stabilising effect of rotation is slightly decreased by the solid rotation of the core, which removes angular momentum from the densest central regions that are dominant for the instability.

The relative contribution by gas and rotation in the determination of the stability limit is made visible by coloured areas. These contributions are estimated by the two terms that appear in the critical condition (81): ${\left(\frac{k_{\text{B}}{T}_c}{\mathrm{\mu }{m}_{\text{H}}{c}^2}{\stackrel{\sim }{I}}_0\right)}^{1/2}$, which represents gas, and $2f^2{\stackrel{\sim }{g}}^{1/2}(T_c){\stackrel{\sim }{I}}_1$, which represents rotation. The four regions, from yellow to red, correspond to the cases where the rotation term represents < 0.1, 0.1−1, 1−10, and > 10 times the gas term. We see that for 0.1%≲f ≲ 1%, both gas and rotation play a role in the determination of the stability limit. But for f ≳ 1%, rotation becomes the only significant stabilising effect (condition 85). In this case, as noticed in Sect. 2.3, the critical masses of the various hylotropes are uniquely given by f, independently of the thermal properties (Eq. (87)). As a consequence, the effect of a change in the central temperature, between Figs. 3 and 4, impacts the mass limits only for f ≲ 1%, while the limits for f ≳ 1% remain unchanged.

Fig. 4. Same as Fig. 3, but for the Pop I case (Tc = 8.5 × 107 K).

The mass limits obtained for constant E_Ω/|W| (Eq. (90)) are shown in Figs. 3 and 4 as red dotted lines for E_Ω/|W| = 0.001−0.01−0.1. Since the ratio E_Ω/|W| scales with f² (Eq. (47)), the range of one order of magnitude in f ∼ 0.1−1% corresponds to a range of two orders of magnitude in E_Ω/|W|∼0.001−0.1. For M_core > 10%M, the ratio E_Ω/|W| = 0.01 corresponds to f ∼ 0.3−0.4%. The fastest rotators shown in these figures (f ∼ 1%) have E_Ω/|W|∼0.1, which corresponds to typical velocities ∼0.3Ω_K (Eq. (48)). We see that for a given value of E_Ω/|W|, hylotropic structures with small M_core/M remain stable up to larger masses than polytropes (M_core = M). More precisely, for a given relative contribution E_Ω/|W|, rotation has a stronger stabilising effect on hylotropes with M_core < M than on polytropes. It reflects the fact that rotational energy is more centralised when the core is smaller (Fig. 2), and thus it has a stronger impact on the dense regions that are relevent for the GR instability.

3.3. Spin parameter

The profile of the spin parameter is given by Eq. (42). For each hylotrope, it is fully determined by the fraction f and the mass scale. As is noticeable in Sect. 2.2, f and the mass scale appear in Eq. (42) as global factors only, so that a change in these quantities affects only the scale of the spin parameter, while the shape of the profile is unique for each hylotrope. Here, we are interested in the profiles of the spin parameter at the limit of stability, where the mass scale is fixed by f, T_c, and μ Eqs. (81) and (82). Since the critical mass scale is not the same for the various hylotropes, the profiles of the spin parameter at the critical limit are rescaled by a different factor when f, T_c, and μ are changed.

The critical profiles obtained for f = 0.2%, T_c = 1.5 × 10⁸ K, and μ = 0.6 are shown in Fig. 5 for three hylotropes in the series. The case f = 1%, T_c = 9 × 10⁷ K, and μ = 0.6 is shown in Fig. 6. The profiles are compared with that of the angular momentum accretion law (25), expressed in Eq. (43). By construction, the profiles match each other in the radiative envelopes, and they are impacted by the stellar structure only in the convective core. By definition, the spin parameter diverges at the centre, but we see that for all the profiles, it decreases below unity at a small fraction of the total mass. The profile set by the accretion law, and kept in the envelope, is monotonously decreasing. Thus, if the spin parameter goes below unity in the deep regions of the envelope, it necessarily remains so in the whole envelope. Only in the convective core, where solid rotation redistributes the angular momentum, can the spin parameter increase outwards, which is in agreement with the polytropic models (Baumgarte & Shapiro 1999a). However, it occurs only for hylotropes with φ_core ≳ 1, that is for a core mass fraction ≳40%.

Fig. 5. Internal profile of the dimensionless spin parameter at the limit of stability for f = 0.2%, Tc = 1.5 × 108 K, and μ = 0.6, as a function of the mass-coordinate and for a choice of hylotropes with indicated φcore. The surface of each hylotrope is indicated by a white square. The black dashed curve shows the angular momentum accretion law (43). The red dashed curve and the red square indicate the limits for black hole formation, according to the criterions discussed in Sect. 4.3.

We notice that in these two cases, f = 0.2% and 1%, the stability limit is set predominantly by gas and by rotation, respectively (Figs. 3 and 4). Condition (85) is already satisfied for f = 1%, so that the critical profiles shown in Fig. 6 actually correspond to the ‘universal’ profiles of Eq. (88), independent of f and given uniquely by the dimensionless structure. As a consequence, these profiles represent an upper limit for each hylotrope, and the increase of the critical masses for f > 1% only translates into a rescaling along the x-axis. In particular, in the case of the polytrope (φ_core = 2.02), we refind the profile of Baumgarte & Shapiro (1999a), with the ‘universal’ value 0.87 at the surface, given by Eq. (89).

Fig. 6. Same as Fig. 5, but for f = 1% and Tc = 9 × 107 K.

4. Discussion

4.1. Final mass as a function of the accretion history

Figures 3 and 4 show that the maximum masses consistent with GR stability sensitively depend on the fraction f of accreted angular momentum, already in the range 0.1–1%. The actual value of f at which a SMS can accrete is set by the ΩΓ-limit, which imposes an upper limit for its surface velocity. For stars with masses ≳10⁵ M_⊙, the surface velocity cannot exceed ∼10% of the Keplerian velocity (Haemmerlé et al. 2018a; Haemmerlé & Meynet 2019). An extrapolation of the models suggests that the limit goes down to ∼3% for ≳10⁷ M_⊙, and to ∼1% for ≳10⁹ M_⊙. As noticed in Sect. 2.2, for a given fraction f of accreted angular momentum, the surface velocities of rapidly accreting SMSs depend on the depth of their convective envelopes, which are not included in the hylotropic models. Models accounting for full stellar evolution show that this envelope covers more than two-thirds of the total stellar radius for rates ≲1 M_⊙ yr⁻¹, which implies f ≲ 0.2−0.3% (Haemmerlé & Meynet 2019). On the other hand, the depth of the envelope is significantly reduced for rates ≳100 M_⊙ yr⁻¹, and it covers a small fraction of the total radius. As long as it remains the case, the SMS can accrete angular momentum at f ∼ 1% up to ∼10⁹ M_⊙ without facing the ΩΓ-limit.

Rapidly accreting SMSs are invoked in two different versions of direct collapse: atomically cooled haloes and galaxy mergers. The accretion rates in atomically cooled haloes are expected to be 0.1−10 M_⊙ yr⁻¹ and most probably ≲1 M_⊙ yr⁻¹ (Latif et al. 2013; Chon et al. 2018; Patrick et al. 2020), which imposes f ≲ 0.2−0.3%. On the other hand, galaxy mergers allow for rates ≳100 M_⊙ yr⁻¹ (Mayer et al. 2010, 2015; Mayer & Bonoli 2019), for which f ≃ 1% remains consistent with the ΩΓ-limit up to large masses. Another difference between these two scenarios is the chemical composition, which is primordial in atomically cooled haloes, while solar metallicities are expected in the case of galaxy mergers. Black hole formation at the collapse of a Pop I SMS requires a mass ≳10⁶ M_⊙ (Montero et al. 2012). In the non-rotating case, such final masses could be reached only for rates ≳1000 M_⊙ yr⁻¹ (Haemmerlé 2020, 2021).

The evolutionary tracks of the GENEC models accounting for full stellar evolution (Haemmerlé et al. 2018b, 2019), computed under accretion at 0.1–1000 M_⊙ yr⁻¹ up to masses of 10⁵ − 10⁶ M_⊙, are plotted in Figs. 3 and 4. These models were built on the assumption of hydrostatic equilibrium, and thus they are insensitive to dynamical instabilities. Their stability with respect to GR perturbations has been addressed in Haemmerlé (2021) for the non-rotating case by the direct use of Eq. (49). As is shown in Haemmerlé (2020), the hylotropic limit reproduces these final masses with a high precision for accretion rates ≥10 M_⊙ yr⁻¹ and remains a good approximation for lower rates. We see that none of these tracks reach the stability limits for f ≳ 0.2%. An extrapolation of the tracks suggests that for conditions of atomically cooled haloes (Pop III, Ṁ ∼ 0.1 − 1 M_⊙ yr⁻¹, f ∼ 0.2−0.3%), rotation allows one to increase the final mass by a factor of ∼2 only. It implies that the mass of SMSs forming in this scenario always remain < 10⁶ M_⊙, and more probably ≲5 × 10⁵ M_⊙. On the other hand, for the conditions of metal-rich galaxy mergers (Pop I, Ṁ ≳ 100 M_⊙ yr⁻¹, f ∼ 1%), rotation could allow SMSs to increase their final mass by several orders of magnitude, leading to an upper mass limit of 10⁸ − 10⁹ M_⊙. Overall, as in the non-rotating case, the final masses of SMSs range in distinct intervals in the two versions of direct collapse: 10⁵ − 10⁶ M_⊙ for atomically cooled haloes and 10⁶ − 10⁹ M_⊙ in galaxy mergers. We notice also that the minimum rate required for final masses ≳10⁶ M_⊙ is slightly reduced by rotation, since for 100 M_⊙ yr⁻¹ already, we can expect f ∼ 1%.

Interestingly, an extrapolation of the tracks at Ṁ ≳ 100 M_⊙ yr⁻¹ in Figs. 3 and 4 suggests that stars accreting at these rates are mostly convective when they reach the GR instability for f ≃ 1%. The same is true for Ṁ ≲ 0.1 M_⊙ yr⁻¹ and f = 0.2%. If the star becomes fully convective before the instability, the instantaneous angular momentum transport affects the surface velocity, and the values of f consistent with the ΩΓ-limit decrease below < 0.1. Thus, either the accreted angular momentum decreases, so that the star rapidly becomes GR unstable, or the surface reaches break-up. In this case, accretion must stop and the star evolves along the mass-shedding limit, in a similar pathway as monolithic models (Fowler 1966; Bisnovatyi-Kogan et al. 1967; Baumgarte & Shapiro 1999a). In this case, only second-order post-Newtonian corrections have been found to trigger the collapse (Bisnovatyi-Kogan et al. 1967). We emphasise, however, that in the accretion scenario considered here, these second-order corrections always remain negligible, and we checked numerically that they do not play a role in the determination of the limits in Figs. 3 and 4. It is also important to notice that the age of a star accreting at 100–1000 M_⊙ yr⁻¹ when it reaches 10⁸ − 10⁹ M_⊙ is ∼10⁶ years, at which H-burning might be achieved (Umeda et al. 2016; Woods et al. 2021). In this case, pair-instability might be reached before the GR instability (Woods et al. 2020).

4.2. Angular momentum accretion law and rotation profiles

The definition of the rotation profiles in Sect. 2.2 relies on the angular momentum accretion law (21) and the mass-radius relation (22). This relation is verified by stellar evolution models for masses ≲10⁵ M_⊙ (Hosokawa et al. 2013; Haemmerlé et al. 2018b), but for larger masses the models depart towards smaller radii. In the absence of stellar evolution models accounting for accretion up to masses ≥10⁶ M_⊙, the mass-radius relation is not known for the most massive objects considered in the present work. According to Eq. (24), smaller radii would imply a slower rotation. However, the radius appears in this equation through its square root, so that a decrease by an order of magnitude in radius would be equivalent to a decrease by a factor of a few only in f. Moreover, the mass-radius relation of the zero-age main-sequence (ZAMS) imposes a lower limit to the radius. Similar to the Hayashi line, the ZAMS of SMSs follows a nearly constant effective temperature (Woods et al. 2020), leading to the same power-law R ∝ M^1/2 as Eq. (22), rescaled down by 2–3 orders of magnitude in radius. Such a shift would be equivalent to a decrease in f by one order of magnitude. Since departures from relation (22) might appear only for M ≳ 10⁶ M_⊙, we can estimate that this effect would translate, in the (M_core, M) diagram of Figs. 3 and 4, into a shift of the mass limits for f ≳ 0.5%, and down to f = 0.4%. We notice however that more recent models of very-massive stars show a shift in the ZAMS towards lower effective temperatures and larger radii (Gräfener 2021). This effect might inhibit a possible contraction of SMSs in the mass range ≳10⁶ M_⊙, allowing for accretion of relatively high angular momentum.

In addition to the mass-radius relation, the definition of the rotation profiles in Sect. 2.2 relies on the assumption of local angular momentum conservation in the envelope, which leads to strong differential rotation. Stellar models (Haemmerlé et al. 2018a) have shown shear diffusion and meridional circulation to be inefficient in rapidly accreting SMSs because of the short timescales and slow rotations. On the other hand, magnetic fields (e.g., a Taylor-Spruit dynamo), which are neglected in the present work, could inhibit significantly differential rotation. In the extreme case of solid rotation, the constraint from the ΩΓ-limit at the surface translates into f < 0.1%, and the effect of rotation on the stability becomes negligible. In other words, for given surface velocities, solid rotation would impose a low angular momentum in the deep regions of the star that are decisive for the GR instability, and this could bring back the stability limit to the non-rotating case.

We also notice that the ratio Ω/Ω_K remains < 0.7 for all the rotation profiles considered (f ≲ 2%), except in the very outer layers of the nearly-polytropic models (M_core/M > 0.95), which play a negligible role in the GR instability. No significant flattening is expected for such slow rotations (Ekström et al. 2008), which justifies the use of the spherical hylotropic structures. Moreover, in the deep regions that are relevent for the GR instability, the ratio $\frac{1}{3}{\mathrm{\Omega }}^2/{\mathrm{\Omega }}_K^2$ typically remains between 0.1–0.2. Thus, we can estimate that post-Newtonian rotation terms in the pulsation Eq. (72), which are neglected in the present method, should remain small compared to the other terms.

4.3. Spin parameter and black hole formation

The spin parameter (41) is key for the direct formation of a supermassive black hole during the collapse of a SMS. Without angular momentum transport, the profile a_r(M_r) at the onset of the instability is conserved during the collapse. For monolithic SMSs, evolving along the mass-shedding limit and fully stabilised by rotation (condition (85)), analytical models in the post-Newtonian limit feature a universal profile of the spin parameter, with a unique value of 0.87 at the surface in particular (Baumgarte & Shapiro 1999a). Except in the central ∼10% of the mass, the spin remains always a_r < 1, the maximum value allowed for the formation of a Kerr black hole, which suggests that the whole stellar mass can be accreted by the black hole in a dynamical time. We notice that since monolithic models are rotating at the Keplerian limit, the conditions for a Newtonian treatment of rotation are not satisfied, and the full GR numerical simulations of Baumgarte & Shapiro (1999a) lead to larger values for the spin parameter (0.97). Only the slow rotations imposed by accretion along the ΩΓ-limit guarantee the validity of the Newtonian treatment of rotation against GR instability.

The hydrodynamical simulations of the collapse of monolithic SMSs have shown that the outer ∼10% of the stellar mass remains in orbit outside the horizon (Shibata & Shapiro 2002; Shapiro & Shibata 2002; Liu et al. 2007; Uchida et al. 2017; Sun et al. 2017, 2018). The resulting non-axisymmetric structure is found to trigger gravitational wave emission and ultra-long gamma-ray bursts, which are currently the main observational signatures of the existence of SMSs proposed in the literature for future detection. Shapiro & Shibata (2002) showed that the mass of the orbiting structure in the simulations can be derived analytically by comparing the specific angular momentum of each layer M_r with that of the innermost stable circular orbit (ISCO) of a black hole with mass M_r and spin a_r: if the angular momentum exceeds the ISCO value in a layer, the layer remains in orbit outside the horizon.

The critical profiles of the spin parameter shown in Figs. 5 and 6 correspond to typical conditions of atomically cooled haloes and galaxy-mergers, respectively. For conditions of galaxy mergers, we refind the profile of Baumgarte & Shapiro (1999a) for the polytrope (φ_core = 2.02), which reflects the fact that condition (85) is already satisfied for f = 1%. In contrast, for the case of atomically cooled haloes (f = 0.2%, Fig. 5), gas plays the dominant role in setting the stability limit (Fig. 3), and the spin parameter departs from the ‘universal’ profiles, as already noticed by Shibata et al. (2016b) and Butler et al. (2018) for monolithic models. By comparing the angular momentum profiles with the ISCO values of the various layers (Eq. (17) of Shapiro & Shibata 2002), we can estimate in which cases an outer envelope rotates fast enough to remain in circular orbit after black hole formation. The limit of this envelope is shown by a red square in the profiles, and we see that it only occurs for nearly polytropic models. For f = 1%, it requires a core mass fraction ≳80% and in the polytropic limit we refind that the outer 10% of the stellar mass can remain in a circular orbit. For f = 0.2%, the core mass fraction must exceed ≳99% and even in the polytropic limit it is only the outer ∼1% of the mass that can stand in orbit.

These results indicate that the angular momentum barrier never prevents the direct formation of a supermassive black hole at the collape of a SMS. They suggest that the star must be mostly convective for a significant mass fraction to remain outside the horizon at the formation of the black hole. Moreover, the mass fraction of such outer structures are found much larger under the conditions of galaxy mergers than under those of atomically cooled haloes. Thus, the conditions of galaxy mergers appear to be more promising than those of atomically cooled haloes for detectable gravitational wave emission and ultra-long gamma-ray bursts during the collapse.

5. Summary and conclusions

We have extended our previous works on the GR instability in rapidly accreting SMSs (Haemmerlé 2020, 2021), based on the relativistic pulsation equation of Chandrasekhar (1964), by including rotation according to the method of Fowler (1966). On the basis of hylotropic models relevent for rapidly accreting SMSs (Begelman 2010; Haemmerlé et al. 2019; Haemmerlé 2020), we defined rotation profiles by assuming local angular momentum conservation in radiative regions, which allows for strong differential rotation. The angular momentum advected by each layer at accretion is assumed to represent a fraction f of the Keplerian angular momentum.

The impact of rotation on the GR instability is captured by a positive term in the pulsation equation, which translates into a stabilising effect, such as gas pressure, against the destabilising GR effects. The stabilising effect of rotation becomes significant as soon as f ≳ 0.1%. For f ∼ 0.2%, the maximum masses consistent with GR stability are increased by a significant factor compared to the non-rotating case. For f ∼ 1%, it is increased by several orders of magnitude, and the GR instability cannot be reached if the mass does not exceed 10⁷ − 10⁸ M_⊙.

The relevent values of f depend on the channel of direct collapse: While f ∼ 0.2% is expected in atomically cooled haloes, because of the deep convective envelope of SMSs accreting at rates ≲1 M_⊙ yr⁻¹, values of f ∼ 1% appear consistent in conditions of galaxy mergers, thanks to the larger rates. It indicates that rotation allows for SMSs forming in atomically cooled haloes to increase their final masses by a factor of a few only. In contrast, SMSs forming in galaxy mergers could increase their masses by several orders of magnitude compared to the non-rotating case, possibly up to ∼10⁹ M_⊙. It follows that the final masses of rapidly accreting SMSs range in distinct intervals in the different versions of direct collapse: 10⁵ − 10⁶ M_⊙ for primordial, atomically cooled haloes and 10⁶ − 10⁹ M_⊙ for metal-rich galaxy mergers.

The rotational properties of rapidly accreting SMSs at the onset of GR instability are not universal, in contrast to the monolithic, non-accreting case. The monolithic profiles are found in the limit of fully convective SMSs with f ≳ 1%, which illustrates the fact that, already for such low f, rotation plays the dominant role in the stability of the star. We estimate that the angular momentum barrier is inefficient to prevent the direct formation of a supermassive black hole during the collapse. It is only in the limit of fully convective SMSs with f ≳ 1% that a significant fraction of the stellar mass (∼10%) can remain in orbit outside the horizon. Thus, the conditions of galaxy mergers appear to be more promising than those of atomically cooled haloes regarding the possibility for detectable gravitational wave emission and ultra-long gamma-ray bursts during the collapse.

Acknowledgments

L.H. has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 833925, project STAREX).

References

Appendix A: Newtonian pulsation analysis with slow rotation

We consider a Newtonian star in slow rotation, so that spherical symmetry can be assumed. Each layer r of the star rotates with a given angular velocity Ω(r), that is to say an average specific angular momentum of

$$\begin{array}{c}\hfill j=⟨r^2{sin}^2\theta \mathrm{\Omega }⟩=\frac{\int r^2{sin}^2\theta \mathrm{\Omega }sin\theta \mathrm{d}\theta \mathrm{d}\phi }{\int sin\theta \mathrm{d}\theta \mathrm{d}\phi }=\frac{2}{3}{r}^2\mathrm{\Omega }.\end{array}$$(A.1)

The average was determined by integration with respect to the solid angle sin θ dθdφ over the whole sphere r, and the additional sin²θ factor that appears in the numerator comes from the rotation radius r sin θ squared. In this case, the average of the radial, centrifugal acceleration is as follows:

$$\begin{array}{c}\hfill ⟨r{sin}^2\theta {\mathrm{\Omega }}^2⟩=\frac{\int r{sin}^2\theta {\mathrm{\Omega }}^{2}sin\theta \mathrm{d}\theta \mathrm{d}\phi }{\int sin\theta \mathrm{d}\theta \mathrm{d}\phi }=\frac{2}{3}r{\mathrm{\Omega }}^2=\frac{3}{2}\frac{j^2}{r^3}.\end{array}$$(A.2)

The two sin θ factors in the averaged quantity come from the rotational radius r sin θ and the projection of the centrifugal force on the radial direction.

We assume an equilibrium state with centrifugal acceleration (A.2) and denote the corresponding quantities with a subscript ₀:

$$\begin{array}{c}\hfill 0=-\frac{G{M}_r}{r_0^2}-4\pi r_0^2\frac{{\mathrm{d}P}_0}{\mathrm{d}{M}_r}+\frac{3}{2}\frac{j^2}{r_0^3}.\end{array}$$(A.3)

In Eq. (A.3), we adopted the Lagrangian mass as a coordinate, so that the subsript could be omitted for M_r. We considered small periodic perturbations to this equilibrium state, with local angular momentum conservation, which justifies the omission of the subscript for j. The dynamical quantities satisfy the Euler equation with centrifugal acceleration (A.2):

$$\begin{array}{c}\hfill \ddot{r}=-\frac{G{M}_r}{r^2}-4\pi r^2\frac{\mathrm{d}P}{\mathrm{d}{M}_r}+\frac{3}{2}\frac{j^2}{r^3}.\end{array}$$(A.4)

The amplitude of the periodic perturbations is expressed as an unknown function ϵ(M_r) ≪ 1, and the pulsation frequency ω is assumed to be constant:

$$\begin{array}{c}\hfill r=r_0(1+ϵe^{i\omega t}).\end{array}$$(A.5)