© The Authors 2024

1. Introduction

The observation of supermassive black holes with masses exceeding 10⁹ M_⊙ at redshifts of z > 6 has prompted an exploration of their formation process (e.g., Rees 1978; Mortlock et al. 2011; Woods et al. 2019). Possible building blocks of supermassive black holes are primordial (Pop III) supermassive stars (SMSs) of 10⁴ − 10⁶ M_⊙, which can be formed by very rapid accretion rates of Ṁ ≳ 0.1 M_⊙ yr⁻¹ (e.g. Hosokawa et al. 2013; Woods et al. 2017; Haemmerlé et al. 2018; Herrington et al. 2023; Nandal et al. 2024). Supermassive star formation with even higher accretion rates has been proposed by Mayer et al. (2010), Mayer & Bonoli (2019), and Zwick et al. (2023), who show that early galaxy collisions can torque gas at up to solar metallicities into supermassive nuclear disks (SNDs). These SNDs may experience radial gravitational instabilities, collapsing at rates of up to 10⁵ M_⊙ yr⁻¹. The central object formed under such extreme conditions is expected to be short-lived and to undergo “dark collapse” at masses above 10⁶ M_⊙ (Mayer & Bonoli 2019; Haemmerlé 2020; Haemmerlé et al. 2021).

The accreting object can become so massive that it encounters the general-relativistic (GR) instability (Iben 1963; Chandrasekhar 1964; Fowler 1966). For masses below 10⁵ M_⊙, the GR instability is unlikely to be encountered during the core hydrogen-burning phase (Woods et al. 2019). However, it may be reached in later evolutionary stages via the pair-instability process (Umeda et al. 2016; Woods et al. 2017, 2020). The collapse of such a star with a mass below 10⁵ M_⊙ may result in an explosion, neutrino emissions, and an ultra-long gamma ray burst, as per studies by Chen et al. (2014) and Nagele et al. (2021, 2022). For Pop III stars exceeding 10⁵ M_⊙, the GR instability is encountered either before or during the core hydrogen-burning phase. Once initiated by the GR instability, the collapse proceeds unimpeded, resulting in the star’s direct transition into a black hole without undergoing an explosion, thereby avoiding mass loss (Fuller et al. 1986; Montero et al. 2012).

Upper mass limits of SMSs due to the GR instability have been studied by different methods. Iben (1963) found the binding energy of static SMSs to decrease with mass and to become negative at a mass between 10⁵ and 10⁶ M_⊙. From the mass dependence of the total energy, Osaki (1966) found that SMSs more massive than 3.5 × 10⁵ M_⊙ encounter the GR instability before the onset of hydrogen burning. Umeda et al. (2016) determined the GR instability critical mass of an accreting SMS at a rapid hydrostatic contraction of the post-Newtonian (PN) hydrostatic structure. On the other hand, using evolution and hydrodynamic codes with PN corrections, Woods et al. (2017) and Herrington et al. (2023) detected the occurrence of GR instability as hydrodynamic collapses of their models. The critical masses obtained from hydrostatic calculations are systematically larger than those from hydrodynamic calculations (see Fig. 11 in Herrington et al. 2023). The latter seem appropriate, because the hydrostatic evolution code tends to average out small changes on short timescales that lead to large changes in the unstable structure. For this reason, GR hydrodynamics are needed to correctly detect the encounter with the GR instability during the SMS evolution.

Another way to find the occurrence of the GR instability is based on linear perturbation analysis, in which the small structure change is expressed as the consequence of a small radial displacement, ξ(r) e^iσct, from the equilibrium position. The displacement, ξ, is governed by a homogeneous differential equation with an eigenvalue, σ² (Chandrasekhar 1964) (Eq. (2) in Sect.3 below). We see the occurrence of the GR instability for a star if we solve the differential equation for the stellar structure and find σ² < 0. Solving the eigenvalue problem, Nagele et al. (2022) found that stars with 2 − 4 × 10⁴ M_⊙ (without accretion) become GR-unstable during or after the core helium-burning stage. On the other hand, based on a variational principle obeyed by the eigenvalue problem, Haemmerlé (2021) derived an integral expression (which consists of only structure variables) for the smallest σ² by using a trial function ξ/r≈ constant. Haemmerlé derived critical masses of 2.29 and 4.37 × 10⁵ M_⊙ for accretion rates of 1 and 10 M_⊙ yr⁻¹, respectively.

