EVOLUTION OF MASSIVE PROTOSTARS WITH HIGH ACCRETION RATES

We employ the method of calculation developed by SPS86 and PS91. In this section, only the outline of modeling is overviewed. A detailed description of the method is presented in Appendix A.

As shown in Figure 1, the whole system can be divided into two parts of different nature, i.e., an accreting envelope and a hydrostatic core. This hydrostatic core is also called a protostar since it eventually grows into a star by accretion. The outer part of the accreting envelope contains dust grains, while the warm (≳2000K) part near the protostar is dust-free as a result of their evaporation. Despite the absence of the dust, the innermost dense part of the envelope becomes optically thick to gas continuum opacity and the photosphere appears outside the accretion shock front in the case of a high accretion rate. Such an inner envelope is called the radiative precursor. We here consider only the protostar and the radiative precursor under a given constant accretion rate (see Figure 1). The dusty outer envelope is not included in our formulation, although stellar feedback exerted on the envelope may be important in the last stage of accretion: the mass accretion can be terminated finally by radiation pressure or other protostellar feedback processes (e.g., stellar wind) exerted on the dusty envelope. For the present, we presume that the protostar continues to grow at a given mass accretion rate and look for solutions of protostars with steady-state accretion. Discussion on feedback that halts the accretion is deferred to Section 4.

Zoom In Zoom Out Reset image size

We calculate protostellar evolution by constructing a time sequence of quasi-steady structures of the protostar and accreting envelope. We solve the stellar structure equations for the protostar. For the accreting envelope, we adopt different treatments depending on whether or not it is opaque to the gas opacity. If the flow remains optically thin and no radiative precursor exists, the free-fall is assumed to be outside the star. For the opaque flow, on the other hand, we solve the structure of the radiative precursor by using the equations for a steady-state flow. At the outer boundary, which is taken to be at the photosphere, the flow is assumed to be in free fall. After solving the protostar and accreting envelope individually, these solutions are connected at the accretion shock front by the radiative shock condition (e.g., SST80a). The shooting method is adopted for solving radial structure. This procedure is repeated until the required boundary conditions are satisfied.

We start the calculation from a very small protostar, typically M_*,0 = 0.01M_☉. The initial models are constructed following SST80b (see Appendix A.2 for details). Although this choice is rather arbitrary, the star converges immediately to a certain structure appropriate for accretion at a given rate. Specific conditions of the initial models do not affect the evolution thereafter. The evolution is followed by increasing the stellar mass owing to accretion step by step. This procedure is repeated until the star reaches the ZAMS phase after the onset of hydrogen burning. In a few runs with very high accretion rates, however, steady accretion becomes impossible before arrival at the ZAMS and the calculation is terminated at this moment.

The calculated runs and their input parameters are listed in Table 1. Cases with a wide range of accretion rates are studied, starting from 10⁻⁶ M_☉ yr^-1 up to 6 × 10⁻³ M_☉ yr^-1. The lowest values of 10⁻⁶ to 10⁻⁵ M_☉ yr^-1 are typical for low-mass star formation, while high rates >10⁻⁴ M_☉ yr^-1 are those envisaged in the accretion scenario of massive star formation. The initial deuterium abundance of [D/H] = 2.5 × 10⁻⁵ is adopted as a fiducial value, following previous work (e.g., SST80a, PS91) for comparison. To assess the role of deuterium burning, we also calculate "noD" runs, where the deuterium is absent. Most of the runs are for the solar metallicity Z_☉(= 0.02). Two runs with 10⁻³ and 10⁻⁵ M_☉ yr^-1 are calculated for the metal-free gas (runs MD3-z0 and MD5-z0) to see the effects of different metallicities. For 10⁻⁵ M_☉ yr^-1, variations in deuterium abundances (runs MD5-dh1 and MD5-dh3) and in initial models (MD5-ps91) are studied and presented in Appendix B for comparison with previous calculations.

Table 1. Calculated Runs and Input Parameters

Run	$\dot{M}_\ast (M_{\odot }\,{\rm yr}^{-1})$a	[D/H](10−5)b	Z c	M*,0(M☉)d	R*,0(R☉)e	Referencef
MD6x3	6 × 10−3	2.5	0.02	0.5	43.9	Section 3.5
MD4x3	4 × 10−3	2.5	0.02	0.3	33.0	Section 3.5
MD3x3	3 × 10−3	2.5	0.02	0.2	27.8	Section 3.5
MD3x3-z0	3 × 10−3	2.5	0.0	0.2	26.4	Section 3.5
MD3	10−3	2.5	0.02	0.05	15.5	Sections 3.1, 3.3 , 3.4, 3.5
MD3-noD	10−3	0.0	0.02	0.05	15.5	Section 3.3
MD3-z0	10−3	2.5	0.0	0.05	13.0	Section 3.4
MD4	10−4	2.5	0.02	0.01	7.9	Section 3.3
MD4-noD	10−4	0.0	0.02	0.01	7.9	Section 3.3
MD5	10−5	2.5	0.02	0.01	3.7	Section 3.2, 3.3 , 3.4
MD5-noD	10−5	0.0	0.02	0.01	3.7	Section 3.2, 3.3
MD5-dh1	10−5	1.0	0.02	0.01	3.7	Appendix B.1
MD5-dh3	10−5	3.0	0.02	0.01	3.7	Appendix B.1
MD5-z0	10−5	2.5	0.0	0.01	2.9	Section 3.4
MD5-ps91g	10−5	2.5	0.02	1.0	4.2	Appendix B.2
MD6	10−6	2.5	0.02	0.01	1.6	Section 3.3
MD6-noD	10−6	0.0	0.02	0.01	1.6	Section 3.3

Note. ^amass accretion rate. ^binitial number fractional abundance of deuterium. ^cmetallicity. ^dmass of initial core model. ^eradius of initial core model. ^fsubsections where numerical results of each run are presented. ^ginitial model is the same as PS91.

Download table as:  ASCIITypeset image

Initial stellar mass in each run is taken to be sufficiently small, typically M_*,0 = 0.01 M_☉. For high accretion rates $\dot{M}_\ast \ge 10^{-3}\,M_{\odot }\, {\rm yr^{-1}}$, somewhat more massive initial stars (see Table 1) are used since convergence of calculation was not achieved for lower mass ones. The radii of initial models are also listed in Table 1.

First, we see a protostar growing at a high accretion rate $\dot{M}_\ast = 10^{-3}\, M_{\odot }\, {\rm yr^{-1}}$ (run MD3), which is relevant to massive star formation. The upper panel of Figure 2 shows evolution of the protostellar radius (thick solid) and the interior structure as a function of the protostellar mass. The mass increases monotonically with time t, $M_{\ast }=M_{\ast,0}+\dot{M}_\ast t$, and can be regarded as a time coordinate. The entire evolution can be divided into the following four phases according to their characteristic features: (I) the adiabatic accretion (for M_* ≲ 6 M_☉), (II) swelling (6 M_☉ ≲ M_* ≲ 10 M_☉), (III) KH contraction (10 M_☉ ≲ M_* ≲ 30 M_☉), and (IV) main-sequence accretion (M_* ≳ 30 M_☉) phases. Below, we see the protostellar evolution in each phase.

Zoom In Zoom Out Reset image size

Adiabatic Accretion Phase. Upper panels of Figure 3 present the evolution of entropy profile within the protostar. This shows that the entropy profiles at different epochs completely overlap: at a fixed mass element, the specific entropy remains constant. The entropy is originally generated at the accretion shock front, and embedded into the stellar interior. Thus, during this earliest phase, the entropy profile inside the star just traces the history of entropy at the postshock point. Owing to a high value of opacity, radiative heat transport is inefficient in the interior of the protostar. This makes the interior luminosity small, then t_KH ≫ t_acc in this phase (see Figure 4 bottom panel, below). Once the accreted matter settles into the opaque interior, the specific entropy is just conserved.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The snapshots of the luminosity profile are also presented in lower panels of Figure 3. Near the surface is a spiky structure where the gradient ∂L/∂M is very large owing to a sharp decrease in opacity there. Although a significant entropy loss occurs in this luminosity spike as indicated by the relation (∂s/∂t)_M ≃ −1/T(∂L/∂M)_t (from Equation A3), this layer is extremely thin and the time for a mass element to pass though it is so short that the entropy hardly changes. This can be confirmed by comparing the two timescales; the local accretion time

