FORMATION OF PRIMORDIAL SUPERMASSIVE STARS BY RAPID MASS ACCRETION

We study the evolution of rapidly accreting stars by numerically solving their interior structure. The governing equations are the usual stellar structure equations, including effects of the mass accretion (e.g., Kippenhahn & Meyer-Hofmeister 1977; Stahler et al. 1980). In this paper, we employ the numerical code developed by Yorke & Bodenheimer (2008), where the Henyey method (Henyey et al. 1964) is employed to solve the stellar interior structure (also see Bodenheimer et al. 2007). The outer boundary conditions for the interior are provided by solving the structure of a gray atmosphere. To do this, we integrate inward the equations of hydrostatic equilibrium and radiative or convective transfer over the atmosphere with a Runge–Kutta procedure. The stellar effective temperature T_eff is defined at the optical depth τ = 2/3 measured from the exterior. The atmosphere integration extends to a connecting point at a temperature T_atm ≃ several × T_eff, which varies slowly with time. We presume that the thermal state of the accreting material is the same as in the outer atmosphere (see below). The material added to the atmosphere is ultimately incorporated into the interior by means of an automatic rezoning procedure.

The current numerical code differs from that used in HOY12 (also see Hosokawa & Omukai 2009), where the same governing equations were solved by a shooting method instead of the Henyey method used here. Inside, an SMS radiation dominates the total pressure, but a small contribution by the gas is crucial in determining its structure, such as the radius. Here we adopt the gas pressure as one of our independent Henyey variables in place of the total pressure. This enables us to accurately solve the structure of an SMS. The outer boundary conditions are also different. Most of the cases examined in HOY12 assume spherical accretion onto a star (e.g., Stahler et al. 1980; Hosokawa & Omukai 2009), whereby steady accretion flow hits the stellar surface forming an accretion shock front. In this case, we also consider the structure of the steady accretion flow falling onto the star. We calculate the radial structure of both the outer accretion flow and stellar interior for satisfying the shock jump conditions.

HOY12 also study a few cases with a different outer boundary, the photospheric boundary conditions:

$$\begin{equation} P_{\rm surf} = \frac{2}{3} \frac{1}{\kappa _{\rm surf}} \frac{{G M}_*}{R_*^2}, \end{equation} \tag{ 2 }$$

and

$$\begin{equation} L_* = 4 \pi R_*^2 \sigma T_{\rm surf}^4 \end{equation} \tag{ 3 }$$

(e.g., Palla & Stahler 1992; Hosokawa et al. 2010), where the suffix "surf" represents quantities at the stellar surface, and σ is the Stefan–Boltzmann constant. This boundary is rather similar to that adopted in the current code, except that the structure of the atmosphere is not solved. These outer boundary conditions assume that the accreting material is added to the star through a geometrically thin disk, whereby most of the stellar surface can still freely radiate (no backwarming).

Although one can see that the photospheric boundary is realistic for the assumed disk accretion, it is important that the outer boundary condition correctly determines the specific entropy of accreting material (e.g., Hartmann et al. 1997; Hosokawa et al. 2011a). With the photospheric boundary, the accreting material has the smallest amount of entropy (so-called cold accretion), which is the value at the stellar surface. In this case, accreting gas would slowly approach the star through the disk, radiating its heat away before reaching the stellar surface. The accreting gas gently lands on the stellar surface with the same entropy as in the outer photosphere. With the spherically symmetric shock boundary condition, on the other hand, the accreting material has a much higher entropy than for the cold accretion case, carrying a larger fraction of the entropy generated in the shock front into the stellar interior. However, even for the case of highly non-spherical accretion through a circumstellar disk, the very rapid mass accretion considered here means that a non-negligible amount of entropy should be advected into the star. For this reason, HOY12 adopt the shock boundary condition for the fiducial cases.

Detailed modeling of the thermal properties of the accretion flow is beyond the scope of spherically symmetric stellar evolution calculations. Realizing that a fraction of the gravitational energy released by the accreting gas will be deposited into the stellar interior (e.g., Siess & Forestini 1996; Siess et al. 1997), we add an additional luminosity

$$\begin{equation} L_{*,{\rm acc}} \equiv \eta L_{\rm acc} = \eta \frac{{G M}_* \dot{M}_*}{R_*} \end{equation} \tag{ 4 }$$

to the bottom of the atmosphere. The fraction of the deposited accretion luminosity η is treated as a free parameter with a value 0 ⩽ η ⩽ 1. The limit of cold accretion corresponds to η = 0. More detailed modeling would be possible if we were to consider how the deposited energy is spatially distributed in the stellar interior (e.g., Siess & Forestini 1996), i.e., not only at the bottom of the atmosphere. However, we find that the evolution does not depend on η for M_* ≳ 100 M_☉.

Table 1 summarizes the cases considered in this paper. We follow the stellar evolution with rapid mass accretion at various constant rates of 0.01, 0.1, 0.3, and 1.0 M_☉ yr⁻¹ beginning with a 2 M_☉ initial model. Cases with different values of η are also examined for each accretion rate. Although the current numerical code enables following the evolution far beyond that considered by HOY12, we ultimately encounter numerical convergence difficulties. Here we present the evolution for ≃ 10⁵ yr in each case. The stellar mass reaches 10⁵ M_☉ with the accretion rate $\dot{M}_*= 1.0 \,M_\odot {\rm \,yr}^{-1}$ (cases a1, see Table 1).

Table 1. Cases Considered

Cases	$\dot{M}_*$	η	M*, f
	(M☉ yr−1)		(M☉)