In this paper, we study the GR instability in rapidly accreting Pop III SMS models obtained by Nandal et al. (2024) under the PN gravity. For these models, we solve the adiabatic linear pulsation equation of Chandrasekhar (1964) to obtain ξ and σ².

In Section 2, we show briefly the properties of the Pop III accreting SMS models to which our stability analysis is applied. In Section 3, we discuss the equations employed in our stability analysis. We present the results in Section 4 and compare them with previous results in Section 5. In the appendix, we discuss briefly the numerical stability analysis applied to the n = 3 PN polytrope to test our PN pulsation code.

2. Primordial (Pop III) stellar evolution models with rapid mass accretion

Nandal et al. (2024) obtained Pop III SMS evolution models rapidly accreting with a range of 10⁻⁶ ≤ Ṁ_acc[M_⊙ yr⁻¹] ≤ 10³, using the Geneva stellar evolution code (GENEC : Eggenberger et al. 2008; Haemmerlé et al. 2016, 2018). It is assumed that the accretion of matter occurs via a geometrically thin cold disk, which implies that the specific entropy of the accreted matter is the same as that of the stellar surface. Furthermore, any excessive buildup of entropy during accretion is assumed to be radiated away by the surface of the star (Palla & Stahler 1992; Hosokawa et al. 2013). The PN gravity is included as per the work of Haemmerlé et al. (2018).

Computations of the models begin by accreting primordial matter on hydrostatic seeds. The composition of accreting matter is the same as that of the hydrostatic seed and consists of a homogeneous distribution of hydrogen with a mass fraction of X = 0.7516, of helium with Y = 0.2484, and deuterium at X(D) = 5 × 10⁻⁵ (Bernasconi & Maeder 1996; Behrend & Maeder 2001; Haemmerlé et al. 2016).

Among the models computed by Nandal et al. (2024), we have applied our stability analysis to models with accretion rates ranging from 0.01 M_⊙ yr⁻¹ to 1000 M_⊙ yr⁻¹. Evolution models with Ṁ_acc ≥ 0.01 M_⊙ yr⁻¹ are shown in Fig. 1 in the HR diagram (left panel) and on the mass-luminosity plane (right panel). Long-dashed lines indicate the relation of the zero-age main sequence (ZAMS) stars that have the same chemical composition. The ZAMS stars, with no accretion energy, are in thermal balance; in other words, the nuclear energy generation rate in the core is exactly balanced with the luminosity at the surface. In accreting stars, accretion heating causes the envelope to expand, resulting in a lower surface temperature compared to a ZAMS model with equivalent luminosity. The lower temperature is limited by the Hayashi limit at which a large range of the stellar envelope is in convective equilibrium. While the surface temperatures (or radii) are very different, the luminosity-to-mass relations of accreting stars are comparable to that of ZAMS.

Fig. 1. Hertzsprung–Russell (HR) diagram (left panel) and mass–luminosity relation (right panel) of evolutionary models for accretion rates between 0.01 and 1000 M⊙ yr−1. The large open circles indicate the positions at which the GR instability occurs, which is discussed in Section 4. The mass–luminosity relationships (right panel) converge toward the Eddington luminosity to mass relation, expressed as LEdd = 4πcGM/κel, where κel denotes the electron-scattering opacity.

The mass–luminosity relations shown in Fig. 1 (right panel) converge to a single relation in the high-luminosity range. The relation corresponds to the Eddington luminosity,

$$\begin{array}{c}\hfill L_{\mathrm{Edd}}=\frac{4\pi cGM}{{\kappa }_{\mathrm{el}}}=6.5\times {10}^4\frac{M/M_{\odot }}{(1+X)}{L}_{\odot },\end{array}$$(1)

where the electron-scattering opacity is given as κ_el = 0.2(1 + X) [cm² g⁻¹]. The occasional excesses of luminosity above L_Edd in a relatively low-luminosity range can be attributed to the effects of surface convective flux.