$$\begin{equation} t_{\rm acc,s} (m) = \frac{m}{\dot{M}_\ast }, \end{equation} \tag{ 8 }$$

and local cooling time

$$\begin{equation} t_{\rm cool,s} (m) = \frac{m \delta s}{\int _0^m \frac{1}{T} \frac{\partial L}{\partial m} dm^{\prime }}, \end{equation} \tag{ 9 }$$

where m = M_* − M is the depth (in mass coordinate) of a mass shell from the accretion shock, and δs denotes arbitrary small entropy change. The former, t_acc,s, is the time for a thin shell of mass m at the surface to be replaced by the newly accreted material, while the latter, t_cool,s, is the time for the same shell to lose the entropy mδs by outward radiation. Comparison between these timescales in the settling layer is presented in Figure 5. where δs = 0.1k_B/m_H is adopted for numerical evaluation. This indicates that t_acc,s is always shorter than t_cool,s: the accreted material is swiftly embedded in the interior before losing the entropy δs by radiation. Therefore, the adiabaticity remains valid in spite of the spike in luminosity profile. Note that the short t_acc,s is a result of the high accretion rate. We will see that the situation changes for much lower accretion rates in Sections 3.2 and 3.3 below.

Zoom In Zoom Out Reset image size

With increasing protostellar mass M_*, the accretion shock strengthens and the postshock entropy increases. Then, the interior entropy increases as well. This causes a gradual expansion of the stellar radius in the adiabatic accretion phase (Figure 2, upper panel). Using the typical density and pressure within a star of mass M_* and radius R_* (e.g., Cox & Giuli 1968);

$$\begin{equation} \rho \sim \frac{M_{\ast }}{R_{\ast }^3}, \qquad P \sim G \frac{M_{\ast }^2}{R_{\ast }^4}. \end{equation} \tag{ 10 }$$

and the expression for specific entropy of ideal monatomic gas

$$\begin{equation} s = \frac{3 \cal {R}}{2 \mu } \ln \left(\frac{P}{\rho ^{5/3}} \right) + const, \end{equation} \tag{ 11 }$$

where ${\cal R}$ is the gas constant and μ is the mean molecular weight, the stellar radius and entropy are related as (Stahler 1988);

$$\begin{equation} R_{\ast } \propto M_{\ast }^{-1/3} \exp \left[ \frac{2 \mu }{3 {\cal R}} s \right], \end{equation} \tag{ 12 }$$

i.e., the stellar radius is larger for the higher entropy within the star. SPS86 have derived the mass–radius relation for protostars in the adiabatic accretion phase;

$$\begin{equation} R_{\ast } \simeq 26\,R_{\odot } \left(\frac{M_{\ast }}{M_{\odot }} \right)^{0.27} \left(\frac{\dot{M}_\ast }{10^{-3}\,M_{\odot }\,{\rm yr}^{-1}} \right)^{0.41}, \end{equation} \tag{ 13 }$$

under the condition that opacity in the radiative precursor is dominated by H⁻ bound-free absorption. As this condition holds in our case, our calculated mass–radius relation is in good agreement with Equation (13). Equation (13) suggests the large radius (>10 R_☉) of the rapidly accreting protostar.

The interior of the protostar remains radiative throughout this phase (Figure 2, upper panel): all the energy transport is via radiation. This is in high contrast with low $\dot{M}_\ast$ (∼10⁻⁵ M_☉ yr^-1) cases, where most of the interior becomes convective owing to deuterium burning for M_* ≳ 0.4 M_☉ (see Section 3.2 below). This can be attributed to a difference in the interior temperature. From Equation (10), typical temperature within the star is

$$\begin{equation} T = \frac{\mu }{\cal {R}} \frac{P}{\rho } \sim \frac{G}{\cal {R}} \frac{\mu M_{\ast }}{R_{\ast }}. \end{equation} \tag{ 14 }$$

Therefore, the large stellar radius leads to the low temperature in the stellar interior. In fact, the maximum temperature in this phase does not exceed the threshold for deuterium burning (lower panel of Figure 2).

Once the accreted matter has passed through the outermost layer with the luminosity spike, the luminosity gradient becomes much milder. The luminosity has a maximum at some radius: outside the luminosity maximum, the gradient ∂L/∂M < 0, while ∂L/∂M > 0 inside. This means that heat is removed from the deep interior and absorbed in the outer layer. This entropy transfer, however, remains small and does not modify the entropy distribution significantly during this phase. Efficiency of the outward entropy transfer is related to the value of opacity. In most of the stellar interior, major sources of opacity are bound–bound and bound-free transitions of heavy elements, which approximately obey Kramers' law, κ ∝ ρT^−3.5 (Hayashi et al. 1962; Clayton 1968). With the increase in stellar mass and then the interior temperature, the opacity decreases. This results in steady increase of the outward heat flux with evolution. In fact, as shown in the lower panels of Figure 3, the amplitude of the luminosity increases with the growth of the protostar. Also the maximum luminosity within the star L_max increases as a power-law function of M_* as seen in the middle panel of Figure 4. This dependence can be understood as follows: for a star with radiative stratification and Kramers' opacity, the luminosity scales as L_rad ∝ M^11/2_*R^−1/2_* (e.g., Cox & Giuli 1968). We have confirmed that L_max roughly obeys,

$$\begin{equation} L_{\rm max} \simeq 0.2\,L_{\odot } \left(\frac{M_{\ast }}{M_{\odot }} \right)^{11/2} \left(\frac{R_{\ast }}{R_{\odot }} \right)^{-1/2} \end{equation} \tag{ 15 }$$

in our calculations. Using Equations (13) and (15), we obtain L_max ∝ M^5.4_* for a constant accretion rate. When opacity becomes sufficiently small, the radiative heat transport becomes significant even in the deep interior. Consequently, the entropy distribution and thus stellar structure are drastically altered, which leads to the subsequent swelling phase.

Swelling Phase. At the mass of about 6 M_☉, the protostar begins to swell up dramatically. The stellar radius suddenly increases by a factor of about three and eventually exceeds 100 R_☉. This swelling is caused by redistribution of entropy within the star. As explained above, the outward heat flow continues to increase, and ceases to be negligible. Figure 3 shows that the entropy distribution changes significantly with time, i.e., with M_*: the entropy decreases in the deep interior and increases in the outer layer. The extent of the inner entropy-losing region becomes larger and larger, and approaches the surface, as more of the interior becomes less opaque.

As seen in the bottom panels of Figure 3, the peak of the luminosity, which is the boundary between the entropy-losing interior and the outer absorbing layer, gradually moves toward the surface with increasing height. This outward propagation of the luminosity peak is called the luminosity wave (SPS86). Only a thin outlying layer absorbs the entropy transported from the deep interior. This rapid increase of the entropy near the surface causes swelling of the radius. Actually, only a small fraction of mass near the surface contributes to this swelling. For example, at the maximum expansion of the radius R_* ≃ 140 R_☉ at M_* ≃ 10 M_☉, only 0.03% of the total mass is contained in the outer layers of R_* > 70 R_☉.

Note that the swelling is not due to the shell burning of deuterium. Palla & Stahler (1990) have presented that a similar swelling occuring in intermediate-mass protostars at the low accretion rate of $\dot{M}_\ast = 10^{-5}\, M_{\odot }\, {\rm yr^{-1}}$, and concluded that this is caused by the shell burning of deuterium. In our case, however, deuterium burning is not important for the swelling. Even in the run without deuterium burning (MD3-noD), the swelling is caused by redistribution of entropy in the same fashion (also see Section 3.3 below). During this phase, the deuterium burning occurs in the inner region (Figure 2, upper panel). Despite the active deuterium burning, most of the stellar interior remains radiative. Because of the higher entropy and thus lower density in the star, the opacity is lower in our case when the active deuterium burning begins. This allows efficient radiative transport of the generated entropy (Stahler 1988). The middle panel of Figure 4 shows that the maximum radiative luminosity L_max always exceeds that by the deuterium burning L_D. This indicates a large capacity of radiative heat transport. Also a thin convective layer near the surface (Figure 2) is not brought about by the deuterium burning. Although this layer receives the transported entropy from inside, the higher opacity due to ionization prevents the outward heat flow. This causes a negative entropy gradient and then convection near the surface.