a1-e0	1.0	0.0	105
a1-e0.01	1.0	0.01	105
a1-e0.1	1.0	0.1	105
a1-e1	1.0	1.0	105
a0.3-e0.1	0.3	0.1	6 × 104
a0.1-e0	0.1	0.0	2 × 104
a0.1-e0.01	0.1	0.01	2 × 104
a0.1-e0.1	0.1	0.1	2 × 104
a0.1-e1	0.1	1.0	2 × 104
a0.01-e0.1	0.01	0.1	2 × 103
a0.01-e0	0.01	0.0	2 × 103

Notes. Column 2: mass accretion rate, Column 3: fraction of the accretion luminosity deposited in the stellar interior (see Equation (4)), Column 4: final stellar mass.

Download table as:  ASCIITypeset image

The cases with the lowest accretion rate of 0.01 M_☉ yr⁻¹ (cases a0.01) are relevant to the normal mode of the primordial star formation, where H₂ molecular cooling operates. We also study these cases as our previous calculations predict that, with such moderately high accretion rates, the stellar radius shows a peculiar behavior, oscillating with increasing stellar mass (e.g., Omukai & Palla 2001, 2003; HOY12). Below, we show that our current numerical code also provides a similar evolution, despite a number of differences, e.g., the outer boundary conditions.

3.1.1. Outline of the Evolution

Figure 1 presents the overall stellar evolution with the accretion rates $\dot{M}_*= 0.1 \,M_\odot {\rm \,yr}^{-1}$ (case a0.1-e0.1) and 1.0 M_☉ yr⁻¹ (case a1-e0.1), where the deposited energy fraction η is set at 0.1. We see that, in both cases, the early evolution of the stellar radius agrees well with our previous results for M_* ≲ 10³ M_☉ (HOY12). The small deviations for M_* ≲ 100 M_☉ are mostly due to the different outer boundary conditions as explained in Section 2.1 (also see Section 3.1.4 below). Our current results show that the star continues to expand for M_* ≳ 10³ M_☉. The stellar radius increases monotonically with mass according to relation (1) for 100 M_☉ ≲ M_* ≲ 10⁴ M_☉, and finally begins to decrease for M_* ≳ 3 × 10⁴ M_☉. Nonetheless, the stellar radius remains very large at R_* ≃ 2 × 10⁴ R_☉ ≃ 100 AU. With the accretion rate of $\dot{M}_*= 1.0 \,M_\odot {\rm \,yr}^{-1}$, we obtain an accreting 10⁵ M_☉ star which is still significantly bloated.

Zoom In Zoom Out Reset image size

To understand why the radius begins to decrease for M_* ≳ 10⁴ M_☉, it is helpful to recall how Equation (1) is derived (also see HOY12). In general, the luminosity of very massive stars is close to the Eddington value,

$$\begin{equation} L_* = 4 \pi R_*^2 \sigma T_{\rm eff}^4 \simeq L_{\rm Edd} \propto M_*. \end{equation} \tag{ 5 }$$

Because of the very strong temperature-dependence of gas opacity with H⁻ bound-free absorption, which dominates near the stellar surface, the effective temperature is locked at several ×10³ K (e.g., Hayashi 1961). Equation (1) is derived by substituting T_eff = 5000 K into Equation (5). The upper panel of Figure 1 also shows that the effective temperature is certainly locked at ≃ 5000 K for M_* ≲ 10⁴ M_☉. However, we see that the effective temperature gradually increases for M_* ≳ 10⁴ M_☉, when the stellar expansion ceases. The deviation of the stellar expansion from relation (1) is due to the rise of effective temperature, which is considered separately in Section 3.1.3 below.

In general, the evolution of an accreting star can be understood by comparing the accretion time

$$\begin{equation} t_{\rm acc} \equiv \frac{M_*}{\dot{M}_*}, \end{equation} \tag{ 6 }$$

the timescale for stellar mass growth, to the Kelvin–Helmholtz (KH) time

$$\begin{equation} t_{\rm KH} \equiv \frac{{G M}_*^2}{R_* L_*}, \end{equation} \tag{ 7 }$$

over which the star radiates away its gravitational energy (e.g., Stahler et al. 1986). The lower panel of Figure 1 shows that the balance between these two timescales is inverted from t_KH > t_acc to t_KH < t_acc during the evolution, e.g., at M_* ≃ 200 M_☉ for 1.0 M_☉ yr⁻¹. As explained in the literature (e.g., Omukai & Palla 2003; Hosokawa & Omukai 2009), this is because the KH timescale decreases as the stellar luminosity L_* rapidly increases with stellar mass. The same timescale inversion is also shown in HOY12 (see their Figure 3).

We see that the accreting stars have a very short KH timescale compared to the accretion timescale for M_* ≳ a few × 10² M_☉. This means that the star efficiently loses its energy via radiation, which usually leads to contraction (so-called KH contraction). However, our results show that the stellar radius continues to increase until M_* ≳ 3 × 10⁴ M_☉. This suggests that the star has a very inhomogeneous structure, e.g., most of the stellar interior contracts while radiating energy away, whereas a surface layer expands as it receives part of the outward heat flux. HOY12 show this core-envelope structure of rapidly accreting stars for M_* ⩽ 10³ M_☉. In Section 3.1.2 below, we show that the same stellar structure is maintained even after the star becomes supermassive, M_* ≳ 10⁴ M_☉.

3.1.2. Stellar Interior Structure

Here we investigate the evolution of the interior structure of accreting SMSs. Figure 2 shows the evolution with the accretion rate $\dot{M}_*= 1.0 \,M_\odot {\rm \,yr}^{-1}$ and input energy fraction η = 0.1 (case a1-e0.1). As explained in Section 3.1.1, the KH timescale falls below the accretion time at M_* ≃ 200 M_☉ in this case.