Although the luminosity-to-mass relation in Equation (1) has no limiting mass or luminosity, the GR instability occurs at a sufficiently large mass (e.g., Iben 1963; Chandrasekhar 1964; Osaki 1966). In this paper, we obtain the critical mass for each accretion rate by solving the equation of the GR adiabatic linear radial pulsation derived by Chandrasekhar (1964). We discuss the method of solving the differential equation in the next section.

3. General-relativistic equations for linear adiabatic radial pulsation

We examined the stability of an SMS model by solving the linear pulsation equation derived by Chandrasekhar (1964). Taking into account the GR effects, Chandrasekhar (1964) has derived the differential equation for infinitesimal radial displacement, ξe^iσct, as

$$\begin{array}{cc}\hfill \frac{\mathrm{d}}{\mathrm{d}r}\left[e^{3a+b}\frac{{\mathrm{\Gamma }}_{1}P}{r^2}\frac{\mathrm{d}}{\mathrm{d}r}(e^{-a}{r}^2\xi )\right]& =e^{2a+b}\xi [\frac{4}{r}\frac{\mathrm{d}P}{\mathrm{d}r}+\frac{8\pi G}{c^4}{e}^{2b}P(P+c^2\rho )\hfill \\ \hfill & -\frac{1}{P+c^2\rho }{\left(\frac{\mathrm{d}P}{\mathrm{d}r}\right)}^2-{\sigma }^2{e}^{2(b-a)}(P+c^2\rho )],\hfill \end{array}$$(2)

where P, ρ, and Γ₁ are, respectively, the pressure, matter density¹, and adiabatic exponent, (dlnP/dlnρ)_ad. Furthermore, a and b are metric coefficients, with which the line element, ds, is given as

$$\begin{array}{c}\hfill d{s}^2=-e^{2a}{c}^{2}d{t}^2+e^{2b}d{r}^2+r^2(d{\theta }^2+{sin}^2\theta d{\varphi }^2).\end{array}$$(3)

Since the static equilibrium values are sufficient for a and b in Eq. (2), they were obtained from the Einstein field equation for a spherically symmetric hydrostatic equilibrium as (Chandrasekhar 1964)

$$\begin{array}{c}\hfill e^{-2b}=1-\frac{2G{M}_r}{r{c}^2},\frac{e^{-2b}}{r}\frac{\mathrm{d}}{\mathrm{d}r}(a+b)=\frac{4\pi G}{c^4}(P+\rho c^2).\end{array}$$(4)

Because Eq. (3) at the stellar surface (r = R) should be equal to Schwarzshild’s metric, we have a relation of

$$\begin{array}{c}\hfill a(R)=\frac{1}{2}ln(1-\frac{2GM}{R{c}^2}).\end{array}$$(5)

Equation (2) is a linear eigenvalue equation with eigenvalue σ². The stellar structure is unstable if a pulsation mode has a negative eigenvalue; that is, σ² < 0.

To solve the second-order differential equation, we introduced, similarly to the Newtonian radial pulsations, two non-dimensional variables,

$$\begin{array}{c}\hfill Y_1\equiv \frac{\xi }{r}\mathrm{and}{Y}_2\equiv \frac{\mathrm{\Delta }P}{P},\end{array}$$(6)

where ΔP is the Lagrangian perturbation of pressure. Using these variables, Eq. (2) can be separated into the two first-order differential equations,

$$\begin{array}{c}\hfill \frac{\mathrm{d}{Y}_1}{\mathrm{d}lnr}=-(3-\frac{\mathrm{d}a}{\mathrm{d}lnr})Y_1-\frac{1}{{\mathrm{\Gamma }}_1}{Y}_2,\end{array}$$(7)

and

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{Y}_2}{\mathrm{d}lnr}& =[4\overline{V}-2\frac{\mathrm{d}(a+b)}{\mathrm{d}lnr}+\overline{V}\frac{\mathrm{d}a}{\mathrm{d}lnr}+{\omega }^{2}W{e}^{2(b-a)}]Y_1\hfill \\ \hfill & +[\overline{V}-\frac{\mathrm{d}(2a+b)}{\mathrm{d}lnr}]Y_2,\hfill \end{array}$$(8)