Kelvin–Helmholtz Contraction. With the approach of the luminosity wave to the surface, an increasing amount of energy flux escapes from the star without being absorbed inside. This results in a significant increase in the interior luminosity L_* for M_* > 5 M_☉ (Figure 4) and thus the shorter KH timescale t_KH (∝1/L_*) of the star. The arrival of the luminosity wave at the stellar surface marks the moment that all parts of the star begin to lose heat. Almost simultaneously, t_KH becomes as short as the accretion timescale t_acc (Figure 4). While maintaining the virial equilibrium against the radiative energy loss, the protostar turns to contraction. This is the so-called KH contraction phase. The epoch of the radial turn-around is roughly estimated by equating t_acc and t_KH,lmax ≡ GM²_*/R_*L_max, where L_max is given by Equation (15),

$$\begin{equation} M_{\rm *,rmax} = 10.8\,M_{\odot } \left(\frac{\dot{M}_\ast }{10^{-3}\,M_{\odot }\,{\rm yr}^{-1}} \right)^{2/9}. \end{equation} \tag{ 16 }$$

In the above, we have omitted a weak dependence on R_* and just substituted R_* ∼ 10 R_☉ for simplicity. During the contraction of the protostar, the KH timescale remains slightly shorter than the accretion timescale. Since the ratio t_KH/t_acc, or equivalently L_*/L_max, is less than unity, the interior luminosity exceeds the accretion luminosity and becomes the dominant source of the total luminosity of the protostar (upper panel of Figure 4).

Since the stellar interior remains radiative in spite of the deuterium burning, the accreted deuterium is not transported to the deep interior by convection. Deuterium in the deep interior is thus soon consumed. Thereafter the deuterium burning continues in a shell-like outer layer (Figure 2, upper panel). This surface burning proceeds at a rate slightly exceeding the steady burning one (e.g., Palla & Stahler 1990),

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

where δ_D is the energy available from the deuterium burning per unit gas mass. Consequently, the mass-averaged deuterium concentration f_d,av significantly decreases during the KH contraction (Figure 2, lower panel).

Arrival at a Main-sequence Star Phase. With the KH contraction, the interior temperature increases gradually. At M_* ≃ 20 M_☉ the maximum temperature reaches 10⁷K and the nuclear fusion of hydrogen begins (Figure 2). Although the hydrogen burning is initially dominated by the pp-chain reactions, energy generation by the CN-cycle reactions immediately overcome this. The energy production by the CN-cycle compensates radiative loss from the surface at M_* ≃ 30 M_☉ (Figure 4, middle), where the KH contraction terminates. A convective core emerges owing to the rapid entropy generation near the center. The stellar radius increase after that, obeying the mass–radius relation of main-sequence stars.

Evolution of the Radiative Precursor. Figure 6 shows the evolution of protostellar (R_*) and photospheric (R_ph) radii. Except in a short duration, the photosphere is formed outside the stellar surface: the radiative precursor persists throughout most of the evolution. In the adiabatic accretion phase, the photospheric radius is slightly outside the protostellar radius and increases gradually with the relation R_ph ≃ 1.4 R_*, which is derived analytically by SPS86. Although the precursor temporarily disappears in the swelling phase around M_* ≃ 10 M_☉, it emerges again in the subsequent KH contraction phase. In the KH contraction phase, the precursor extends spatially and becomes more optically thick. Due to the extreme temperature sensitivity of H⁻ bound-free absorption opacity, the region in the accreting envelope with temperature ≳6000 K becomes optically thick. Thus, the photospheric temperature T_ph remains at 6000–7000 K throughout the evolution (Figure 6, lower panel). Since the photospheric radius is related to the total luminosity by

$$\begin{equation} R_{\rm ph}=\left(L_{\rm tot}\big/4\pi \sigma T_{\rm ph}^4 \right)^{1/2}, \end{equation} \tag{ 18 }$$

for the constant T_ph, the rapid increase in the luminosity in the contraction phase causes the expansion of the photosphere.

Zoom In Zoom Out Reset image size

Next, we revisit the evolution of an accreting protostar under the lower accretion rate of $\dot{M}_\ast = 10^{-5}\,M_{\odot }\, {\rm yr^{-1}}$ (run MD5). Although this case has been extensively studied by previous authors, a brief presentation should be helpful to underline the effects of different accretion rates on protostellar evolution. Different properties of the protostar even at the same accretion rate owing to updates from the previous works, e.g, initial models and opacity tables, will also be described. For a thorough comparison between our and some previous calculations, see Appendix B.

The lower accretion rate means the longer accretion timescale $t_{\rm acc}\propto 1/\dot{M}_\ast$: even at the same protostellar mass M_*, the evolutionary timescale is longer in the case of low $\dot{M}_\ast$. Therefore, the protostar has ample time to lose heat before gaining more mass. This results in the lower entropy and thus a smaller radius at the same protostellar mass in the low $\dot{M}_\ast$ case as shown in the upper panel of Figure 7. Although with the smaller value, the overall evolutionary features of the high and low $\dot{M}_\ast$ radii are similar.

Zoom In Zoom Out Reset image size

Again, the protostellar evolution can be divided into four characteristic stages, i.e., (I) convection (M_* ≲ 3 M_☉), (II) swelling (3 M_☉ ≲ M_* ≲ 4 M_☉), (III) KH contraction (4 M_☉ ≲ M_* ≲ 7 M_☉), and (IV) main-sequence accretion (M_* ≳ 7 M_☉) phases. Note that the first phase is now the convection phase instead of the adiabatic accretion. In the low $\dot{M}_\ast$ case, the deuterium burning drastically affects the protostellar evolution in early phases (e.g., Stahler 1988). Since the evolution in phases (III) and (IV) is similar to the previous high-$\dot{M}_\ast$ case (Section 3.1), we emphasize the two early phases (I) and (II) in the following subsection.

Convection Phase. In the low $\dot{M}_\ast$ case, the D burning begins in an earlier phase and has a more profound effect on the evolution than in the high-$\dot{M}_\ast$ case (run MD3). As a result of the smaller radius, the temperature, as well as density, is higher in the low $\dot{M}_\ast$ case (see lower panels of Figure 2 and Figure 7). The maximum temperature reaches 10⁶ K as early as M_* ≃ 0.2 M_☉ and active deuterium burning begins subsequently (Figure 7), while in the high-$\dot{M}_\ast$ case the D ignition does not occur until M_* ≃ 6 M_☉ (Figure 2). Unlike in the high-$\dot{M}_\ast$ case, the deuterium burning makes most of the interior convective. This is a result of high density and thus high opacity at the ignition of deuterium burning, occurring at a fixed temperature of ≃10⁶ K. Since high opacity hinders the radiative heat transport, entropy generated by deuterium burning cannot be transported radiatively. This necessarily gives rise to convection.

Figure 7 shows that the extent of the convective layers experiences a complicated history. At M_* ≃ 0.3 M_☉, a convective layer appears for the first time slightly off-center. Inside the convective layer, entropy is homogenized and its value increases owing to the D burning (Figure 8). With this increased entropy, the outer layer is gradually incorporated into the convective region. Fresh deuterium is supplied from the newly incorporated layer and spreads over the convective region, and finally is consumed by nuclear burning at the bottom of the region. However, the expansion of the convective region cannot be maintained if the supply of deuterium becomes insufficient. Without the fuel, the deuterium concentration in the convective region f_d,cv drops and the nuclear burning ceases: the convective region disappears. This occurs at M_* ≃ 0.6 M_☉ (Figure 7, upper). At M_* ≃ 0.7 M_☉, another episode of deuterium burning begins just outside the first convective region, where fresh deuterium still remains. The second convective region develops there and immediately reaches the stellar surface. Although this complicated evolution of the convective regions is different from that reported in SST80b, this difference can be attributable to high sensitivity of the evolution on adopted initial deuterium abundance. We discuss this in greater details in Appendix B.1.

Zoom In Zoom Out Reset image size

Since the energy generation by D burning is very sensitive to the temperature, temperature remains constant at ≃1.5 × 10⁶ K during the active deuterium burning (Figure 7, lower panel). This is the so-called thermostat effect (e.g., Stahler 1988). If this works, M_*/R_* ∝ T ≃ const. from Equation (14): the protostellar radius R_* increases as ∝M_*. This linear increase in R_* with M_* can be observed in the range 0.7 M_☉ ≲ M_* ≲ 1 M_☉ (Figure 7, upper panel).