Zoom In Zoom Out Reset image size

Although deuterium burning begins soon after the timescale inversion, this hardly influences the subsequent evolution (also see HOY12). At this stage, deuterium burns at nearly a steady rate, because it is replenished by accretion as quickly as it is consumed by burning,

$$\begin{equation} L_{\rm D,st} \equiv \dot{M}_*\delta _{\rm D} = 1.5 \times 10^6 \,L_{\odot }\left(\frac{\dot{M}_*}{1 \,M_\odot {\rm \,yr}^{-1}} \right) \left(\frac{[{\rm D}/{\rm H}]}{2.5 \times 10^{-5}} \right), \end{equation} \tag{ 8 }$$

where δ_D is the energy produced via deuterium burning per unit mass. The energy released by deuterium burning is negligible in comparison to the local luminosity, approximately given by the Eddington value. The heat generated by nuclear fusion is rapidly transported via radiation without causing convection. This is not the case with extremely cold mass accretion, η = 0 (see below).

Figure 2 shows that for M_* ≳ 200 M_☉ the star has an inhomogeneous mass distribution; a surface layer, which contains only a tiny fraction of the total stellar mass, significantly inflates. Moreover, we see that the star becomes increasingly centrally condensed as the stellar mass increases. This is because most of the stellar interior is contracting in this stage, as explained in Section 3.1.1. The gravitational energy released by the contracting interior is radiatively transported outward. However, part of the outward flux is absorbed in a surface layer where the opacity is higher than in the interior. As a result, the surface layer has a high specific entropy and largely inflates. Figure 2 shows that the energy is partly transported by convection in this semi-opaque surface layer. Figure 3(a) presents the characteristic profiles of the specific entropy in the stellar interior, showing that the entropy increases outward to reach a maximum at the bottom of the surface convective layer.

Zoom In Zoom Out Reset image size

Figure 4(a) shows that, as in the ordinary KH contraction stage, the stellar central temperature increases with mass. Hydrogen burning begins at M_* ≃ 4 × 10³ M_☉, after which the central temperature is fixed at ≃ 1.5 × 10⁸ K by the thermostat effect, due to the strong temperature-dependence of energy production through the CNO cycle. Although it is initially at zero metallicity, a small amount of carbon is produced in the center through the 3-α process during KH contraction, which allows efficient hydrogen nuclear burning via the CNO cycle. Figures 2 and 4 show that a convective hydrogen-burning core grows with increasing stellar mass. The ignition of hydrogen burning halts the KH contraction within and near the convective core. In fact, Figure 4(b) shows that the ratio of the central density to the average density decreases with increasing mass after hydrogen burning begins. However, most of the stellar interior is still contracting, as can be inferred from Figure 2. As the mass increases from 10⁴ M_☉ to 10⁵ M_☉ the radii of 80% and 60% of the total mass continue to decrease, and the radii of 40% and 20% increase only slightly, though much less than the growth in mass. Finally, at M_* ≃ 10⁵ M_☉ the outer 99% of the star in radius contains only 25% of its total mass (see Figure 4(c)). Figures 3(b) and (c) show that the average interior density and temperature still rise with increasing stellar mass as a result of the ongoing contraction for M_* ≳ 10⁴ M_☉. The gravitational energy released by the contracting envelope compensates most of the stellar surface luminosity, even after the ignition of hydrogen burning. Most of the heat generated by nuclear fusion is absorbed in the convective core, whose entropy is consequently elevated as shown in Figure 4(a).

Zoom In Zoom Out Reset image size

Finally, we consider how the above evolution changes with a lower accretion rate $\dot{M}_*= 0.1 \,M_\odot {\rm \,yr}^{-1}$. As can be inferred from Figure 1, the evolution is qualitatively the same as in the previously discussed cases with $\dot{M}_*= 1.0 \,M_\odot {\rm \,yr}^{-1}$. There is a quantitative difference, however. At the lower accretion rate the accretion timescale is correspondingly much longer than the KH timescale for M_* ≳ 100 M_☉. The star has a more centrally concentrated mass distribution for a given stellar mass, because it has had more time to contract and has radiated more energy ∫L_* dt, which has to be covered by gravitational contraction. Figure 4(a) shows that, as a result, the central temperature is higher at comparable stellar masses, so that hydrogen burning begins before the stellar mass reaches 10³ M_☉. The outer envelope of the star still continues to contract after the ignition of hydrogen, as previously described. Figures 3(b) and (c) show that, for a 10⁴ M_☉ star, the average interior temperature and density for $\dot{M}_*= 0.1 \,M_\odot {\rm \,yr}^{-1}$ are much higher than those with $\dot{M}_*= 1.0 \,M_\odot {\rm \,yr}^{-1}$. It is interesting that, regardless of the different interior structure, the evolution of the stellar radius, while obeying relation (1), is almost independent of the accretion rate.

3.1.3. Stellar Surface Structure and Termination of the Expansion

Our results show that rapidly accreting SMSs have a large radius ∼100 AU, even after the stellar mass exceeds 10⁴ M_☉. However, the stellar expansion gradually slows down, no longer following the analytic mass–radius relation (1). Figure 1 shows that, especially with $\dot{M}_*= 1.0 \,M_\odot {\rm \,yr}^{-1}$, the stellar radius finally turns around and decreases slightly for M_* ≳ 3 × 10⁴ M_☉. Here we consider why the stellar expansion finally ceases.