where

$$\begin{array}{c}\hfill \overline{V}\equiv -\frac{\mathrm{d}lnP}{\mathrm{d}lnr}\mathrm{and}{\omega }^{2}W\equiv {\sigma }^2{r}^2(1+\frac{c^2\rho }{P}),\end{array}$$(9)

with the square of non-dimensional pulsation frequency,

$$\begin{array}{c}\hfill {\omega }^2\equiv {\sigma }^2{c}^2{\left(\frac{\mathit{GM}}{R^3}\right)}^{-1}.\end{array}$$(10)

The boundary conditions of being finite at the center for Y₁ and at the surface for Y₂ were imposed. In addition, we adopted a normalization, Y₁(R) = 1, at the surface. The linear homogeneous differential equations (7) and (8) with an eigenvalue of ω² can be solved by the Henyey-type relaxation method often used in the analysis for Newtonian stellar pulsations (e.g., Unno et al. 1989).

We have obtained ω² for several stellar structures with accretion rates ranging from 0.01 to 1000 M_⊙ yr⁻¹. For a stellar model, we obtained a series of ω² depending on the number of nodes in ξ(r). The fundamental mode has no node in ξ and has the smallest ω². The sign of ω² of the fundamental mode determines the stability of the structure. If ω² < 0, we have pure imaginary frequencies, ±i|ω|, so that the GR instability occurs with a growth rate given by $|\omega |\sqrt{GM/R^3}$ in physical dimensions. If ω² of the fundamental mode is positive, the structure is stable, because perturbations cause only small-amplitude pulsations given as the real part of $[\xi (r)exp(i\omega \sqrt{GM/R^3}t)]$. In this case, $\omega \sqrt{GM/R^3}$ represents an angular frequency of adiabatic radial pulsation.

4. Results of the stability analysis

Figure 2 shows the stability of the fundamental modes in selected models transiting from stable to unstable structures (while models with Ṁ_acc = 0.01 M_⊙ yr⁻¹ remain stable). For an accretion rate of 0.05 M_⊙ yr⁻¹ or higher, the fundamental mode becomes unstable when the mass becomes sufficiently large, which is recognized by the steep increase in the growth rate from a negative value. The critical mass of the GR instability is larger in the models with a higher accretion rate, while the increase in the instability growth rate is steeper in the cases of lower accretion rates (and hence smaller masses).

Fig. 2. Growth rate of the GR instability versus stellar mass for each model. The negative part of the vertical axis corresponds to the angular frequency of the stable fundamental pulsation mode. The horizontal dashed line indicates the boundary between stable and unstable regions. A steep increase in the growth rate indicates that a strong GR instability occurs as soon as the mass exceeds the critical mass, depending on the accretion rate.

The stable side of Fig. 2 shows that the angular frequency, $\omega \sqrt{GM/R^3}$, is nearly constant in a wide range of stellar masses. We can understand this from the model properties seen in Fig. 1,

$$\begin{array}{c}\hfill L\approx L_{\mathrm{Edd}}\propto M\mathrm{and}L\propto R^2{T}_{\mathrm{eff}}^4\propto R^2.\end{array}$$

The second relation is the Stefan-Boltzmann law with an assumption of T_eff ∼ constant. Combining these relations and using ω ∼ 1 for the adiabatic fundamental mode, we obtained

$$\begin{array}{c}\hfill \omega \sqrt{GM/R^3}\propto M^{-\frac{1}{4}},\end{array}$$

which indicates that the adiabatic pulsation frequency has only a weak dependence on the mass. (The explanation would be improved if we took into account the slight luminosity dependence of T_eff).