When the deuterium concentration in the convective region drops to a few percent, the thermostat effect ceases to work. The maximum temperature starts increasing again and the radial expansion slows down after that. Owing to the high temperature, opacity decreases in the deep interior with increasing stellar mass. This allows efficient radiative energy transport and renders a large portion of the interior radiative. The radiative core gradually extends to the outer layer, and the convective layer moves outward. Since fresh deuterium is no longer supplied to the radiative interior by the mixing, deuterium is soon exhausted there. The expansion of the radiative interior occurs in a later phase in the calculation of PS91. This difference can be attributed to different initial models: whereas PS91 assumed a fully convective star of 1 M_☉ as an initial model, we started calculation from the initial mass M_*,0 = 0.01 M_☉. At M_* = 1 M_☉, the protostar is not fully convective in our case. The outer layers of M > 0.25 M_☉ are convective, where the D burning occurs at the bottom, while the remaining inner core is radiative at this moment. We confirmed that, starting from the same initial model, the evolution becomes similar to that of PS91 as discussed in Appendix B.2.

Because of the early ignition of deuterium burning, entropy and luminosity profiles in this phase look rather different from those in the high $\dot{M}_\ast$ case (see Figure 3 and Figure 8). For easier comparison, let us see the case without deuterium burning, where similar profiles are in fact reproduced (Figure 9). However, one important difference exists, the spiky structure in the entropy profiles near the surface. This shows significant entropy loss near the surface in the low $\dot{M}_\ast$ case, which leads to the low entropy within the star. The reason can be seen clearly in Figure 10, where the comparison between the local accretion timescale t_acc,s and cooling timescale t_cool,s (Equation (8) and (9)) is presented. In contrast to the case of high $\dot{M}_\ast$, in most of the settling layer, t_acc,s is shorter than t_cool,s: accreted matter loses entropy radiatively near the surface before settling into the adiabatic interior.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Swelling Phase. In the short period of 3 M_☉ ≲ M_* ≲ 4 M_☉, the stellar radius jumps up by a factor of two. At the same time, the whole interior becomes radiative (Figure 7, upper). Although deuterium burning continues near the surface, it does not drive convection any more.

The cause of this swelling has been attributed to the shell burning of deuterium, which also occurs in our calculation (Palla & Stahler 1990). However, we conclude that the main driver of the swelling is not the D shell burning, but the outward entropy transportation as in the case of the high accretion rate. The reasons are as follows. First, evolutionary features of entropy and luminosity profiles in this phase are similar to those in the swelling phase in the high $\dot{M}_\ast$ case (run MD3; see Figures 3 and 8): in both cases, entropy decreases in the deep interior and increases significantly near the stellar surface; simultaneously, the peak of the luminosity distribution moves to the surface as the luminosity wave. Second, even in the run without deuterium burning (run MD5-noD), this swelling occurs as shown in Figure 11: the protostar swells up at M_* ≃ 3 M_☉.

Zoom In Zoom Out Reset image size

The interior luminosity L_* jumps up at M_* ≃ 0.7 M_☉ with the arrival of the convective region at the surface as a result of high efficiency of convective heat transfer (Figure 12, upper panel). Although this rise in L_* decreases t_KH, t_KH remains much longer than t_acc (Figure 12, lower panel). It is not until another sudden rise in L_* at ≃3 − 4 M_☉ due to the arrival of the luminosity wave that t_KH becomes as short as t_acc. The maximum radius is attained at this epoch since the swelling is also caused by the luminosity wave, as seen above. Then, the analytical argument leading to Equation (16) also holds in this case: for $\dot{M}_\ast = 10^{-5}\, M_{\odot }\, {\rm yr^{-1}}$, this leads to the maximum radius at M_*,rmax ≃ 3.9 M_☉, which agrees with our numerical result.

Zoom In Zoom Out Reset image size

Simultaneous returning to the radiative structure with the swelling has been explained by Palla & Stahler (1990). The luminosity from the steady deuterium burning is comparable to the accretion luminosity by chance:

$$\begin{equation} L_{\rm D, st} = \dot{M}_\ast \delta _{\rm D} \sim \dot{M}_\ast \frac{G M_{\ast }}{R_{\ast }} = L_{\rm acc}. \end{equation} \tag{ 19 }$$

The epoch of the radial turn-around obeys Equation (16), which is derived by equating t_acc and t_KH,lmax = GM²_*/R_*L_max. This means that L_max becomes comparable to L_acc and thus to L_D,st in the swelling phase. Hence, the radiative luminosity now has a capacity to transport the energy generated by deuterium burning: L_max ∼ L_D,st. The deuterium burning no longer gives rise to any convective layer after the swelling.

Later Phases. Subsequently, the protostar enters the KH contraction phase. The evolution thereafter is similar to the previous high-$\dot{M}_\ast$ case. The maximum temperature within the star increases with the contraction, and hydrogen burning begins at M_* ≃ 7 M_☉.

Evolution of the Radiative Precursor. Figure 13 shows that the accretion flow is optically thin except in the early phase of M_* < 1 M_☉ at the low accretion rate of 10⁻⁵ M_☉ yr^-1. With the expansion at the active deuterium burning, the stellar surface emerges above the photosphere at M_* ≃ 1 M_☉ and the radiative precursor disappears thereafter. Without the precursor, radiation streams out from the stellar surface without hindrance with effective temperature rising remarkably in the KH contraction phase, i.e., for M_* > 4 M_☉. This is in stark contrast to the case of high accretion rate 10⁻³ M_☉ yr^-1, where the effective temperature remains several thousand K throughout the evolution (Figure 6).

Zoom In Zoom Out Reset image size

Here, we summarize the effects of the different accretion rates on the protostellar evolution. The mass–radius relations of protostars with accretion rates $\dot{M}_\ast =10^{-6}$, 10⁻⁵, 10⁻⁴, and 10⁻³ M_☉ yr^-1 are shown in the upper panel of Figure 14.

Zoom In Zoom Out Reset image size

First, the protostellar radius is larger for the higher accretion rate at the same protostellar mass, as discussed in Section 3.1. For example, at M_* = 1 M_☉, the radius R_* ≃ 25 R_☉ for 10⁻³ M_☉ yr^-1, while it is only 2.5 R_☉ for 10⁻⁶ M_☉ yr^-1. For easier comparison, the mass–radius relations in the no-deuterium cases are also shown by thin lines in Figure 14. The deuterium burning enhances the stellar radius, in particular, in low $\dot{M}_\ast$ cases (see in this section below). Nevertheless, the trend of larger R_* for higher $\dot{M}_\ast$ remains valid even without the deuterium burning. For simplicity, we do not include the effect of deuterium burning in the following argument. The origin of this trend is higher specific entropy for higher accretion rate. Recall that the stellar radius is larger with the high entropy in the stellar interior as shown by Equation (12). The entropy distributions at M_* = 1 M_☉ for the no-D runs are shown in Figure 15, which indeed demonstrates that the entropy is higher for the higher $\dot{M}_\ast$ cases. The interior entropy is set in the postshock settling layer where radiative cooling can reduce the entropy from the initial postshock value as observed in spiky structures near the surface (e.g., the case of $\dot{M}_\ast \lesssim 10^{-5}\, M_{\odot }\, {\rm yr^{-1}}$ in Figure 15). For efficient radiative cooling, the local cooling time t_cool,s in the settling layer must be shorter than the local accretion time t_acc,s, as explained in Sections 3.1 and 3.2. Figure 16 shows that the ratio t_cool,s/t_acc,s is smaller for lower $\dot{M}_\ast$ and becomes less than unity for ≲10⁻⁵ M_☉ yr^-1. Thus the accreted gas cools radiatively in the settling layer in low $\dot{M}_\ast$ (≲10⁻⁵ M_☉ yr^-1) cases, while in the higher $\dot{M}_\ast$ cases, the gas is embedded into the stellar interior without losing the postshock entropy. The difference in entropy among high $\dot{M}_\ast$ (≳10⁻⁴ M_☉ yr^-1) cases comes from the higher postshock entropy for higher $\dot{M}_\ast$. In those cases, the protostar is enshrouded with the radiative precursor, whose optical depth is larger for the higher accretion rate (see Figures 6 and 13, bottom). Consequently, more heat is trapped in the precursor, which leads to the higher postshock entropy.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 14 shows not only the swelling of the stellar radius occurs even in the "no-D" runs, but also that their timings are almost the same as in the cases with D. This supports our view that the cause of the swelling is outward entropy transport by the luminosity wave rather than the deuterium shell-burning (Section 3.1 and 3.2). Also, in all cases, the epoch of the swelling obeys Equation (16), which is derived by equating the accretion time t_acc and the KH time t_KH,lmax.