As explained in Section 3.1.2, the stellar structure is not homogeneous during expansion. Only a surface layer which has a tiny fraction of the total stellar mass inflates. Whereas most of the stellar interior contracts and radiates the energy away, the outer layer absorbs a part of the outward flux and gains a relatively high specific entropy, resulting in the surface-layer bloating. The physical state of the outermost part of a star is thus a key to understanding the termination of the expansion. Figure 5 presents several snapshots of the outermost radial profiles of pressure, density, and opacity as functions of the interior temperature T for case a1-e0.1 ($\dot{M}_*= 1.0 \,M_\odot {\rm \,yr}^{-1}$ and η = 0.1). Note that higher interior temperatures represent deeper stellar layers, because the interior temperature increases monotonically toward the center. For M_* ≲ 10⁴ M_☉ while the stellar radius is still increasing with mass according to Equation (1), the star has a characteristic feature for T ≲ 3 × 10⁴ K. Panel (a) shows that the gas pressure dominates over the radiation pressure there, shaping the almost flat profile within the photosphere. We also see that the density accordingly increases outward (panel (b)), i.e., it is not monotonic. This structure is called a density inversion. Such density inversions occasionally occur in the outermost layers where the radiative flux locally exceeds the Eddington value. The hydrostatic balance is achieved with the inward gas pressure gradient (around T ≃ 2 × 10⁴ K in panel (a)), which helps gravity balance the outward radiation force. The gas opacity has a peak at T ≃ 10⁴ K because of the bound-free absorption of H and H⁻ (e.g., Kippenhahn et al. 2013). In the outer layers where T ≲ 10⁴ K, the opacity sharply drops with decreasing temperature, a general feature of H⁻ absorption, because at lower temperatures there are fewer free electrons available for forming H⁻ ions.

Zoom In Zoom Out Reset image size

However, this characteristic feature at the surface changes when the stellar mass exceeds M_* > 10⁴ M_☉. Radiation pressure starts to dominate the total pressure everywhere in the interior, and the density inversion almost disappears by the epoch of M_* ≃ 10⁵ M_☉. This behavior is due to the density-dependence of gas opacity. Figure 5(c) shows that, at the low density ρ ≲ 10⁻¹¹ g cm⁻³, opacity in the surface layers drops to the Thomson scattering value as H⁻ absorption becomes inefficient (also see, e.g., Mayer & Duschl 2005). This is inevitable because collisional events between neutral hydrogen atoms and free electrons are required to form H⁻ ions and the H⁻ opacity is therefore proportional to n_Hn_e.

Recall that Equation (1) is derived assuming that the effective temperature is locked at T_eff ≃ 5000 K because of the strong temperature-dependence of H⁻ opacity. At the point when H⁻ ions no longer dominate the opacity in surface layers, these surface layers no longer efficiently absorb the heat flux coming from the contracting interior. We thus predict that the stellar radius will continue to decrease after the stellar mass exceeds 10⁵ M_☉, in contrast to analytical considerations by Schleicher et al. (2013), who do not consider the detailed structure in the outermost part of the star, where H⁻ opacity is important. We stress that the density-dependence of H⁻ opacity is a key to understanding the termination of the stellar expansion.

3.1.4. Evolution with Different Input Energy Fractions η

As described in Section 2.1, we parameterize the gravitational energy deposited into the stellar interior with the free parameter η. Here, we investigate how the stellar evolution varies with different values of η. We naively expect that, at a fixed stellar mass, the stellar radius would be larger with higher η, as the specific entropy of accreting materials would be enhanced. Figure 6 shows this is only true for early evolutionary stages when M_* ≲ 200 M_☉ (30 M_☉) for the accretion rate $\dot{M}_*= 1.0 \,M_\odot {\rm \,yr}^{-1}$ (0.1 M_☉ yr⁻¹, respectively), which corresponds to the epoch when the accretion timescale is shorter than the KH timescale (see Section 3.1.1). With this timescale imbalance, gas in the stellar interior retains the entropy originally obtained when accreted (Stahler et al. 1986), which in our case is determined by η. Differences of η thus influence the stellar structure during this stage. After the timescale inversion to t_acc > t_KH, the accretion luminosity L_acc becomes much smaller than the stellar luminosity L_*. Note that the timescale ratio t_acc/t_KH is equal to the luminosity ratio L_*/L_acc by definition. Thus, the additional heat input into the atmosphere ηL_acc is negligible.

Zoom In Zoom Out Reset image size

We find that the previous results of HOY12 with the shock boundary condition are similar to the current cases with η = 0.1. For the shock outer boundary, the specific entropy of accreting material is determined by solving the flow structure across the accretion shock front satisfying jump conditions (e.g., Stahler et al. 1980). We have no free parameters in this case. For the cold-accretion limit with η = 0, the stellar radius initially decreases with increasing stellar mass, until the star abruptly expands, e.g., at M_* ≃ 50 M_☉ for $\dot{M}_*= 1 \,M_\odot {\rm \,yr}^{-1}$. HOY12 also find similar evolutionary behavior with the photospheric boundary conditions (2) and (3). This abrupt expansion in the cold-accretion limit is related to the extension of a deuterium-burning layer in the stellar interior, which is studied in detail in Hosokawa et al. (2010). We shall also examine the interior structure in the context of the current results in Section 3.1.2.