Table 1 lists the obtained critical mass and luminosity as well as the central hydrogen abundance (X_c) at the occurrence of the GR instability for each accretion rate. (For an accretion rate of 0.01 M_⊙ yr⁻¹, the quantities of the final model are shown because we have never obtained the GR instability in this case). For the models with 0.05 M_⊙ yr⁻¹ ≲ Ṁ_acc ≲ 1000 M_⊙ yr⁻¹, the GR instability occurs during the core hydrogen-burning stage, while in the cases of Ṁ_acc > 1000 M_⊙ yr⁻¹, the GR instability would occur before the beginning of the core hydrogen burning. In contrast, hydrogen at the center of the models with Ṁ_acc ≲ 10⁻² M_⊙ yr⁻¹ is exhausted before the mass is increased enough for the GR instability. While we did not pursue further evolution, Nagele et al. (2022) found the GR instability to occur in (non-accreting) stars of 2 ∼ 3 × 10⁴ M_⊙ during or after the core helium-burning stage.

Figure 3 shows radial displacements, ξ (normalized as unity at the surface), of the fundamental mode as a function of the normalized distance from the center, r/R, for models of different masses accreting at 10 M_⊙ yr⁻¹. Models with M ≲ 4 × 10⁵ M_⊙ are stable and the displacement, ξ, of the fundamental mode is roughly proportional to r; that is, ξ/r≈ constant in the outer layers irrespective of the stellar mass. Models with M ≳ 4.5 × 10⁵ M_⊙ are unstable and the displacement, ξ, of the fundamental mode (ω² < 0) has a large peak in the interior where the GR effect is strongest. The radial distribution of ξ indicates that the GR instability leads to the collapse of the central part first. We note, however, that in the mass coordinate, the ξ peaks are located at M_r/M ≈ 0.92 and 0.98, respectively, for the models of M/M_⊙ = 6.7 × 10⁵ and 5.0 × 10⁵, indicating that the GR instability involves most of the mass in the SMS.

Fig. 3. Displacement, ξ, of the radial fundamental mode as a function of the distance from the center in some models with an accretion rate of 10 M⊙/yr. The displacement is normalised to ξ = 1 at the surface. For unstable modes (ω2 < 0), ξ is scaled by a factor to fit it in the figure.

Chandrasekhar (1964) obtained the condition for the occurrence of the GR instability in a polytrope of n = 3 (under the PN approximation) as

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_1-\frac{4}{3}<1.1245\frac{2GM}{c^{2}R}.\end{array}$$(11)

We have confirmed this relation numerically by applying our stability analysis (Sect. 3) to n = 3 PN polytropes for various values of M/R (see appendix). However, we cannot directly compare inequality (11) with the occurrence of the GR instability in SMS models, because Γ₁ varies in the interiors of SMSs, as is shown in the left panel of Fig. 4. Furthermore, because the outer layers of an accreting SMS model are extended compared with that of the n = 3 polytrope, 2GM/(c²R) at the surface of a GR unstable SMS (Table 1) is too small to compare with the inequality (11).

Fig. 4. Left: difference of the adiabatic exponent, Γ1, from 4/3 as a function of Mr/M for selected models. Solid and dashed lines indicate unstable and stable models against the GR instability, respectively. Accretion rates are color-coded, as is indicated. Right: quantity, Q, defined in Eq. (12) is plotted as a function of Mr/M for the structures of SMSs with various accretion rates. Unstable and stable structures are plotted with solid and short-dashed lines, respectively. These curves are compared with the long-dashed black line representing Qcrit (Eq. (13)) that corresponds to the Q for the critical n = 3 polytrope.

Taking into account the variable Γ₁, we considered the quantity, Q, given as

$$\begin{array}{c}\hfill Q\equiv \frac{R_{\mathrm{S}}}{r}\frac{1}{({\mathrm{\Gamma }}_1-4/3)}\mathrm{with}{R}_{\mathrm{S}}\equiv \frac{2G{M}_r}{c^2}.\end{array}$$(12)

The right panel of Fig. 4 shows Q as a function of the normalized mass coordinate, M_r/M, for selected GR unstable (solid lines) and stable models (short-dashed line), where different accretion rates are color-coded, as is indicated. As the evolution proceeds with an accretion rate, Q(M_r/M) increases as the total mass increases. When the mass becomes sufficiently large (depending on the accretion rate), the GR instability occurs so that the curve of Q(M_r/M) is switched to a solid line in Fig. 4. The stable-unstable transitions can be seen by very closely separated dashed and solid lines. It is interesting to note that Q of the unstable model (solid line) is always larger than that of the stable model (dashed line) throughout the stellar interior, even when the difference is very small, indicating that the maximum value of Q (which depends on the accretion rate) does indeed determine the mass at the first encounter with the GR stability.