Another clear tendency is that the deuterium burning begins at higher stellar mass for higher accretion rate. Similarly, the onset of hydrogen burning and subsequent arrival at the ZAMS are postponed until higher protostellar mass for the higher accretion rate. For example, the ZAMS arrival is at M_* ≃ 4 M_☉ (40 M_☉) for 10⁻⁶ (10⁻³, respectively) M_☉ yr⁻¹. These delayed nuclear ignitions reflect the lower temperature within the protostar for the higher $\dot{M}_\ast$ at the same M_* (Figure 14, lower panel). Since R_* is larger with higher $\dot{M}_\ast$ at a given M_*, typical temperature is lower within the star. Substituting numerical values in Equation (14), we obtain

$$\begin{equation} T_{\rm max} \sim 10^7\,{\rm K} \left(\frac{M_{\ast }}{M_{\odot }} \right) \left(\frac{R_{\ast }}{R_{\odot }} \right)^{-1}, \end{equation} \tag{ 20 }$$

which is a good approximation for our numerical results.

Finally, we consider the significance of deuterium burning for protostellar evolution. Comparing with "no-D" runs, we find that the effect of deuterium burning not only appears later, but also becomes weaker with increasing accretion rates. At the low accretion rates $\dot{M}_\ast \le 10^{-5}\, M_{\odot }\, {\rm yr^{-1}}$, convection spreads into a wide portion of the star owing to the vigorous deuterium burning. In addition, the thermostat effect works and the stellar radius increases linearly during this period. At the higher accretion rates $\dot{M}_\ast \ge 10^{-4}\, M_{\odot }\, {\rm yr^{-1}}$, on the other hand, deuterium burning affects the protostellar evolution only slightly. Convective regions do not appear even after the beginning of deuterium burning for $\dot{M}_\ast = 10^{-3}\, M_{\odot }\, {\rm yr^{-1}}$ (Section 3.1).

A key quantity for understanding this variation is the maximum radiative luminosity at the ignition of deuterium burning. The maximum radiative luminosity within the star L_max at deuterium ignition can be evaluated by eliminating R_* and M_* in Equation (15) with (13), (20), and setting T_max ∼ 10⁶ K:

$$\begin{equation} L_{\rm max} \sim 10^3\,L_{\odot } \left(\frac{\dot{M}_\ast }{10^{-3}\, M_{\odot }\,{\rm yr}^{-1}} \right)^{2.8}. \end{equation} \tag{ 21 }$$

On the other hand, the luminosity from deuterium burning is (Equation 17) $L_{\rm D,st} \sim 10^3\,L_{\odot } (\dot{M}_\ast /10^{-3}\, M_{\odot }\,{\rm yr}^{-1})$. For an accretion rate $\dot{M}_\ast \ge 10^{-3}\, M_{\odot }\, {\rm yr^{-1}}$, L_max > L_max,D: radiation is able to carry away all the entropy generated by the deuterium burning. Therefore, the protostellar interior remains radiative in the high $\dot{M}_\ast$ cases (e.g., run MD3, see Section 3.1), while convection is excited for the lower $\dot{M}_\ast$ (e.g., run MD5, see Section 3.2).

Besides the massive star formation in the local universe, high accretion rates have been invoked in the primordial star formation in the early universe. Without important metal coolants, the temperature remains at several hundred K in the primordial gas (e.g., Saslaw & Zipoy 1967; Matsuda et al. 1969; Palla et al. 1983). The accretion rate, which is given by

$$\begin{equation} \dot{M}_\ast \sim \frac{c_{\rm s}^3}{G} \simeq 10^{-3}\,M_{\odot }/{\rm yr} \left(\frac{T}{500\,{\rm K}} \right)^{3/2}, \end{equation} \tag{ 22 }$$

where c_s and T are the sound speed and temperature, respectively, in star-forming clouds, becomes as high as 10⁻³ M_☉ yr⁻¹ even without turbulence (SPS86). As a result of similar accretion rates, protostellar evolution of the primordial stars well resembles with those of the high $\dot{M}_\ast$ cases studied in this paper (SPS86; Omukai & Palla 2001, 2003). The mass–radius relations for the primordial (Z = 0) protostars are presented in Figure 17, along with their solar-metallicity (Z = 0.02) counterparts for comparison. This indicates that basic evolutionary features are similar despite differences in metallicity: as in the present day protostars, primordial protostars pass through the same four characteristic phases (see Section 3.1 and 3.2). The following two differences, however, still exist. First, the swelling of the protostar and subsequent transition to the KH contraction occur earlier in the Z = 0 cases. For example, for $\dot{M}_\ast =10^{-3}\,M_{\odot }\,{\rm yr}^{-1}$, the maximum radius is attained and the KH contraction begins at 10 M_☉ in the Z = 0.02 case, while it is 7 M_☉ in the Z = 0 case. Recall that the transition to those phases is caused by the decrease in opacity with increasing M_* (Section 3.1). Owing to the lack of heavy elements, the opacity is lower in the primordial gas at the same density and temperature. Efficient radiative heat transport begins earlier, which results in the earlier transitions in the evolutionary phases. Second, the stellar radius in the main-sequence phase is smaller as a result of higher interior temperature (≃10⁸ K) for the primordial stars. Due to the lack of C and N, sufficient energy production by the CN cycle is not available at a temperature similar to that of the solar metallicity case (≃10⁷ K). The KH contraction is not halted until the temperature reaches ≃10⁸ K, where a small amount of carbon produced by the He burning enables the operation of the CN cycle (Ezer & Cameron 1971; Omukai & Palla 2003). Since more contraction is needed, the arrival at the ZAMS can be later in the Z = 0 case, despite the earlier beginning of the KH contraction: for example, it is 30 M_☉ (50 M_☉) in the Z_☉ (Z = 0, respectively) case.

Zoom In Zoom Out Reset image size

The higher the mass accretion rate, the more massive the protostar grows before the arrival at the ZAMS phase (Section 3.1 and 3.3). However, this trend does not continue unlimitedly: for a rate exceeding a threshold value, stars have an upper mass limit set by the radiation pressure onto the radiative precursor.

The evolution of the protostellar radius is shown in Figure 18 (upper panel) for cases with very high accretion rates >10⁻³ M_☉ yr^-1. Up to a certain moment in the KH contraction phase, the evolution in all cases proceeds in a qualitatively similar way although the radius is larger and the onset of swelling and the KH contraction are later for higher $\dot{M}_\ast$ cases. However, for accretion rates $\dot{M}_\ast > 3 \times 10^{-3}\, M_{\odot }\, {\rm yr^{-1}}$, the contraction is halted at ≃40M_☉, followed by abrupt expansion at ≃50 M_☉. In such cases, further steady accretion onto the star is not possible owing to strong radiative pressure exerted onto the radiative precursor. This phenomenon has been found for the primordial star formation by Omukai & Palla (2001, 2003). They found that abrupt expansion occurs when the total luminosity from a protostar, $L_{\rm tot} = L_{\ast } + L_{\rm acc} =L_{\ast } + G M_{\ast } \dot{M}_{\ast }/R_{\ast }$, becomes close to the Eddington luminosity L_Edd. The strong radiation pressure decelerates the accretion flow and then reduces the ram pressure onto the stellar surface. With reduced inward ram pressure, an outward thrust by the interior radiation blows a thin surface layer away.

Zoom In Zoom Out Reset image size

In the KH contraction phase, the interior luminosity L_* exceeds the accretion luminosity L_acc by a factor of a few (see Section 3). However, the interior luminosity L_* alone never reaches the Eddington luminosity. Additional contribution by the accretion luminosity is crucial to boost the total luminosity to the Eddington limit. We define the critical accretion rate $\dot{M}_{\rm \ast, cr}$, above which the abrupt stellar expansion occurs in the KH contraction phase. By using the condition that L_tot reaches L_Edd before the protostar reaches the ZAMS, (Omukai & Palla 2003) derived the critical rate as follows:

$$\begin{equation} L_{\rm ZAMS} + \frac{G M_{\rm ZAMS} \dot{M}_{\rm \ast, cr}}{R_{\rm ZAMS}} = L_{\rm Edd}, \end{equation} \tag{ 23 }$$

that is,

$$\begin{equation} \dot{M}_{\rm \ast, cr} = \frac{4 \pi c R_{\rm ZAMS}}{\kappa _{\rm esc}} \left(1 - \frac{L_{\rm ZAMS}}{L_{\rm Edd}} \right), \end{equation} \tag{ 24 }$$