The lower panel of Figure 6 depicts the evolution of the stellar radius at the lowest accretion rate considered here $\dot{M}_*= 10^{-2} \,M_\odot {\rm \,yr}^{-1}$ with the input energy fractions η = 0.1 and 0.0 (cases a0.01-e0.1 and a0.01-e0). In these cases, the timescale inversion to t_KH < t_acc occurs around the first local maximum of the radius at M_* ≃ 15 M_☉. The star contracts for a while after this point, radiating the energy away (KH contraction). As briefly mentioned in Section 2.2, however, we see a peculiar evolution for M_* ≳ 100 M_☉, whereby the stellar radius slowly oscillates as the stellar mass increases. Omukai & Palla (2003) show that a similar evolution occurs with $\dot{M}_*\gtrsim 4 \times 10^{-3} \,M_\odot {\rm \,yr}^{-1}$ as the total stellar luminosity L_tot ≡ L_* + L_acc approaches the Eddington value during the KH contraction. Whereas Omukai & Palla (2003) only show this evolution for shock outer boundary conditions, we show here that this oscillatory behavior occurs even in the extreme case of cold accretion, η = 0. Our results suggest that this characteristic evolutionary stage should emerge only with the moderately high accretion rates $\dot{M}_*\sim 10^{-2} \,M_\odot {\rm \,yr}^{-1}$, regardless of different thermal properties of the mass accretion.

The evolution of the interior structure differs significantly when the mass accretion is extremely cold. Figure 7 shows the evolution of the stellar interior structure for the extreme case of an input energy fraction η = 0 but with the same accretion rate (case a1-e0). The evolution of stellar structure deviates strongly from the case shown in Figure 2, particularly in the early stages for M_* ≲ 100 M_☉. Cold mass accretion leads to a much lower average entropy in the stellar interior, which results in a smaller stellar radius for M_* ≲ 30 M_☉. The central temperature is higher with the smaller radius for a given stellar mass, as shown in Figure 4(a); the central temperature exceeds 10⁶ K before the stellar mass reaches 10 M_☉. Deuterium burning begins much earlier than in the case with η = 0.1, even before the timescale inversion, due to decreasing opacity in the stellar interior (see Section 3.1.1 above). Figure 2 also shows that the protostar is fully convective in this early stage, as the energy produced in the nearly opaque interior is transported mostly via convection.

Zoom In Zoom Out Reset image size

Convective mixing transports the recently accreted deuterium inward to the stellar center, where nuclear fusion takes place. Central deuterium burning ceases soon after a radiative layer emerges in the stellar interior, which prevents convective mixing (this is a so-called radiative barrier; e.g., Palla & Stahler 1992). Deuterium burning resumes at the bottom of the surface convective layer, when the local temperature exceeds 1.5 × 10⁶ K. Deuterium shell burning injects heat near the surface, which causes the surface layer to abruptly inflate at M_* ≃ 50 M_☉.⁶ The later evolution for M_* ≳ 10³ M_☉, in particular, looks similar to the fiducial case shown in Figure 2. The evolution of accreting SMSs is thus insensitive to potential variations of the thermal efficiency of mass accretion, which is modeled with varying η in this paper.

3.1.5. Evolution in the H-R Diagram

Figure 8 shows the evolutionary tracks in the Hertzsprung–Russell (H-R) diagram for several representative cases. As shown in HOY12 for cases with rapid mass accretion $\dot{M}_*\ge 0.1 \,M_{\odot }$, the tracks are almost vertical as the stellar mass, luminosity, and radius increase for M_* ≳ 100 M_☉. The tracks for $\dot{M}_*= 1.0 \,M_{\odot }$ show a gradual increase of the effective temperature for M_* ≳ 10⁴ M_☉, finally reaching ≃ 10⁴ K at M_* = 10⁵ M_☉ (see Figure 1). Only during the early phases before the stellar mass (luminosity) exceeds 100 M_☉ (10⁶ L_☉), does the effective temperature depend on the input energy fraction η, i.e., lower η results in higher effective temperature. The evolutionary track for $\dot{M}_*= 10^{-2} \,M_\odot {\rm \,yr}^{-1}$ lies between the above tracks and the zero-age main sequence (ZAMS) line. Despite the apparent differences among these tracks, the stellar luminosities for a given mass are nearly the same in all cases for M_* ≳ 1000 M_☉, namely the Eddington value. It is clear that SMSs with extremely high mass accretion rates will have low effective temperatures compared to non-accreting or slowly accreting SMSs. This fact has a number of consequences on star and BH formation in the early universe, as discussed in Section 4 below.

Zoom In Zoom Out Reset image size

Next we study the pulsation stability of accreting SMSs and estimate the resulting mass-loss rates for unstable cases. Following our previous work (IHO13; Inayoshi et al. 2013b), we apply linear stability analysis to the stellar models with different masses. IHO13 show that an SMS is more unstable with higher accretion rates and at higher stellar masses. We thus focus on the case with the highest accretion rate $\dot{M}_*=1.0 \,M_\odot {\rm \,yr}^{-1}$ for examining whether the mass loss caused by the stellar pulsation can limit stellar growth via mass accretion. Since accreting SMSs are most unstable against radial perturbations for M_* ≲ 10³ M_☉ (IHO13), we also analyze the stability against radial perturbations. Our current results allow us to extend the analysis to the high-mass range of 10³ M_☉ ⩽ M_* ⩽ 10⁵ M_☉.

The basic procedure of our linear stability analysis is briefly summarized as follows (see Appendix in IHO13 for more details). First, we linearize the four time-dependent stellar structure equations (i.e., the continuity, momentum, energy, and energy transfer equations) and the Poisson equation, applying radial perturbations to the physical quantities. Here we consider the perturbations to be proportional to e^iσt, where σ = σ_R + iσ_I, and σ_R and |σ_I| are the frequency and growth (or damping) rate of the pulsation. For instance, one of the perturbed quantities is the radial displacement of a fluid element from the equilibrium position, ξ(r, t) ≡ ξ_r(r)e^iσt. Retaining first order results in a system of ordinary derivative equations for the perturbed quantities (e.g., ξ_r(r)) and σ. We solve these equations with appropriate boundary conditions, finding solutions for particular eigenvalues of σ. The radial profiles of the perturbations are given as the eigenfunctions of this system of equations.