For the critical (or neutrally stable) n = 3 PN polytrope, Q can be written as

$$\begin{array}{c}\hfill Q_{\mathrm{crit}}=\frac{1}{1.1245}\frac{M_r}{M}\frac{R}{r},\end{array}$$(13)

which is derived from Eq. (12) using Eq. (11) after replacing the inequality (< ) with the equality (because of the neutral condition). The polytropic Q_crit is shown by a long-dashed black line in the right panel of Fig. 4, for comparison with Q(M_r/M) of accreting SMS models. It is interesting to note that the maximum of Q_crit ( ≈ 1.7) of the n = 3 polytrope is the lower bound for the max(Q) of the GR unstable accreting SMSs. In other words, the polytropic relation provides the necessary condition, max(Q) > 1.7, for accreting SMS models to be GR unstable. The max(Q) of the unstable less massive SMS for Ṁ_acc = 0.05 − 0.1 M_⊙ yr⁻¹ is comparable with max(Q_crit), while the max(Q) of a more massive SMS with a higher accretion rate (Ṁ_acc ≥ 1 M_⊙ yr⁻¹) is considerably higher than the polytropic limit. This property is consistent with the fact (left panel of Fig. 4) that Γ₁ − 4/3 is nearly constant throughout the interior in the less massive SMS (i.e., structure is comparable with a polytrope), while it varies considerably in more massive SMSs. We note, in passing, that Fig. 4 clearly indicates that the 2 × 10⁴ M_⊙ model for the accretion rate of 0.01 M_⊙ yr⁻¹ does not satisfy the necessary condition of the GR instability.

5. Comparison with previous results

Figure 5 shows our GR instability-critical masses (open squares) versus accretion rates in comparison with the results of previous works. Umeda et al. (2016) determined the GR instability critical mass as the mass when the central temperature has increased to 10^9.2 K due to hydrostatic contraction. Their masses are systematically larger than those of our models, indicating that GR instability collapse should occur earlier if the hydrodynamic effects included. This is consistent with the fact that the critical masses from hydrodynamic models obtained by Herrington et al. (2023) (triangles) are smaller than those obtained from their hydrostatic models (not shown in Fig. 5), which are comparable with those of Umeda et al. (2016). On the other hand, the critical masses based on the hydrodynamic calculations of Herrington et al. (2023) (triangles in Fig. 5) are comparable with our critical masses.

Fig. 5. Critical masses of GR instability for models with various accretion rates (horizontal axis) are compared with results in the literature. For Herrington et al. (2023), results from hydrodynamics models (read from their Fig.11) are shown.

Filled green circles in Fig. 5 show the GR instability critical masses obtained by Woods et al. (2017) using another hydrodynamic evolution code with PN corrections. Despite the use of a hydrodynamic code, the green dots (for Ṁ ώ 1 M_⊙ yr⁻¹ in particular) deviate from the relation of Herrington et al. (2023) (filled blue triangles) as well as from our relation (open squares). The cause of the deviation is not clear.

The results of Haemmerlé (2021) for the models with accretion rates of 1 M_⊙ yr⁻¹ and 10 M_⊙ yr⁻¹ (diamond symbols in Fig. 5) agree well with ours. This is reasonable because his stability analysis is based on a variational property of Eq. (2) in Sect. 3 with a trial function ξ/r≈constant, which is a reasonable approximation before the occurrence of the GR instability, as is seen in Fig. 3. In addition, his evolution models were obtained by GENEC code, which was also used for the models in Nandal et al. (2024).

6. Conclusions

We have examined when the GR instability occurs in a rapidly accreting Pop III SMS by solving the linear adiabatic radial pulsation equation derived by Chandrasekhar (1964). We find that the displacement, ξ, of the fundamental mode of pulsation monotonically increases from the center to the surface in a stable model, while it has a large peak between the center and the surface in an unstable model.