where κ_esc is the electron-scattering opacity, and the suffix "ZAMS" means quantities of the ZAMS stars. Although the critical rate $\dot{M}_{\rm \ast, cr}$ by Equation (24) depends on the stellar mass, the dependence is very weak and $\dot{M}_{\rm \ast, cr} \simeq 4 \times 10^{-3}\, M_{\odot }\, {\rm yr^{-1}}$ for the primordial stars. Since the ZAMS radius increases with metallicity (see Section 3.4), $\dot{M}_{\rm \ast, cr}$ is expected to increase with metallicity. For example, for Z = 10⁻⁴, $\dot{M}_{\rm \ast, cr}=9 \times 10^{-3}\, M_{\odot }\, {\rm yr^{-1}}$ by Equation (24), and the evolutionary calculation demonstrated that abrupt stellar expansion occurs in fact for higher accretion rates (Omukai & Palla 2003). Our calculations, however, indicate that the critical accretion rate in the solar metallicity case is ≃3 × 10⁻³ M_☉ yr⁻¹, slightly smaller than that of primordial protostars. This apparent contradiction comes from the fact that, in the solar metallicity case, opacity in the accreting envelope is higher than that of the electron scattering one κ_es postulated in the above derivation for the critical accretion rate. In fact, the violent stellar expansion is found to take place with a total luminosity of about a half of the Eddington luminosity, as shown in the lower panel of Figure 18. A marginal case is Run MD3x3, where the contraction is about to be reversed, but the star barely reaches the ZAMS after some loitering. In this case, the Eddington ratio f_Edd = L_tot/L_Edd slightly exceeds 0.5 at M_* ≃ 45 M_☉. For comparison, the case of the primordial star is shown for the same accretion rate, $\dot{M}_\ast = 3 \times 10^{-3}\, M_{\odot }\, {\rm yr^{-1}}$, which is slightly lower than the critical accretion rate (run MD3x3-z0). Because of the earlier transition to the KH contraction phase, the Eddington ratio also begins to increase earlier than the present-day counterpart (run MD3x3). Despite high luminosity almost reaching the Eddington limit, the KH contraction is not reversed and the protostar can grow further. Taking into account the higher opacity in the envelope, we add the Eddington factor of 0.5 to Equation (23):

$$\begin{equation} L_{\rm ZAMS} + \frac{G M_{\rm ZAMS} \dot{M}_{\rm \ast, cr}}{R_{\rm ZAMS}} \simeq 0.5 L_{\rm Edd} \end{equation} \tag{ 25 }$$

for present day protostars. The critical accretion rate becomes $\dot{M}_{\rm \ast, cr} \simeq 3 \times 10^{-3}\, M_{\odot }\, {\rm yr^{-1}}$, which agrees well with our numerical results. Although this critical rate is a function of the stellar mass, its dependence is weak as in the primordial case.

Our results suggest an upper mass limit of pre-main-sequence (PMS) stars at ≃60 M_☉, which corresponds to the mass where the protostar reaches the ZAMS at the critical accretion rate of M_cr ∼ 3 × 10⁻³ M_☉ yr⁻¹. More massive PMS stars cannot be formed by steady mass accretion. With higher accretion rates, the steady accretion is halted by the abrupt expansion at M_* ≃ 50 M_☉. What happens after the abrupt expansion is speculated upon in the next section.

A.1.1. Protostar

For protostellar structure, the following four stellar structure equations are solved:

$$\begin{equation} \left(\frac{\partial r}{\partial M} \right)_t = \frac{1}{4 \pi \rho r^2}, \end{equation} \tag{ A1 }$$

$$\begin{equation} \left(\frac{\partial P}{\partial M} \right)_t = - \frac{GM}{4 \pi r^4}, \end{equation} \tag{ A2 }$$

$$\begin{equation} \left(\frac{\partial L}{\partial M} \right)_t = \epsilon - T \left(\frac{\partial s}{\partial t} \right)_M, \end{equation} \tag{ A3 }$$

$$\begin{equation} \left(\frac{\partial s}{\partial M} \right)_t = \frac{G M}{4 \pi r^4} \left(\frac{\partial s}{\partial p} \right)_T \left(\frac{L}{L_s} - 1 \right) C, \end{equation} \tag{ A4 }$$

with a spatial variable of Lagrangian mass coordinate M. In the above, is the energy production rate by nuclear fusion, s is the entropy per unit mass, L_s is the radiative luminosity with adiabatic temperature gradient, and other quantities have ordinary meanings. The coefficient, C in Equation (A4) is unity if L < L_s (i.e., in radiative layers), and given by the mixing-length theory if L > L_s (i.e., in convective layers). The momentum Equation (A2) simply states the hydrostatic equilibrium. This is valid because the stellar mechanical timescale is generally much shorter than the accretion timescale. In energy Equation (A3), on the other hand, the nonequilibrium term −T(∂s/∂t)_M should be included. This is because the KH timescale can be longer than or comparable to the accretion timescale in our calculations. Thus, the thermal equilibrium is not generally satisfied in accreting protostars. Nuclear reactions included in are deuterium burning and hydrogen burning via the pp-chain and CN-cycle. The elements are homogeneously distributed in convective layers by the convective mixing.

There are two alternative schemes to integrate Equations (A3) and (A4), which are the explicit and implicit schemes:

Explicit Scheme. In the explicit scheme, Equation (A3) is considered as a time development equation of the entropy. In order to update the entropy at each spatial grid from a time step t − Δt to t, we use (∂L/∂M)_t evaluated at the previous time step t − Δt and calculate (∂s/∂t)_M by Equation (A3). Equation (A4) is transformed as,

$$\begin{equation} L = L_s \left[ 1 + \frac{4 \pi r^4}{G M} \left(\frac{\partial P}{\partial s} \right)_T \left(\frac{\partial s}{\partial M} \right)_t \right], \end{equation} \tag{ A5 }$$

with C = 1 (radiative). We substitute the updated entropy gradient (∂s/∂M)_t to Equation (A5) and calculate the luminosity profile L(M,t). This scheme enables the precise integration where the heat transport among mass cells is negligible, i.e., (∂L/∂M)_t is much smaller than other terms in Equation (A3). For example, this quasiadiabatic situation is created in the very opaque interior of the protostar.

Implicit Scheme. In the implicit scheme, on the other hand, the time-derivative term −T(∂s/∂t)_M in Equation (A3) is regarded as a source term to calculate the luminosity profile by spatial integration. We simultaneously integrate Equations (A3) and (A4), evaluating (∂s/∂t)_M at every increment of spatial grids. This scheme works effectively where the heat transport among mass cells works to redistribute entropy. For example, this is the case near the stellar surface. We devote special attention to discretizing the time-derivative (∂s/∂t)_M in Equation (A3). For example, consider the gas element $\dot{M}_\ast \Delta t$ measured from the surface in a stellar model. This element accretes to the star from a time step t − Δt to t, and is still in the accreting envelope at the previous time step t − Δt. The accreting matter plunges to the star through the accretion shock front, which is treated as a mathematical discontinuity in our formulation. Describing the derivative across the discontinuity is a problem. We avoid this difficulty adopting a coordinate conversion (SST80b),

$$\begin{equation} \left(\frac{\partial s}{\partial t} \right)_M = \left(\frac{\partial s}{\partial t} \right)_m + \dot{M}_\ast \left(\frac{\partial s}{\partial m} \right)_t, \end{equation} \tag{ A6 }$$

where m is the inverse mass coordinate m ≡ M_* − M. The term (∂s/∂t)_m can be obtained only with the entropy distribution within the shock front. This treatment is fairly valid, because the structure of the settling layer hardly changes over a few time steps, observing from the shock front.

In our calculations, we properly adopt different schemes in different evolutionary stages. In the early phase of the evolution, we use both the explicit and implicit schemes to construct one stellar model (e.g., SST80b); the explicit scheme is used in the deep interior and is switched to the implicit scheme near the stellar surface. First, we solve the interior with the explicit scheme; we integrate Equations (A1) and (A2) outward with a guess of central pressure P_c. When the explicitly solved region is radiative, we calculate L(M,t) with Equation (A3). Once active deuterium burning begins, however, the explicit update of entropy sometimes makes negative entropy gradients (∂s/∂M)_t < 0. Such a layer should be regarded as a convective layer, and Equation (A5) is not valid any more. We temporarily change the integration method in such a layer. We correct the entropy distribution as (∂s/∂M)_t = 0 there after the explicit update, following the "equal-area rule" developed by SST80b. Once the entropy profile is corrected, we calculate (∂s/∂t)_M and obtain the luminosity distribution by integrating Equation (A3). As the integration approaches the surface, the term (∂L/∂M)_t in Equation (A3) gradually grows. We switch to the implicit method where (∂L/∂M)_t attains $0.5 \,\dot{M}_\ast T (\partial s/\partial m)_t$. We integrate all Equations (A1)–(A4) with a guess of the entropy at the switching point s_sw. When the integration reaches the stellar surface, we check the differences from the required boundary conditions (see Section A.1.2 below). We improve the guessed values of P_c and s_sw, and iterate this procedure until the model satisfies the outer boundary conditions.