Figure 9 shows the evolution of the growth/damping rates of pulsation as a function of the accreted mass. The mass-loss rates are estimated assuming that all the pulsation energy is converted into kinetic energy of the outflowing gas:

$$\begin{equation} \frac{\dot{M}_{\rm loss}}{2}v_{\rm esc}^2=2|\sigma _{\rm I}|E_{\rm W}, \end{equation} \tag{ 9 }$$

where v_esc = (2GM_*/R_*)^1/2 is the escape velocity from the stellar surface,

$$\begin{equation} E_{\rm W}=\frac{\sigma _{\rm R}^2}{2}\int ^{M_\ast }_0 |\xi _r|^2dM_r \end{equation} \tag{ 10 }$$

is the pulsation energy, and M_r the enclosed mass. Since pulsation energy is partially lost through radiative dissipation, our estimate of the mass-loss rate should be considered an upper limit (Papaloizou 1973a, 1973b).

Zoom In Zoom Out Reset image size

In agreement with our previous analysis, we find that the instability is excited by the κ mechanism due to the opacity bump in the He⁺ ionization layer. The accreting SMS is always unstable, except within the narrow mass range 10⁴ M_☉ ≲ M_* ≲ 3 × 10⁴ M_☉. The mass-loss rate rises with increasing mass for M_* ≲ 6 × 10³ M_☉. The derived mass-loss rate at M_* ∼ 10³ M_☉ is somewhat lower than the value shown in IHO13, which we attribute to the different numerical code employed in this paper. The current stellar model has a slightly higher temperature in the surface layers than the model analyzed in IOH13. The layers with the opacity bump, where the pulsation is excited, contain a smaller fraction of the stellar mass, which results in a lower mass-loss rate. The maximum mass-loss rate is 5 × 10⁻³ M_☉ yr⁻¹ at M_* ≃ 6 × 10³ M_☉. After this, the mass-loss rate hardly rises in spite of the growth rate increasing with stellar mass. This is primarily because the escape velocity increases with stellar mass (see Equation (9)), as the stellar expansion slows down for M_* ≳ 10⁴ M_☉ (Figure 1(a)). Since the mass-loss rates are much lower than the assumed accretion rate, we conclude that the pulsation instability does not prevent the formation of an SMS via rapid mass accretion.

When the stellar mass exceeds 6 × 10³ M_☉, the pulsation is stabilized and the pulsation-driven outflows are also suppressed. To understand this stabilization, we consider the work integral defined by

$$\begin{equation} W(M_r)=\frac{\pi }{\sigma _{\rm R}}\int ^{M_r}_0 \Re \left[ \frac{\delta T^*}{T} \left(\delta \epsilon - \frac{d}{dM_r}\delta L_{\rm rad} \right) \right] dM_r \,, \end{equation} \tag{ 11 }$$

where is the nuclear energy production rate per unit mass, L_rad the radiative luminosity, the symbols with δ represent the Lagrange perturbations, and ℜ[] denotes taking the real part of a quantity. The work integral physically represents the pulsation energy gained during one period within the mass coordinate M_r. The stellar pulsation is excited (damped) in the layer where dW/dM_r is positive (negative). If the work integral at the stellar surface W(M_*) is positive (negative), the protostar is unstable (stable).

Figure 10 shows the radial profile of the work integral near the surface of a 3 × 10⁴ M_☉ star. We see that the work integral increases outward within the He⁺ ionization layer, but decreases beyond it because of radiative diffusion. In the outermost layers with T ≲ 18, 000 K, where the radiative cooling time is shorter than the pulsation period, the work integral remains constant without excitation or damping (so-called the non-adiabatic region). In this case, the work integral at the stellar surface is negative and the protostar is stable. On the other hand, in the most unstable case when M_* = 6 × 10³ M_☉, radiative cooling in the outermost part is more efficient and the non-adiabatic region extends further inward. The transition point to the non-adiabatic region is located just outside of the He⁺ ionization layer. Therefore, the stellar pulsation excited in the He⁺ ionization layer does not suffer from damping via radiative diffusion outside the ionization layer. Such a profile of the work integral is displayed in Figure 3 of IHO13, where by simply extrapolating the available results for M_* < 10³ M_☉ without considering the damping effect outside the He⁺ ionization layer, we had speculated that the mass-loss rates would increase with stellar mass for M_* > 10³ M_☉. By contrast, our current results show that the pulsation instability is strongly suppressed for M_* ≫ 10³ M_☉ and mass loss is much lower than predicted by IHO13.

Zoom In Zoom Out Reset image size

Note that the pulsational instability considered above qualitatively differs from that of non-accreting massive primordial stars, which is caused by the -mechanism (e.g., Baraffe et al. 2001; Sonoi & Umeda 2012). The pulsation driven by the -mechanism could work when an SMS contracts and approaches the ZAMS stage after the mass accretion ceases (also see Section 4.3).