For each accretion rate (≳0.05 M_⊙ yr⁻¹), there is a critical mass above which the GR instability occurs and the star begins to collapse. The critical mass of an accreting SMS is larger for a higher accretion rate. We find that 8 × 10⁴ M_⊙ (at the end of the core hydrogen burning) for Ṁ> = 0.05 M_⊙ yr⁻¹ and 10⁶ M_⊙ (at the beginning of the core hydrogen burning) for 1000 M_⊙ yr⁻¹.

Comparisons with previous results in the literature show that our results agree well with Haemmerlé (2021)’s results for 1 and 10 M_⊙ yr⁻¹, and hydrodynamic calculations of Herrington et al. (2023) for 0.1 to 1 M_⊙ yr⁻¹. Our present work has extended critical masses for the GR instability from those of the previous works to 7 × 10⁵ and 10⁶ M_⊙ for the accretion rates of 100 and 1000 M_⊙ yr⁻¹, respectively.

Acknowledgments

We are grateful to the referee, Dr. Chris Nagele for very useful comments and suggestions, which have greatly improved the manuscript. We also thank Lionel Haemmerlé for useful comments. The authors of this paper have received supports from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement No 833925, project STAREX).

References

Appendix A: Stability analysis for Post-Newtonian Polytrope

To test our radial pulsation code, we have applied it to an n = 3 polytrope under the Post-Newtonian approximation for which the critical condition for the GR instability has been derived by Chandrasekhar (1964). In this appendix we summarize equations needed in the adiabatic radial pulsation analysis for a polytrope under the Post-Newtonian approximation. We refer the reader Chandrasekhar (1964) for details.

A.1. Equilibrium Post-Newtonian Polytrope

The pressure P and the density ρ are expressed as

$$\begin{array}{c}\hfill P=P_{\mathrm{c}}{\mathrm{\Theta }}^{n+1}\mathrm{and}\rho ={\rho }_{\mathrm{c}}{\mathrm{\Theta }}^n,\end{array}$$(A.1)

where the subscript ’c’ indicates the value at the center, and n is the polytropic index. In the Post-Newtonian (PN) approximation, Θ is expressed as

$$\begin{array}{c}\hfill \mathrm{\Theta }=\theta +q\varphi \mathrm{with}q=\frac{P_{\mathrm{c}}}{{\rho }_{\mathrm{c}}{c}^2}=\frac{{\eta }_1}{2(n+1)v_1}\frac{2GM}{R{c}^2},\end{array}$$(A.2)

where θ is the Lane-Emden function of η; dimension-less distance from the center defined as

$$\begin{array}{c}\hfill \eta =r/\alpha \mathrm{with}\alpha \equiv {\left[\frac{(n+1)q{c}^2}{4\pi G{\rho }_{\mathrm{c}}}\right]}^{1/2}.\end{array}$$(A.3)

The subscript ‘₁’ indicates the value at first zero of θ. The classical Lane-Emden function θ and an additional function ϕ representing the PN effects satisfy the following equations:

$$\begin{array}{c}\hfill \frac{1}{{\eta }^2}\frac{\mathit{dv}}{d\eta }={\theta }^n\mathrm{with}v=-{\eta }^2\frac{d\theta }{d\eta }\end{array}$$(A.4)

and

$$\begin{array}{c}\hfill {\eta }^2\frac{d\varphi }{d\eta }=-[\theta +2(n+1)\frac{v}{\eta }]v-{\eta }^3{\theta }^{n+1}-w\mathrm{with}\frac{\mathit{dw}}{d\eta }=n{\eta }^2{\theta }^{n-1}\varphi .\end{array}$$(A.5)

We have obtained functions θ(η), v(η), ϕ(η) and w(η) by integrating above differential equations from a point very near the center to the first zero poing of θ. For the initial values at the starting point where η ≪ 1, we have used the following relations:

$$\begin{array}{c}\hfill \theta \approx 1-\frac{\eta }{6},\varphi \approx -\frac{2}{3}{\eta }^2,\mathrm{and}w\approx -\frac{2n}{15}{\eta }^5.\end{array}$$(A.6)

A.2. Adiabatic radial pulsation analysis

Having polytropic functions, θ, v, ϕ, and w, we can evaluate metrics a and b as