The above method successfully works in the early phase of evolution. However, only in the later stage of evolution, where the entropy redistribution among mass cells works even in the deep stellar interior, we need to use another method with the implicit scheme. This occurs when a widespread convective layer appears during active deuterium burning, or when radiative heat transport becomes efficient due to the opacity decrease in the deep interior (see Section 3 for details). In the fully implicit method, we integrate the Equations (A1)–(A4) outward from the center with a guess of the central pressure P_c and entropy s_c. If the integration reaches the stellar surface, we advance to the same fitting process at the surface as in the above method. If the integration diverges and fails before reaching the surface, we adopt another shooting method to a halfway fitting point. We additionally guess the radius R_* and luminosity L_* at the stellar surface, and integrate backward to the fitting point. We adjust boundary values to match the physical quantities at the fitting point, which are obtained by the forward and backward integrations.

The outer boundary condition of the protostar is given by constructing the structure of an accreting gas envelope. The boundary layer connecting the accreting envelope and protostar actually has detailed structure. The acreted gas initially hits a standing shock front and is heated up by compression there. The shocked gas next enters the relaxation layer behind the shock front and is cooled down by radiative loss. In our formulation, this detailed structure is not solved, but treated as a mathematical discontinuity (SST80). We define the following two points across the discontinuity: a postshock point, where gas and radiation temperature first become equal behind the relaxation layer, and a postshock point ahead of the shock front. In our calculations, the former is the outer boundary of the protostar, and the latter is the inner boundary of the accreting envelope. Physical states at each point are related by a radiative shock jump condition.

Here, we focus on cases where the postshock quantities are given in advance by the trial outward integration of the protostar (see Section A.1.1). When the outward integration fails and the shooting to a halfway fitting point is needed, we only slightly modify the same procedure (e.g., see PS91). First, we judge if the accreting envelope is optically thin or thick. Temporarily assuming that the accretion flow is optically thin and in the free fall, velocity and mass density at the preshock point are,

$$\begin{equation} u_{\rm pre} = - \sqrt{\frac{2 G M_{\ast }}{R_{\ast }} }, \end{equation} \tag{ A7 }$$

$$\begin{equation} \rho _{\rm pre} = \frac{\dot{M}_\ast }{4 \pi R_{\ast }^2 u_{\rm pre}}. \end{equation} \tag{ A8 }$$

Ahead of the shock front, outward radiation comes from both the stellar interior and relaxation layer behind the shock front. Therefore, radiation temperature at the preshock point is written as,

$$\begin{equation} T_{\rm rad, pre} = \left(\frac{L_{\ast } + L_{\rm acc}}{4 \pi R_{\ast }^2 \sigma } \right)^{1/4}, \end{equation} \tag{ A9 }$$

where σ is the Stefan–Boltzmann constant, and L_acc is the accretion luminosity defined by Equation (7) Using Equations (A7)–(A9), optical thickness of the accreting envelope is estimated as,

$$\begin{equation} \tau _{\rm pre} = R_{\ast } \rho _{\rm pre} \kappa (\rho _{\rm pre}, T_{\rm pre}), \end{equation} \tag{ A10 }$$

where κ is the Rosseland mean opacity. If τ_pre < 1, the accreting envelope is optically thin and initial adoption of the thin envelope is verified. The jump condition at the accretion shock front is provided by considering the flow structure across the shock front (SST80b). In the optically thin case, the postshock temperature is related to the stellar radius and luminosity (PS91),

$$\begin{equation} T_{\rm post} = \left[ \frac{1}{4 \pi R_{\ast }^2 \sigma } \left(\frac{1}{2} L_{\ast } + \frac{3}{4} L_{\rm acc} \right) \right]^{1/4}. \end{equation} \tag{ A11 }$$

This is the outer boundary condition of the protostar. The protostellar models are constructed for the postshock quantities to satisfy Equation (A11).

If τ_pre > 1, on the other hand, the accreting envelope should actually be optically thick. In this case, we solve the structure of the flow inside the photosphere, i.e., radiative precursor. First, we determine the locus of the photosphere. Assuming the free-fall flow outside the photosphere, optical depth to the photosphere τ_ph can be written following Equations (A7)–(A10) by substituting physical quantities at the photosphere. The only unknown quantity is the luminosity at the photosphere L_ph, which is needed to calculate temperature in Equation (A9). We calculate this photospheric luminosity using the conservation law of the net energy outflow rate L_e through the precursor,

$$\begin{equation} L_e = L_{\rm ph} - \dot{M}_\ast \left(h_{\rm ph} + \frac{1}{2} u_{\rm ph}^2 - \frac{G M_{\ast }}{R_{\rm ph}} \right), \end{equation} \tag{ A12 }$$

where h_ph = h(ρ_ph, T_ph) is the enthalpy of both gas and radiation at the photosphere. Using the known postshock quantities, we can give the numerical value of L_e at the bottom of the accretion flow,

$$\begin{equation} L_e = L_{\ast } - \dot{M}_\ast \left(h_{\rm post} + \frac{1}{2} u_{\rm pre}^2 - \frac{G M_{\ast }}{R_{\ast }} \right). \end{equation} \tag{ A13 }$$

Assigning this value to L_e in Equation (A12), the photospheric luminosity is written with the infall velocity and thermal state at r = R_ph. Consequently, T_ph and τ_ph are given as implicit functions of R_ph. Therefore, we can fix R_ph in an iterative fashion by requiring τ_ph = 1. Once the photospheric radius is fixed, we solve the structure of the precursor by integrating the following Equations inward from the photosphere,

$$\begin{equation} \rho = - \frac{\dot{M}_\ast }{4 \pi r^2 u}, \end{equation} \tag{ A14 }$$

$$\begin{equation} u \frac{\partial u}{\partial r} = - \frac{G M_{\ast }}{r^2} - \frac{1}{\rho } \frac{\partial P}{\partial r} + \frac{\kappa }{c} F, \end{equation} \tag{ A15 }$$

$$\begin{equation} \frac{\partial T}{\partial r} = - \frac{3 \rho \kappa }{16 \sigma T^3} F, \end{equation} \tag{ A16 }$$

where F is the radiative flux. Using Equations (A12) and (A14), F is written as a function of ρ, u, and T,

$$\begin{equation} F = - \rho u \left(\frac{L_e}{\dot{M}_\ast } + h + \frac{1}{2} u^2 - \frac{G M_{\ast }}{r} \right). \end{equation} \tag{ A17 }$$

When the integration reaches the preshock point, a jump condition bridges the preshock and postshock points. An isothermal jump condition is applied for the optically thick radiative shock front,

$$\begin{equation} T_{\rm post} = T_{\rm pre}. \end{equation} \tag{ A18 }$$

This is another outer boundary condition in addition to Equation (A11).

We construct initial models by solving their structure. The method of calculation is based on that of evolutionary calculations, but with slight modifications. First, we assume the entropy distribution in the stellar interior as an increasing function with M,

$$\begin{equation} s(M) = s_{c,0} + \beta \frac{k_{\rm B}}{m_{\rm H}} \frac{M}{M_{*,0}}, \end{equation} \tag{ A19 }$$

where s_c,0 is specific entropy at the stellar center, k_B is Boltzmann constant, m_H is atomic mass unit, and β(>0) is a free parameter. In the core interior, we integrate Equations (A1) and (A2) outward beginning with guessed s_c,0 and central pressure, p_c. According to Equation (A5), the luminosity profile there is calculated using an entropy gradient obtained from (A19). Whereas this procedure is valid in the adiabatic interior, entropy and luminosity profiles should be consistently calculated in place of Equation (A19) near the surface (also see Section A.1.1). The switching point is defined as follows. Using Equations (A3) and (A6), the energy equation can be written as,