As mentioned in Section 1, for a more typical case of primordial star formation, stellar UV feedback becomes so strong that it is able to shut off the mass accretion onto the star (e.g., Hosokawa et al. 2011b, 2012b). By contrast, the stellar evolution calculations presented here predict that this is not the case when forming primordial SMSs via very rapid mass accretion. Figure 11 shows that with extremely high accretion rates $\dot{M}_*\ge 0.1 \,M_\odot {\rm \,yr}^{-1}$, the stellar output of hydrogen-ionizing photons remains very low even after the stellar mass exceeds 10⁴ M_☉. For accretion at 1.0 M_☉ yr⁻¹ the hydrogen-ionizing luminosity for M_* ≃ 10⁵ M_☉ is still comparable to that of M_* ≲ 100 M_☉ primordial ZAMS stars. This is a consequence of the low effective temperature resulting from the large stellar radius. For $\dot{M}_*= 0.1 \,M_\odot {\rm \,yr}^{-1}$, the accretion rate of hydrogen atoms is ≃ 3 × 10⁴⁸ s⁻¹. An H ii region cannot grow if the stellar hydrogen-ionizing photon emissivity is lower than this hydrogen accretion rate. Only for the extreme cases with low-entropy accretion η = 0, would the ionizing luminosity temporarily increase in an early stage when M_* < 100 M_☉, visible as a spike in Figure 11. This occurs shortly before the abrupt expansion of the star (Figures 6 and 7), when the stellar total luminosity also sharply rises. Even if this "flash" of ionizing radiation would occur, however, its duration of <100 yr is too short to significantly disturb the accretion flow. It is therefore unlikely that the mass accretion onto rapidly accreting SMSs is hindered by stellar UV feedback, at least for M_* ≲ 10⁵ M_☉.

Zoom In Zoom Out Reset image size

After the stellar mass exceeds 10⁵ M_☉, however, the stellar ionizing flux will continue to increase so that an H ii region might finally emerge. The UV feedback caused by the dynamical expansion of the H ii region would shut off the mass accretion as expected in the ordinary cases of the primordial star formation (e.g., McKee & Tan 2008; Hosokawa et al. 2011b). The maximum stellar mass by the UV feedback would be higher with more rapid mass accretion, as suggested in our radiation hydrodynamic numerical simulations (Hosokawa et al. 2011b, 2012b).

By contrast, stellar UV feedback at a lower accretion rate $\dot{M}_*= 10^{-2} \,M_\odot {\rm \,yr}^{-1}$ will slow and perhaps stop the accretion flow. Although this rate is too low to form an SMS, such moderately rapid accretion can occasionally be expected in the ordinary mode of primordial star formation (e.g., Hirano et al. 2013). Figures 8 and 11(a) show that, in this case, the effective temperature always lies between the ZAMS and the rapidly accreting ($\dot{M}_*\ge 0.1 \,M_\odot {\rm \,yr}^{-1}$) values. However, Figure 11(b) shows that the evolution of the stellar ionizing luminosity is almost the same as in the ZAMS case. This is because the larger radius (Figure 6(b)) compensates for the lower effective temperature, resulting in a comparable UV luminosity. Although stellar UV feedback is certainly expected in this case, it is somewhat weaker than for the ZAMS stars because of the lower average energy of ionizing photons.

Although the UV luminosity of an accreting SMS is much weaker than its ZAMS counterpart, the bolometric luminosity could exceed 10⁹ L_☉, almost linearly increasing with stellar mass. With such a high luminosity, individual SMSs could be detectable with future observational facilities such as the James Webb Space Telescope (JWST). Figure 12 demonstrates this, showing that bloated SMSs at z ∼ 10 are detectable by JWST with exposure times of 10–100 hr. The cosmological parameters for the standard ΛCDM (Komatsu et al. 2011) and black-body spectrum with the effective temperature T_eff are adopted here.

Zoom In Zoom Out Reset image size

Actually, the observational features of accreting SMSs are very similar to those of hypothetical supermassive dark stars (SMDSs), which are powered by energy production via annihilation of dark matter (DM) rather than nuclear fusion (e.g., Spolyar et al. 2008; Umeda et al. 2009). This is because SMDSs could have a very large radius R_* ≳ 10 AU and a low effective temperature T_eff ≲ 10⁴ K, similar to the accreting SMSs presented in this paper. Recent studies that examined the observational signatures and detectability of SMDSs (e.g., Zackrisson et al. 2010; Freese et al. 2010; Ilie et al. 2012), are thus applicable to accreting SMSs as well. For instance, SMDS and SMS colors obtained from multi-band survey data should be unusual compared to ordinary galaxies and active galactic nuclei because of their very low effective temperature (Zackrisson et al. 2010). Note that the evolution and apparent colors of SMDSs are strongly dependent on supply rates of annihilating DM (Freese et al. 2010). Once the SMDSs are no longer supplied with sufficient DM, they contract, increasing their effective temperature. Thus, it is only in extreme cases that a 10⁵ M_☉ SMDS would still have a low effective temperature of T_eff ≲ 10⁴ K. By contrast, rapid mass accretion keeps SMSs inflated, resulting in a low T_eff without DM annihilation. The brightest sources with unusually red colors might be a signature of rapidly accreting SMSs rather than of SMDSs.

The detection rate of accreting SMSs depends on their formation rate per unit redshift and typical lifetime. Agarwal et al. (2012) and Johnson et al. (2013a) show that the formation of the SMSs via rapid mass accretion should be more ubiquitous than previously thought. They show that, within the lifetime of 10⁶ yr, at least a few accreting SMSs should lie in the survey area of ∼100 arcmin². The exact lifetime of an SMS depends on the final stage of its evolution, which is affected by the mass accretion history (also see Section 4.3).

Our results suggest that rapid mass accretion onto SMSs drastically changes their stellar structure, causing significant bloating. This could also influence their ultimate fate. As described in Section 1, SMSs exceeding 10⁵ M_☉ are thought to directly collapse to form massive BHs via a GR instability. However, previous studies assume that the SMSs are polytropic stars with index n = 3, i.e., fully convective stars with homogeneous entropy distributions. Our calculations show that this is not the case for rapid mass accretion. As Figure 3(a) shows, in the stellar interior, the specific entropy increases outward until reaching the surface convective layer. The homogeneous convective core gradually extends after hydrogen burning is ignited at the stellar center.