$$\begin{array}{c}\hfill e^{-2b}=1-2(n+1)q\frac{v}{\eta }\mathrm{and}{e}^{2a}=1-2(n+1)q(\theta +\frac{v_1}{{\eta }_1}).\end{array}$$(A.7)

Using these relations we express quantities needed in the differential equations (Eqs. 7 and 8) for Y₁ and Y₂ as

$$\begin{array}{c}\hfill \frac{\mathit{da}}{dln\eta }=(n+1)q\frac{v}{\eta }\mathrm{and}\frac{d(a+b)}{dln\eta }=(n+1)q{\eta }^2{\theta }^n,\end{array}$$(A.8)

$$\begin{array}{c}\hfill \overline{V}=\frac{(n+1)}{\theta +q\varphi }[\frac{v}{\eta }+q(\theta +2(n+1)\frac{v}{\eta })\frac{v}{\eta }+q{\eta }^2{\theta }^{n+1}+q\frac{w}{\eta }]\end{array}$$(A.9)

and

$$\begin{array}{c}\hfill W{e}^{2(b-a)}=\frac{(n+1)v_1{\eta }^2(q+{\theta }^{-1})}{{\eta }_1^3[1-2(n+1)q(\theta +v/\eta +v_1/{\eta }_1)]}.\end{array}$$(A.10)

For the numerical test we have adopted n = 3. For each assumed values of M/R’s we determined the critical value of $({\mathrm{\Gamma }}_1-\frac{4}{3})$ for the neutral stability, and compared with the theoretical limit 1.1245 × 2GM/(c²R). As shown in Fig. A.1, we have obtained a very good agreement, which give us confidence in the discussions presented in the main text of the present paper.

Fig. A.1. Numerically determined critical ${\mathrm{\Gamma }}_1-\frac{4}{3}$ for the GR instability for PN n = 3 polytropes versus analytical prediction $1.1245\frac{2GM}{c^{2}R}$.

All Tables

All Figures

Fig. 1. Hertzsprung–Russell (HR) diagram (left panel) and mass–luminosity relation (right panel) of evolutionary models for accretion rates between 0.01 and 1000 M⊙ yr−1. The large open circles indicate the positions at which the GR instability occurs, which is discussed in Section 4. The mass–luminosity relationships (right panel) converge toward the Eddington luminosity to mass relation, expressed as LEdd = 4πcGM/κel, where κel denotes the electron-scattering opacity.

In the text

Fig. 2. Growth rate of the GR instability versus stellar mass for each model. The negative part of the vertical axis corresponds to the angular frequency of the stable fundamental pulsation mode. The horizontal dashed line indicates the boundary between stable and unstable regions. A steep increase in the growth rate indicates that a strong GR instability occurs as soon as the mass exceeds the critical mass, depending on the accretion rate.

In the text

	Fig. 3. Displacement, ξ, of the radial fundamental mode as a function of the distance from the center in some models with an accretion rate of 10 M⊙/yr. The displacement is normalised to ξ = 1 at the surface. For unstable modes (ω2 < 0), ξ is scaled by a factor to fit it in the figure.
In the text

Fig. 4. Left: difference of the adiabatic exponent, Γ1, from 4/3 as a function of Mr/M for selected models. Solid and dashed lines indicate unstable and stable models against the GR instability, respectively. Accretion rates are color-coded, as is indicated. Right: quantity, Q, defined in Eq. (12) is plotted as a function of Mr/M for the structures of SMSs with various accretion rates. Unstable and stable structures are plotted with solid and short-dashed lines, respectively. These curves are compared with the long-dashed black line representing Qcrit (Eq. (13)) that corresponds to the Q for the critical n = 3 polytrope.

In the text

	Fig. 5. Critical masses of GR instability for models with various accretion rates (horizontal axis) are compared with results in the literature. For Herrington et al. (2023), results from hydrodynamics models (read from their Fig.11) are shown.
In the text

	Fig. A.1. Numerically determined critical ${\mathrm{\Gamma }}_1-\frac{4}{3}$ for the GR instability for PN n = 3 polytropes versus analytical prediction $1.1245\frac{2GM}{c^{2}R}$.
In the text