$$\begin{equation} \left(\frac{\partial L}{\partial M} \right)_t = \dot{M}_\ast T \left(\frac{\partial s}{\partial M} \right)_t, \end{equation} \tag{ A20 }$$

omitting the nuclear energy production and time-derivative terms. Once the luminosity gradient (∂L/∂M)_t becomes comparable to $\dot{M}_\ast T (\partial s/\partial M)_t$ calculated with Equation (A19), we solve the full equations including Equation (A20). The guessed s_c,0 and p_c are adjusted so that the core satisfies outer boundary conditions, as explained in Section A.1.2. We have confirmed that some variations of initial models only affect subsequent protostellar evolution in a very early phase. In this paper, we adopt β = 1 as the entropy slope in Equation (A19).

For the opacity, we adopt OPAL tables (e.g., Igresias & Rogers 1996) with the composition given by Grevesse & Noels (1993) for temperature T > 7000 K. For T < 7000 K, we use other opacity tables based on calculations by Alexander & Ferguson (1994). Contribution from grains is excluded in the adopted composition (Asplund et al. 2005; Cunha et al. 2006). In all runs, we adopt X = 0.7, Y = 1 − X − Z, and relative abundances of heavy elements following Grevesse & Noels (1993) composition.

The effect of different deuterium abundances on the stellar interior structure in early phases is shown in Figure 22 for cases with [D/H] = 3 × 10⁻⁵, 2.5 × 10⁻⁵, and 1 × 10⁻⁵. The middle panel is the same as Figure 7 (upper), but reproduced for comparison. As described in Section 3.2, our run MD5 exhibits somewhat complex evolution of convective layers in the early phase. A convective layer first appears at M_* ≃ 0.3 M_☉, just after the deuterium burning begins. This convective layer disappears at M_* ≃ 0.6 M_☉. The second convective layer emerges just outside the former outer boundary of the first one at M_* ≃ 0.7 M_☉, and quickly reaches the surface. SST80a calculated the same model, but did not observe such a complex evolution. In their calculation, the first convective layer survives and eventually incorporates most of the stellar interior.

Zoom In Zoom Out Reset image size

This difference comes from the fact that the evolution of convective layers sensitively depends on the initial deuterium abundance. For example, in the case with the deuterium abundance [D/H] = 3 × 10⁻⁵ (Figure 22, upper; run MD5-dh3), slightly higher than the fiducial value, the convective layer first appears at the same epoch as the fiducial case with [D/H] = 2.5 × 10⁻⁵, and extends outward. In this case, another convective layer soon emerges at the surface and begins to extend inward. These two convective layers finally merge at M_* ≃ 0.6 M_☉, and render most of the interior convective. At M_* ≃ 1 M_☉, the convective region contains about 97% of the total mass. These features are very similar to those found by SST80a, who terminated their calculation at M_* = 1 M_☉. Our calculation shows that the subsequent evolution in this case is almost the same as in our fiducial run MD5. The central radiative core gradually expands as the stellar mass increases, and finally occupies most of the interior at M_* ≃ 4 M_☉.

This sensitive dependence on the deuterium abundance is explained as follows. The expansion of convective layers is caused by increase in the specific entropy within it, which is generated by the nuclear burning at the bottom of the layer. For continuing expansion of the convective layer, deuterium needs to be supplied continuously. Fresh deuterium becomes available when the convective layer newly incorporates the radiative region. This deuterium immediately spreads over the convective layer, and is consumed for further nuclear fusion. For sufficiently high deuterium abundance, acquired deuterium is high enough to maintain the active nuclear burning. In this case, the convective layer continuously grows as in the case with [D/H] = 3 × 10⁻⁵ (run MD5-dh3). On the other hand, for low deuterium abundances, the nuclear burning becomes too weak to maintain development of the convective layer. The convective layer ceases to grow and disappears as in the case with [D/H] = 2.5 × 10⁻⁵ (run MD5). For even smaller deuterium abundance [D/H] = 1 × 10⁻⁵ (Figure 22, lower; run MD5-dh1), no convective layer appears even after ignition of the deuterium burning until at M_* ∼ 0.8 M_☉, when the surface convective layer finally emerges. In this case, the subsequent mass–radius relation is almost same as in the case without deuterium burning (run MD5-noD) although the surface convective layer pesists until M_* ∼ 3 M_☉. Our fiducial run MD5 is just intermediate between the deuterium-rich (run MD5-dh3) and -poor (MD5-dh1) cases presented in Figure 22. Therefore, a slight difference in input physics affects the internal structure significantly.

Evolution of accreting intermediate-mass (1 M_☉ < M_* < 8 M_☉) protostars has been studied in some previous work (e.g., PS91). Our results agree well with the previous work, but still show some differences even with the same accretion rate, which can be ascribed to different initial models. Whereas our initial model is a tiny radiative protostar with ≪1 M_☉ (see Section A.2), the initial model of PS91, for example, is a 1 M_☉ fully convective protostar. We confirmed that our code reproduces almost the same evolution as that of PS91, if we begin the calculation with the same initial model.

The initial model of PS91 is constructed in a manner similar to that described in Section A.2. The homogeneous entropy distribution is adopted in place of Equation (A19), and its value is taken from the semianalytic model of Stahler (1988). We have reconstructed the model of Stahler (1988), and obtained s_c,0 = −4.12k_B/m_H for a 1 M_☉ star with our adopted opacity tables. As in Section A.2, we separately solve the core interior and the surface superadiabatic layer. First, we integrate the core interior outward from the center with a guessed value of the central pressure. We also guess the luminosity at the bottom of the surface superadiabatic layer in advance, and record the entropy gradient ∂s/∂M using Equation (A4). Once the recorded ∂s/∂M attains a small finite value, we switch the scheme and integrate all Equations (A1)–(A4). The guessed quantities are adjusted for the core to satisfy the outer boundary conditions. The constructed model is a fully convective star, whose mass and radius are 1 M_☉ and 4.2 R_☉.

We calculate the subsequent evolution beginning with this initial model (run MD5-ps91). Figure 23 shows the evolution of such a protostar. Our calculation reproduces the basic features of a calculation done by PS91, which adopted the same accretion rate and accretion-shock outer boundary. In an early phase, the star remains fully convective and deuterium burning occurs around the stellar center. Deuterium concentration significantly falls in this phase. This is because the central temperature continuously increases, and total burning rate L_D slightly exceeds the steady burning rate L_D,st (e.g., Palla & Stahler 1993). These features are somewhat different from our fiducial run MD5, where a radiative core persists throughout the evolution and gradually expands with the increase in M_* by accretion (see Figure 7). These differences come from different initial models. Figure 7 also shows that mass-averaged deuterium concentration f_d,av does not fall much below 0.01 in our run MD5. This is because the deuterium burning layer gradually moves to the outer layer, where temperature is always around 10⁶ K and L_D ∼ L_D,st. At M_* = 2 M_☉, for example, f_d,av ∼ 0.04 in run MD5, but only f_d,av ∼ 10⁻⁵ in run MD5-ps91. Therefore, our calculation predicts that intermediate-mass protostars have a larger amount of deuterium than assumed by PS91. Figure 23 shows that the fully convective phase ends at M_* ≃ 2.6 M_☉, when a thin radiative layer, i.e., radiative barrier, suddenly appears within the star. The radiative barrier blocks inward convective transport of the accreted deuterium. Consequently, the inner region, where deuterium has run out, immediately returns to being radiative. After that, the convective layer remains only above the radiative core, and deuterium burning occurs at the bottom of the convective layer. As in our fiducial run MD5, the protostar swells up at M_* ≃ 3.4 M_☉. PS91 have concluded that this swelling is caused by the shell-burning of deuterium. To see the role of the deuterium shell-burning in the swelling, we also calculate evolution with cessation of the deuterium burning after the appearance of the radiative barrier. We find that the swelling occurs similarly, although it is slightly delayed, even without deuterium burning. This is caused by the outward transport of embedded entropy, which is observed as propagation of a luminosity wave (see Sections 3.1 and 3.2). Although PS91 have attributed the swelling only to the shell-burning of deuterium, we conclude that the main driver of the swelling is always the propagation of the luminosity wave. After the turn-around of the radius, the star enters the subsequent KH contraction phase. The calculated mass–radius relation thereafter gradually converges with that in our fiducial run MD5. The epoch of the Hydrogen burning, M_* ≃ 7 M_☉, also agrees with that in the fuducial run.

Zoom In Zoom Out Reset image size