Whereas GR effects are not taken into account in the current calculations, we have also performed a few test runs with a gravitational constant G_eff modified by post-Newton corrections

$$\begin{equation} G_{\rm eff} = G \left(1 + \frac{P}{\rho c^2} + \frac{4 \pi r^3 P}{M c^2} + \frac{2 {G M}}{c^2 r} \right) \end{equation} \tag{ 12 }$$

(e.g., Fuller et al. 1986) and see no significant differences from the results presented above. This is demonstrated in Figure 13, which shows the stability of SMSs accreting at the rates $\dot{M}_*= 1.0 \,M_\odot {\rm \,yr}^{-1}$ and 0.3 M_☉ yr⁻¹. The SMSs are still stable as far as we have followed their evolution. Although the entire stellar structure is not well approximated by a n = 3 polytrope, note that the unstable area in the figure is applicable to the central convective core.

Zoom In Zoom Out Reset image size

Extrapolating the current results, we expect that the star becomes unstable when the stellar mass reaches a few × 10⁵ M_☉ for $\dot{M}_*= 1.0 \,M_\odot {\rm \,yr}^{-1}$. We expect that not the entire star but only its inner core might collapse due to the GR instability. This is because the outer part broadly extends with high specific entropy obtained with very rapid mass accretion. If a BH forms as a result of the collapse, an accretion disk should also form and grow as the surrounding gas falls toward the center. A quasi-star powered by the energy production in the accretion disk instead of nuclear fusion (Begelman et al. 2008) might finally emerge. The figure also suggests that the GR instability could work for a slightly lower mass assuming a lower accretion rate, because the mass distribution in the interior is more centrally condensed (also see Section 3.1.2). With significantly lower accretion rates $\dot{M}_*\lesssim 0.1 \,M_\odot {\rm \,yr}^{-1}$, however, the star would exhaust its nuclear fuel and collapse directly to a BH before the GR instability sets in.

Note that the above scenario assumes ongoing mass accretion when the GR instability sets in. As discussed in Section 4.1, however, stellar UV feedback might halt the mass accretion before the GR instability sets in. Without mass accretion, the bloated SMS should contract on a KH timescale and finally become a ZAMS star, which has a homogeneous entropy distribution because of efficient convective mixing. For such n ≃ 3 polytrope stars, the GR instability should cause most of the stellar matter to collapse to a massive BH (e.g., Fricke 1973; Fuller et al. 1986). In this case, the stellar lifetime might be somewhat longer than in the above case, where the SMS collapses before the mass accretion ceases.

Even if a part of the star becomes unstable, this does not necessarily mean that a BH forms from the collapse. It is conceivable that a very energetic supernova explosion results from runaway nuclear fusion. Although this is less likely with the lower metallicities (e.g., Montero et al. 2012), recent studies suggest that this might occur in the early universe (e.g., Whalen et al. 2013a; Johnson et al. 2013b). Moreover, the evolution of an SMS should also depend on other parameters, e.g., the interior rotation and magnetic fields. Indeed, Yoon et al. (2012) show that the ultimate fates of M_* ⩽ 10³ M_☉ primordial stars vary with these parameters. The dynamics of the final collapse should also depend on the stellar rotation. Reisswig et al. (2013) show that a massive binary BH might form through the collapse of a rapidly rotating SMS. Thus, the final fates of SMSs should be studied further under consideration of all the potential effects of very rapid mass accretion, truncation of accretion, runaway nuclear fusion, the GR instability, and stellar rotation and magnetic fields.

Besides following the evolution of individual accreting SMSs and the resulting formation of massive BHs, it is also important to consider how often this occurs in the early universe. As described in Section 1, SMSs should form in massive dark halos where H₂ molecules are destroyed. Previous studies have mostly considered that H₂ molecules are destroyed by photodissociating radiation from neighboring stars, showing that the environments for forming the SMSs are extremely rare (e.g., Dijkstra et al. 2008). Regarding this issue, Agarwal et al. (2012) demonstrate that the occurrence rate of SMS formation should be elevated when considering the soft UV radiation from Pop II stars, and that the number density of massive BHs at z ≳ 6 expected from current observations (e.g., Treister et al. 2011; Mortlock et al. 2011) could agree with this scenario. However, it is still uncertain whether H₂ photodissociation is the dominant process which causes the SMS formation or not. Inayoshi & Omukai (2012), for instance, show that H₂ molecules are easily collisionally dissociated in dense shocks emerging in colliding cold accretion flows at the galaxy formation epoch. This could be more common than the strong photodissociating radiation field and should be examined in detail in future studies.

Recent numerical simulations are beginning to reveal the dynamical evolution in an early stage of SMS formation, whereby a massive circumstellar disk forms as gas falls toward the center of a gas cloud (e.g., Regan & Haehnelt 2009; Latif et al. 2013a; Choi et al. 2013). However, several assumptions of the simulations are somewhat arbitrary for enabling SMS formation, e.g., turning off H₂ molecular cooling or elevating the photodissociating background field by hand. We need more realistic simulations beginning with the initial conditions naturally provided by modern cosmology.

These previous simulations do not show that an SMS with M_* ≫ 1000 M_☉ actually forms via rapid mass accretion. We have shown, for the first time, the formation of a M_* ≃ 10⁵ M_☉ SMS with the large accretion rates predicted by previous studies. Our results also suggest that stellar radiative feedback could be effective after the stellar mass exceeds 10⁵ M_☉, a prediction which needs to be examined with more detailed simulations. It would be necessary to consistently follow the growth of the supergiant protostar and the evolution of the accretion flow, by including stellar radiative feedback (e.g., Hosokawa et al. 2011b).