FORMATION OF MASSIVE PRIMORDIAL STARS: INTERMITTENT UV FEEDBACK WITH EPISODIC MASS ACCRETION

We use a modified version of the public multi-dimensional magneto-hydrodynamics code PLUTO 3.0 (Mignone et al. 2007). The original code had been modified to study present-day high-mass star formation by implementing the sink cell philosophy with the embedded SE code outlined above, as well as self-gravity and radiation transfer modules (e.g., Kuiper et al. 2010a, 2010b, 2011; Kuiper & Klessen 2013). For this investigation we added several physics modules necessary to study primordial star formation.

The governing equations of hydrodynamics are

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}+{\rm{\nabla }}\cdot (\rho {\boldsymbol{v}})=0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial (\rho {\boldsymbol{v}})}{\partial t}+{\rm{\nabla }}\cdot (\rho {\boldsymbol{v}}\otimes {\boldsymbol{v}})=-\rho {\rm{\nabla }}{\rm{\Phi }}-{\rm{\nabla }}p,\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial e}{\partial t}+{\rm{\nabla }}\cdot (e{\boldsymbol{v}})=-p{\rm{\nabla }}\cdot {\boldsymbol{v}}+{\rm{\Gamma }}-{\rm{\Lambda }},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&p=(\gamma -1)e,\end{eqnarray} \tag{ 4 }$$

where ρ and p are the gas mass density and pressure, ${\boldsymbol{v}}$ the 3D velocity vector, Φ the gravitational potential, e the gas thermal energy density, Γ and Λ the heating and cooling rates per volume, respectively. We allow the adiabatic exponent γ to vary depending on the chemical abundances and gas temperature (e.g., Omukai & Nishi 1998). We calculate the energy source term ${\rm{\Gamma }}-{\rm{\Lambda }}$ with the same thermal processes as considered in HS11 (Table S1 in HS11, but see Section 2.1.3 below for some differences). Note that Equation (2) does not include the viscosity terms which enabled angular momentum transport under 2D axial symmetry (e.g., α-viscosity, Shakura & Sunyaev 1973). Instead, gravitational torques induced by the non-axisymmetric disk structure (e.g., spiral arms, clumps) operate in 3D.

As in Kuiper et al. (2011), we solve the above equations in spherical coordinates. The spatial grids for the polar (θ) and azimuthal (ϕ) directions are homogeneously distributed over 0 ≤ θ ≤ π and 0 ≤ ϕ ≤ 2π. The angular resolutions are thus Δθ = π/N_θ and Δϕ = 2π/N_ϕ, where N_θ and N_ϕ are the grid numbers. We set Δθ = Δϕ by taking N_ϕ = 2N_θ (Table 1). The radial grid size Δr logarithmically increases with increasing r as

$$\begin{eqnarray}&&{\rm{\Delta }}r(r)=r({10}^{f}-1),\end{eqnarray} \tag{ 5 }$$

where $f\equiv \mathrm{log}({r}_{{\rm{max}}}/{r}_{{\rm{min}}})/{N}_{r}$, N_r is the number of radial grid cells, and r_min (r_max) is the radial inner (outer) boundary of the computational domain.

Table 1.  Cosmological Cases Considered

Cases	NR · Nθ · Nϕ	rmax (104 au)	Mg,tot (${M}_{\odot }$)	Δt (104 year)	M*,3D (${M}_{\odot }$)	M*,2D (${M}_{\odot }$)
A	96 · 32 · 64	150	1.0 × 104	7	462.4a	751.3
B	96 · 32 · 64	150	1.3 × 104	7	598.8a	283.9
B-NFb	96 · 32 · 64	150	⋯	7	⋯	⋯
B-NF-HR2c-m0d	128 · 64 · 128	4	1.1 × 103	0.3	⋯	⋯
B-NF-HR2-m40	128 · 64 · 128	4	⋯	0.3	⋯	⋯
B-NF-HR2-m70	128 · 64 · 128	4	⋯	0.3	⋯	⋯
B-NF-HR2-m120	128 · 64 · 128	4	⋯	0.3	⋯	⋯
B-NF-HR4c-m20	192 · 128 · 256	0.67	⋯	0.1	⋯	⋯
B-NF-HR4-m20-fce	192 · 128 · 256	0.67	⋯	0.1	⋯	⋯
B-NF-HR4-m20-lskf	192 · 128 · 256	0.67	⋯	0.1	⋯	⋯
C	64 · 32 · 64	4	1.0 × 103	12	286.0	160.3
C-HR2-m0	128 · 64 · 128	4	⋯	5.5	⋯	⋯
D	64 · 32 · 64	4	490	8	67.4	55.5
D-NF	64 · 32 · 64	4	⋯	8	⋯	⋯
D-DFb	64 · 32 · 64	4	⋯	8	⋯	⋯
E	96 · 32 · 64	150	2.7 × 103	6	14.3	24.4

Notes. Column 2: cell numbers, Column 3: position of the radial outer boundary, Column 4: total mass contained in the computational domain, Column 5: time duration followed Column 6: final stellar masses in the current work, Column 7: final stellar masses obtained in 2D SE-RHD simulations (Hirano et al. 2014).

^aStellar masses at the end of the simulations, t = 7 × 10⁴ years. ^bSuffix "NF" represent cases with no UV feedback and "DF" with only dissociation (FUV) feedback. ^cSuffix "HR2" represents cases with doubled and "HR4" with quadrupled spatial resolutions. ^dSuffix "mX" represents the points for the restart after doubling the resolution when M_* ≃ X ${M}_{\odot }$. ^eStringent limit for the cooling with f_limit = 24 (see Section 2.1.2). ^fLarger sink cell with r_min ≃ 52 au (r_min = 30 au for all other cases).

Download table as:  ASCIITypeset image

In order to avoid very short timesteps required for the dense gas near and within the star, we do not resolve the central structure. We instead assume a fixed central spherical sink cell of radius r_min, within which a single star grows. The sink cell is assumed to be semi-permeable, meaning accretion flow into the sink is allowed, but there is no outflow. We evaluate mass accretion rates onto the star by measuring the gas inflow rates into the sink (also see Section 2.2). We set the sink radius r_min = 30 au for our fiducial cases. Although the radius of an accreting protostar is generally much smaller than this sink size when ${\dot{M}}_{\ast }\lesssim {10}^{-2}\;{M}_{\odot }\;{{\rm{yr}}}^{-1}$ (e.g., Omukai & Palla 2003), an accreting star expands rapidly when ${\dot{M}}_{\ast }\gtrsim {10}^{-2}\;{M}_{\odot }\;{{\rm{yr}}}^{-1}$, and inflates up to red supergiant radii. Hosokawa et al. (2012a) derive the mass–radius relationship for such a "supergiant protostar" accreting material above the critical rate,

$$\begin{eqnarray}&&{R}_{\ast }({M}_{\ast })\simeq 11.14\;{\rm{au}}{\left(\displaystyle \frac{{M}_{\ast }}{100{M}_{\odot }}\right)}^{1/2},\end{eqnarray} \tag{ 6 }$$

which is comparable to our sink-cell size.

In the following, we will see that accreting protostars enter the supergiant stage for many of the examined cases. Even with the sink size r_min = 30 au, our spherical coordinate system offers relatively high spatial resolution in the vicinity of the sink, (Δr × rΔθ)_min ≃ 3.57 au × 3.12 au for fiducial cases. This is even better than the finest resolution used in HS11 ∼10 au, so that the disk scale height is always resolved with several grid cells even in the innermost regions.

With this sink-cell technique, the gravitational potential Φ has two components: one arising from the central star Φ_* and one from the gas' self-gravity Φ_g,

$$\begin{eqnarray}&&{\rm{\Phi }}={{\rm{\Phi }}}_{\ast }+{{\rm{\Phi }}}_{g}.\end{eqnarray} \tag{ 7 }$$

The specific gravitational force due to the central star is

$$\begin{eqnarray}&&-{\rm{\nabla }}{{\rm{\Phi }}}_{\ast }=-\displaystyle \frac{{{GM}}_{\ast }}{{r}^{2}}{{\boldsymbol{e}}}_{r},\end{eqnarray} \tag{ 8 }$$

where ${{\boldsymbol{e}}}_{r}$ is the radial unit vector. The potential for the gas' self-gravity Φ_g is obtained by solving Poisson's equation

$$\begin{eqnarray}&&{\rm{\Delta }}{{\rm{\Phi }}}_{g}=4\pi G\rho .\end{eqnarray} \tag{ 9 }$$

We make use of the Poisson solver developed in Kuiper et al. (2010a), which has been well tested and used with the PLUTO code.

Truelove et al. (1997) have suggested that, in order to avoid artificial fragmentation in numerical simulations, the Jeans length has to be spatially resolved by at least four cells. Recent studies argue that even higher resolutions are needed to correctly capture the gravity-driven turbulence that generally appears in the star formation process (e.g., Federrath et al. 2011; Turk et al. 2012; Meece et al. 2014), suggesting that the Jeans length should be resolved with 32–64 cells. Unfortunately, it is computationally prohibitive to consider the long-term (∼10⁵ years) evolution while fully satisfying this condition. We adopt the following ansatz instead:

For each cell and each timestep, we calculate the local Jeans length λ_J and monitor the ratio of λ_J to the greater of Δr and rΔθ. When this ratio falls below a certain fixed number f_limit, we artificially squelch the cooling function in the energy equation, i.e.,

$$\begin{eqnarray}&&{\rm{\Lambda }}\to 0\quad {\rm{for}}\quad \displaystyle \frac{{\lambda }_{{\rm{J}}}}{\mathrm{max}({\rm{\Delta }}r,r{\rm{\Delta }}\theta )}\lt {f}_{{\rm{limit}}}.\end{eqnarray} \tag{ 10 }$$

We normally take f_limit = 12 as our fiducial value. In practice, for a smooth transition, we evaluate the ratio $\xi \equiv {f}_{{\rm{limit}}}\mathrm{max}({\rm{\Delta }}r,r{\rm{\Delta }}\theta )/{\lambda }_{{\rm{J}}}$ and multiply the cooling rate by a suppressing factor when ξ exceeds unity,

$$\begin{eqnarray}&&{C}_{{\rm{limit}}}=\mathrm{exp}\left[-{\left(\displaystyle \frac{\xi -1}{0.1}\right)}^{2}\right],\end{eqnarray} \tag{ 11 }$$

which falls below 10⁻⁴ for ξ ≳ 1.3.

With this procedure, the limit of radiative cooling varies across the grid due to differing spatial resolutions. Choosing a finer grid resolution $\mathrm{max}({\rm{\Delta }}r,r{\rm{\Delta }}\theta )$ overall means that the condition in Equation (10) is satisfied for fewer cells and cooling operates without suppression for a larger number of cells. Because the disk becomes more unstable with more efficient cooling, disk fragmentation occurs more readily in simulations at higher grid resolution, which is also seen in other studies (e.g., Machida & Doi 2013). We study the effects of varying the spatial resolution in Section 4.3 and show that, contrary to the predictions of previous studies, more fragmentation does not necessarily reduce the final stellar masses, but rather facilitates the formation of massive stars by inducing intermittent UV feedback more often.

Stellar UV radiation. We consider different kinds of radiation which play different roles in our simulations. The stellar UV flux is the primary source of UV feedback which potentially halts the mass accretion onto stars and fixes the final stellar masses. We solve the transfer of ionizing (EUV) photons emitted by a protostar with a frequency-averaged approximation,

$$\begin{eqnarray}&&\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial ({r}^{2}{F}_{{\rm{EUV}}})}{\partial r}=-n(1-x){\sigma }_{{\rm{EUV}}}{F}_{{\rm{EUV}}},\end{eqnarray} \tag{ 12 }$$

where F_EUV is the photon number flux above the Lyman limit, n the particle density, x the degree of ionization, and σ_EUV the photoionization cross section per particle, which is a function of the stellar effective temperature T_eff. Although the frequency-dependent UV transport can somewhat broaden primordial ionization fronts, which can enhance the thin-shell instabilities in the fronts (Whalen & Norman 2008), we do not expect that such a effect largely modifies the dynamics especially with dense media in the vicinity of the star. Equation (12) is easily integrated

$$\begin{eqnarray}&&{F}_{{\rm{EUV}}}=\displaystyle \frac{{S}_{{\rm{EUV}}}}{4\pi {r}^{2}}\mathrm{exp}(-{\tau }_{{\rm{EUV}}}),\end{eqnarray} \tag{ 13 }$$

where S_EUV is the stellar EUV emissivity, and τ_EUV the optical depth given by

$$\begin{eqnarray}&&{\tau }_{{\rm{EUV}}}={\int }_{{r}_{{\rm{min}}}}^{r}n(1-x){\sigma }_{{\rm{EUV}}}\;{{dr}}^{\prime }\end{eqnarray} \tag{ 14 }$$

(We implicitly assume no UV absorption within the sink cell).

We solve Equations (13) and (14) making use of the adopted spherical grid (Section 2.1.1), i.e., along the radial cells having the same angular coordinate (θ_j, ϕ_k). The obtained EUV flux F_EUV is used to provide the photoionization rate in the chemical rate equations and the associated heating rate in the energy equation. We consider the ionization of hydrogen only, treating helium as "hydrogenic," i.e., singly ionized within an H ii region. For the purposes of our investigation, this approximation is not critical, but certainly could influence the precise numerical values of our results.

Within an H ii region there is a component of diffuse EUV radiation emitted by the recombining gas, namely recombination of hydrogen directly into the ground state and recombinations of helium into the lowest states. HS11 have separately solved the transfer of diffuse hydrogen recombination photons using the flux-limited diffusion (FLD) approximation (Levermore & Pomraning 1981; Richling & Yorke 1997). However, UV feedback on the accreting material is mostly caused by the direct EUV stellar flux (also see Tanaka et al. 2013). Thus, we only solve for the transfer of the direct stellar EUV component and use the "on-the-spot" approximation (each recombination of hydrogen directly into the ground state causes the immediate ionization of a hydrogen atom in the vicinity). In praxis this is done by using a hydrogen recombination coefficient that excludes recombinations directly into the ground state.

As for the FUV radiation transfer, we only consider the stellar direct component, by exactly the same method using the self-shielding function (Draine & Bertoldi 1996) as in HS11. We recognize that by excluding the diffuse components of EUV and FUV, one will encounter "shadowing" effects caused by UV-opaque clumps and the accretion disk itself. Our previous comparison studies with the 2D SE-RHD code have shown that this simplification does not significantly change the accretion history or final masses of the stars formed, which is the primary goal of this investigation.

Molecular hydrogen line emission. We follow the method used in HS11 for calculating the cooling rate via H₂ line emission; we sum up the contributions from all "allowed" quadrupole transitions between the rotational and vibrational levels of J = 0–20 and v = 0–2. To estimate the trapping effect for each transition, we approximate its optical depth τ by the local Jeans length λ_J, i.e., τ = αλ_J, where α is the transition's absorption coefficient. This is a crude approximation compared to recently developed novel (and much more computationally expensive) methods (e.g., Clark et al. 2012; Hirano & Yoshida 2013; Greif 2014; Hartwig et al. 2015). In our current simulations, however, it is not reasonable to model the detailed thermal disk structure, because of our limitations in following the exact thermal state of the densest gas (Section 2.1.2). For a given spatial resolution, our treatment underestimates the cooling and thus leads to less gravitationally unstable disks. It therefore affects the quantitative behavior of disk fragmentation, such as the exact number and sizes of the emerging fragments. To assess the impact of the details of disk fragmentation on our final results, we vary our spatial resolution (see Section 4.3).

Continuum emission. As described in HS11, continuum cooling via H⁻ free-bound emission is important for the gas infalling onto the disk surface. HS11 have implemented this cooling process by solving continuum radiation transfer using the FLD approximation. In their 2D simulations, however, the H⁻ continuum emission is always optically thin throughout the evolution. This knowledge allows us to use a simpler approximation for the current 3D simulations; we simply subtract the H⁻ bounding energy 0.755 eV for each formation reaction H + e $\to \quad {{\rm{H}}}^{-}+\gamma$ (R10 in Table 2) from the gas' internal energy. We have also implemented this method in our 2D SE-RHD code and confirmed the validity of this approximation.

Table 2.  Chemical Reactions Included

No.	Reactions	References
R1	H + e $\to$  H+ + 2 e	1
R2	H+ + e $\to$  H + γ	2
R3a	H− + H $\to$  H2 + e	3
R4	H2 + ${{\rm{H}}}^{+}\;\to$ ${{{\rm{H}}}_{2}}^{+}$ + H	4
R5	H2 + e $\to$  2 H + e	4
R6a	H2 + H $\to$  3 H	5
R7a	3 H $\to$  H2 + H	6
R8	2 H + ${{\rm{H}}}_{2}\;\to$ 2 H2	7
R9	2 H2 $\to$  2 H + H2	7
R10	H + e $\to$  H− + γ	4
R11	H + $\gamma \;\to$ H+ + e	2
R12	H2 + $\gamma \;\to$ 2 H	8, 9
R13	2 H $\to$  H+ + e + H	7
R14	H− + e $\to$  H + 2 e	1
R15	H− + ${{\rm{H}}}^{+}\;\to$ 2 H	4
R16b	H− + ${{\rm{H}}}^{+}\;\to$ ${{\rm{H}}}_{2}^{+}$ + e	4
R17b	H− + $\gamma \;\to$ H + e	10
R18b	H + ${{\rm{H}}}^{+}\;\to$ ${{{\rm{H}}}_{2}}^{+}$ + γ	4
R19b	${{{\rm{H}}}_{2}}^{+}$ + H $\to$ H2 + H+	4
R20b	${{{\rm{H}}}_{2}}^{+}$ + e $\to$  2 H	4
R21b	${{{\rm{H}}}_{2}}^{+}$ + ${{\rm{H}}}^{-}\;\to$ H2 + H	11

Notes.

^aReaction rates updated with respect to HR14. ^bAdditional reactions added with respect to HR14.

References. (1) Abel et al. (1997), (2) Osterbrock & Ferland (2006), (3) Kreckel et al. (2010), (4) Galli & Palla (1998), (5) Martin et al. (1996), (6) Forrey (2013), (7) Palla et al. (1983), (8) Tielens & Hollenbach (1985), (9) Draine & Bertoldi (1996), (10) John (1988), (11) Millar (1991).

Download table as:  ASCIITypeset image

We do not include the cooling via H₂ collision-induced continuum emission (CIE, Omukai & Nishi 1998). This is because CIE only becomes important for the dense molecular gas with n ≳ 10¹³ ${{\rm{cm}}}^{-3}$, which only occasionally occurs in our simulations. Such dense gas normally lies close to the central star, primarily within the sink cell. Only with the enhanced spatial resolution considered in Section 4.3 below, however, does the density attain ∼10¹⁴ ${{\rm{cm}}}^{-3}$ in the densest parts of disks. Note that we do not pursue the exact thermal state of such tiny dense parts anyway, because of our inability to follow the long-term evolution at sufficiently high spatial resolution (see discussion in Section 2.1.2).

We have implemented an updated version of the primordial chemistry network used in HS11 and HR14 to the PLUTO code. Here, we solve for the non-equilibrium chemical abundances of H, H₂, H⁺, e, H⁻, and ${{{\rm{H}}}_{2}}^{+}$ without the approximation of using the equilibrium value for H⁻ (as employed by e.g., Abel et al. 1997, HS11). The reactions included are summarized in Table 2.

As described in Section 3.1, HD molecules will form during the early collapse in some cases. The resulting HD molecular cooling reduces the temperature and density in the accretion envelope, which in turn leads to a somewhat slower initial protostellar growth (e.g., Yoshida et al. 2007, HS12). However, this effect ends during the early collapse by the time n ∼ 10⁸ ${{\rm{cm}}}^{-3}$ and is unimportant thereafter. While we include deuterium chemistry for cosmological simulations following the early collapse stage (see Section 3.1 below), we do not in the current network for the accretion stage. We have checked with our previous 2D SE-RHD code with the same setting, starting after the density exceeds ∼10¹⁰ ${{\rm{cm}}}^{-3}$ as in HS12, that omitting the deuterium chemistry hardly changes the results.

We first describe in detail the evolution of cases D and E, which produced ordinary massive stars with tens of solar masses. As shown below, the evolution for these cases is similar to that found in our previous 2D SE-RHD simulations (HS11; HS12).

Case D. Figure 7 presents the evolution of the temperature structure for case D as an H₂ PDR and an H ii region expand into the accretion envelope. After the protostar enters the KH contraction stage (see Figure 6), the stellar FUV emissivity rises and forms a bipolar H₂ PDR (Figure 7(a)). At the same time, the stellar EUV emissivity increases and slightly later creates a bipolar H ii region within the larger PDR (Figure 7(b)). The PDR and H ii region continue to grow and expand as shown for the epoch t = 5 × 10⁴ years (Figure 7(c)), but these regions are still centrally "pinched" because a circumstellar disk blocks the stellar UV radiation close to the disk plane. However, even the gas protected from the UV photons is disturbed by shocks that propagate into the envelope ahead of the ionization fronts. Comparing Figures 7(b) and (c), we see that even the gas outside the PDR is heated up by shock compression. HS11 have shown that the mass influx from the envelope onto the disk is hindered by this effect. Moreover, the disk loses mass via photoevaporation, being exposed to the stellar EUV radiation. The H ii region is filled by the photoevaporating outflow launched from the disk surface. Ultimately, the circumstellar disk is destroyed due to photoevaporation, allowing the H ii region to expand in equatorial directions (Figure 7(d)), and accretion onto the star is effectively shut off for M_* ≃ 67.4 ${M}_{\odot }$ at t = 8 × 10⁴ years.

Zoom In Zoom Out Reset image size

In Figure 8, we plot the evolution of the stellar mass for case D (solid black line). For comparison, we show the stellar mass growth for case D-NF (blue dashed line), for which all UV feedback effects had been turned off. We see that mass accretion onto the protostar is significantly reduced by UV feedback effects. This figure also displays the stellar mass growth for a test case (case D-DF) that only included photodissociation (FUV) feedback (solid magenta line). Comparing these mass growth histories we conclude that FUV feedback does hinder mass accretion somewhat, but its impact is much weaker than EUV feedback.

Zoom In Zoom Out Reset image size

Figure 9 shows the variation of the envelope temperature–density–velocity structure with differing UV feedback effects. For case D-NF with neither EUV nor FUV feedback (Figure 9(a)), the disk is surrounded by a warm (T ≃ 300–500 K) envelope, where compressional heating due to gas infall is balanced by H₂ molecular line cooling. The gas temperature is higher (T ≃ 5000 K) in a central region r ≲ 10³ au, where H⁻ free-bound emission is the main cooling process. For case D-DF with only FUV feedback (Figure 9(b)), the warm envelope is larger, because FUV radiation reduces the abundances of efficient coolants via H₂ photodissociation and H⁻ photodetachment. However, this only occurs in polar regions. Because the bulk of stellar FUV radiation is blocked by the disk, the envelope structure is hardly affected in equatorial directions, where the disk mostly accretes gas from the envelope. As seen in Figure 9(c) EUV feedback dramatically changes the envelope structure. In addition to the bipolar photoevaporation flow within the H ii region, portions of the envelope outside the H ii region are expelled from the computational grid as shocks propagate through the envelope and spread the pressure excess of ionized gas.

Zoom In Zoom Out Reset image size

Our results demonstrate that FUV feedback does reduce mass accretion rates, which qualitatively agrees with 3D RHD simulations performed by Susa (2013). However, the FUV feedback seen in our simulations seems to be weaker than in Susa (2013), who claims that FUV feedback is largely responsible for terminating mass accretion onto primordial protostars. We discuss possible causes of this difference in Section 4.2.3 below.

Case E. The evolution of case E is qualitatively similar to that for case D, except that HD molecular line emission also contributes to the cooling in the early collapse stage. During collapse, in general, the thermal evolution of the central homogeneous core strongly affects the radial structure of the accretion envelope. Since the radius of the core is well approximated by the Jeans length, the envelope temperature and density at a given radius r obey r = λ_J ∝ T^1/2ρ^−1/2, i.e., a lower temperature implies a lower density. In case E the accretion envelope thus has relatively low temperatures and densities, which can been seen in Figure 4. KH contraction begins at lower stellar masses with correspondingly lower accretion rates (Figure 6). The final stellar mass is also lower than for case D, M_* ≃ 14.3 ${M}_{\odot }$.

Next we describe the results for cases A–C, for which very massive (M_* ≳ 250 ${M}_{\odot }$) stars form that pass through a supergiant evolutionary stage.

Case B. Figures 10 and 11 show that for case B, variable mass accretion leads to abrupt and radical changes of the stellar radius and the associated recurrent extinction and re-formation of H ii regions. A bipolar H ii region, which have formed prior to t ≃ 4.5 × 10⁴ years (Figure 10(a)), begins to retreat in a direction away from the star a few hundred years later (Figures 10(b), (c)), shortly after the accretion rates rise above ∼10⁻² ${M}_{\odot }\;{{\rm{yr}}}^{-1}$ and the protostar rapidly inflated to a supergiant (Figure 11(a)). Due to the protostar's expansion the stellar EUV emissivity accordingly drops by almost ten orders of magnitude (Figure 11(b)), initiating the extinction of the H ii region.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 12 depicts the evolution of the stellar interior structure for this case B. Although the star initially begins KH contraction when M_* ≃ 30 ${M}_{\odot }$, the contraction is interrupted by a phase of rapid accretion at M_* ≃ 45 ${M}_{\odot }$. The star quickly and substantially bloats up again, because the accretion rate exceeds ∼10⁻² ${M}_{\odot }\;{{\rm{yr}}}^{-1}$, and only begins to contract after the accretion rate falls during a subsequent quiescent phase (also see Figure 11(a)). The abrupt expansion depicted in detail in Figures 10 and 11 occurs for M_* ≃ 274 ${M}_{\odot }$; it is the third of such events. Figure 12 shows that, even during the supergiant stage, an inner core with a large fraction of the total mass continues to contract while releasing gravitational energy. Only a surface layer with a tiny fraction of the stellar mass participates in the bloating up, trapping a part of the energy transferred from the contracting core. We note that hydrogen burning occurs in the stellar center for M_* ≳ 180 ${M}_{\odot }$ and conclude that rapid mass accretion will even trigger the bloating of an H-burning star and extinguish its stellar UV emissivity.

Zoom In Zoom Out Reset image size

The H ii region begins to recombine each time it loses its exciting UV source. Because the recombination timescale varies as ${t}_{{\rm{rec}}}\propto 1/{n}_{{\rm{H}}{\rm{II}}}$, hydrogen recombination occurs in the densest region first. Within the bipolar H ii region, the radial density gradient is established by the photoevporation flow with geometrical dilution, which roughly follows

$$\begin{eqnarray}{n}_{{\rm{H}}{\rm{II}}} & \sim & {10}^{4}\;{{\rm{cm}}}^{-3}\times {\left(\displaystyle \frac{r}{{10}^{4}{\rm{au}}}\right)}^{-2}\\ & & {\left(\displaystyle \frac{{v}_{{\rm{evp}}}}{100\mathrm{km}{{\rm{s}}}^{-1}}\right)}^{-1}\left(\displaystyle \frac{{\dot{M}}_{{\rm{evp}}}}{{10}^{-3}\;{M}_{\odot }\;{{\rm{yr}}}^{-1}}\right),\end{eqnarray} \tag{ 25 }$$

where ${\dot{M}}_{{\rm{evp}}}$ and v_evp are the photoevaporation rate and flow velocity, normalized with typical values (e.g., Hollenbach et al. 1994; Tanaka et al. 2013). The structure of recombining ionized gas thus resembles a "recombination front," starting at the dense inner regions and propagating outward through the bipolar H ii region (Figures 10(b), (c)). The recombination timescale for the current case is

$$\begin{eqnarray}&&{t}_{{\rm{rec}}}\sim 100\;{\rm{years}}{\left(\displaystyle \frac{{n}_{{\rm{H}}{\rm{II}}}}{{10}^{4}{{\rm{cm}}}^{-3}}\right)}^{-1},\end{eqnarray} \tag{ 26 }$$

which matches the fast propagation of the recombination front shown in Figures 10(b)–(c).

Even after the ionized gas has recombined, the hot atomic gas still remains as a "remnant" of the H ii region (or fossil H ii region) for a long time (Figures 10(c), (d)). This recombined atomic gas still flows outward through inertia, but it is no longer replenished via photoevaporation of the disk. The temperature of the expanding recombined atomic gas gradually decreases due to adiabatic cooling and Lyα cooling. About 1.2 × 10⁴ years after the extinction of the UV source (Figure 10(d)), for instance, the remnant neutral gas has ≃1000–5000 K over ∼pc scales. The central regions r ≲ 10⁵ au have the coldest (T ≲ 100 K) gas because of strong adiabatic cooling.

Whereas the signature of a fossil H ii region remains over ∼pc scales for many thousands of years, influences of the stellar UV radiation in the stellar vicinity r ≲ 10⁴ au disappear much quicker. Figure 13(a) shows that the pressure structure of the accretion envelope is greatly modified by a bipolar H ii region if present. The radial pressure gradient created by the photoevaporation outflow spreads over the envelope even beyond the H ii region (also see Section 4.2.1). However, after the stellar UV emissivity shuts off, these pressure structures quickly disappear. At the epoch depicted in Figure 13(c), almost no signatures of UV feedback are present. The M_* ≃ 500 ${M}_{\odot }$ protostar quiescently accretes the gas through the circumstellar disk with negligible stellar radiative feedback.

Zoom In Zoom Out Reset image size

Figure 14 shows the stellar mass growth history for this case. We see that UV feedback does hinder the mass accretion for 2 ≲ (t/10⁴ years) ≲ 4.5, when a bipolar H ii region emerges (also see Figures 11 and 13(a)). However, the accretion flow is no longer halted once the H ii region disappears at t ≃ 4 × 10⁴ years. The gas accumulated in the envelope thus far begins to fall back toward the star. As a result, the stellar mass rapidly increases during the period 4.5 ≲ (t/10⁴ years) ≲ 6, during which the protostar remains in the supergiant stage emitting few ionizing photons.

Zoom In Zoom Out Reset image size

After the star has accreted most of the gas accumulated in the envelope, the accretion rate gradually decreases for t ≳ 6 × 10⁴ years. The protostar then begins to contract, and the stellar UV emissivity dramatically increases again (Figure 11). Figure 10(f) shows that a bipolar H ii region begins to grow again within the hot (T ≃ 5000 K) atomic envelope, the remnant of the previous H ii region that disappeared after the star had inflated. Figure 13(d) shows that the photoevaporation outflow refills the H ii region. The pressure structure of the accretion envelope is accordingly modified by the expanding H ii region. At this point UV feedback begins to regulate the mass accretion as before.

For this case B, however, UV feedback only intermittently operates during the 7 × 10⁴ years considered, so that the stellar mass growth is not efficiently reduced. Figure 14 shows that, in fact, UV feedback only temporarily postpones stellar mass growth. At the end of the simulation, the stellar mass is slightly higher than for the comparison case B-NF,⁷ for which UV feedback has been turned off. This slightly accelerated stellar mass growth can be understood as follows. As described above, the gas accumulated in the envelope outside the disk rapidly falls back toward the star–disk system once the H ii region disappears for t ≳ 4 × 10⁴ years. With the resulting rapid mass supply from the envelope, the disk becomes highly gravitationally unstable. This promotes the gravitational torque operating in the disk, which consequently helps the stellar mass growth.

Case A: Intermittent UV feedback is also seen in case A. Figure 15 shows that variable accretion causes the abrupt expansion of the stellar radius and rapid drop of UV emissivities as in case B. This case exemplifies a characteristic feature of SE with accretion bursts. For instance, let us focus on the short accretion burst of duration ∼10² years occurring at t ≃ 10⁴ years. The high rate of accretion for this burst causes the protostar to inflate, terminating KH contraction. The protostar remains swollen in the supergiant stage for ∼10⁴ years, long after the accretion rate had dipped below ${\dot{M}}_{\ast }\lesssim {10}^{-2}\;{M}_{\odot }\;{{\rm{yr}}}^{-1}$. This is because, whereas the protostar can react almost immediately to a high accretion rate and quickly expands, it can contract only on a KH timescale even if the accretion drops to zero (e.g., Sakurai et al. 2015). Thus, even short episodes of accretion bursts can have a long-lasting impact on the star's evolution. A similar feature is also seen for the second burst event at t ≃ 3 × 10⁴ years. This effect can only be correctly simulated by solving for the stellar interior structure simultaneously with time-dependent accretion histories, as provided by SE-RHD simulations.

Zoom In Zoom Out Reset image size

In our simulations, circumstellar disks normally form perpendicular to the polar axis of the spherical coordinate system because of our choice of the polar axis orientation (see Section 3.1). The orientation of the disk occasionally fluctuates with time, however, especially for cases A and B, for which the spin parameters measured at the cloud scale are relatively low (Figure 1). For instance, Figure 16 shows the density contours of the disk and ionization contours of the bipolar H ii region at different epochs for case A. Although viewing from the same angle, the orientations of the H ii region as well as the disk changes with time. Since the primordial cloud forms through cosmological structure formation and is not necessarily perfectly axisymmetric, material falling onto the disk has different mean angular momentum vectors at different times. Similar phenomena have been also reported for present-day star formation, whereby turbulent motions in molecular clouds easily twist the distribution of the angular momentum (e.g., Bate et al. 2010; Fielding et al. 2015). Due to this disk "precession," the stellar EUV radiation can ionize the envelope gas over rather wide opening angles. The orientation of the disk at any particular time also depends on whether or not UV feedback is included in the simulations. This is not surprising because UV feedback blows away the gas in the polar regions and modifies the pressure structure even near the equator. As a result, UV feedback can modify the accretion history of both angular momentum and mass from the envelope onto the disk.

Zoom In Zoom Out Reset image size

Case C: For our case C, which also shows the intermittent UV feedback in an early stage, we follow the much longer evolution than the above. For cases A and B, mass accretion is still on-going under intermittent UV feedback for 7 × 10⁴ years, because episodes of rapid mass accretion are able to interrupt the KH contraction sufficiently often. Although it is uncertain how long this evolution can continue for these cases, UV feedback will likely fix the final stellar masses, if KH contraction remains uninterrupted for a sufficiently long time. Our case C exemplifies this evolution. After the protostar enters the supergiant stage for M_* ≳ 50 ${M}_{\odot }$, the star again begins KH contraction, after the stellar mass reaches ∼130 ${M}_{\odot }$ at t ≃ 10⁴ years (Figure 6). The bipolar H ii region quickly expands as the star approaches the ZAMS. Mass accretion is not immediately shut off by UV feedback, however, but continues at a much reduced rate. As a result, the star's mass almost doubles by t ≃ 1.2 × 10⁵ years, when the accretion rate falls below 10⁻⁵ ${M}_{\odot }\;{{\rm{yr}}}^{-1}$.

As mentioned in Section 3.1, HR14 have previously studied the evolution of protostellar accretion for the same primordial clouds considered above, but under the assumption of 2D axial symmetry. Figure 17 presents the stellar mass growth histories for both the current 3D and the previous HR14 2D results in comparison. We see that the 3D results roughly match the overall trends of the 2D results, in spite of a number of technical differences between the adopted numerical codes. Comparing the 2D and 3D results for each case, the stellar masses at a given epoch differ by a factor of a few at most. Some differences come from the fact that, for a 2D simulation, the original 3D structure of a natal cloud is modified to fit the 2D axisymmetric grid by averaging physical quantities over the azimuthal direction. Moreover, the 2D results also depend on the somewhat ad-hoc use of an α-viscosity to model the gravitational torque operating in self-gravitating disks.

Zoom In Zoom Out Reset image size

Figure 18 shows the protostellar mass and UV emissivities as a function of stellar mass for several representative cases. For case D, the evolution is quite similar between 3D and 2D cases. The qualitative differences appear for cases A and B, where both the stellar radius and EUV emissivity abruptly change for ${M}_{\odot }$ ≳ 100 ${M}_{\odot }$ because of the variable accretion encountered in 3D (also see Section 4.2.2). In our previous 2D cases, by contrast, there was an intermediate stage between the supergiant and ZAMS stage, where the stellar radius gradually increases with the mass (named the "P2" path in HR14). SE calculations demonstrate that this behavior is generally present, when the accretion rates are in the range of $4\times {10}^{-3}\;{M}_{\odot }\;{{\rm{yr}}}^{-1}\lesssim {\dot{M}}_{\ast }\lesssim {10}^{-2}\;{M}_{\odot }\;{{\rm{yr}}}^{-1}$ (e.g., Omukai & Palla 2003). For 2D cases, this stage appears because the accretion rates smoothly evolved with time. In 3D, however, the highly variable accretion rates do not remain in this narrow range very long. The intermediate P2 stage thus does not occur in our 3D cases. Instead, the protostar evolves back and forth between the supergiant and ZAMS stages depending on details of the variable accretion history. This effect reduces the stellar EUV emissivity for M_* ≳ 100 ${M}_{\odot }$ in comparison with the 2D cases, as episodes of rapid accretion bring the star back to the supergiant stage several times.

Zoom In Zoom Out Reset image size

Figure 19(a) displays the stellar masses at the end of the 3D simulations as a function of the infall rates measured at the cloud (Jeans) scale (see Section 3.1), demonstrating that the 3D results still show a similar correlation as found in 2D. Considering the fact that the masses for cases A and B are only lower limits, the correlation found in 3D appears to have a steeper slope than for the 2D results (Equation (24)), albeit with much poorer statistical significance (five cases in 3D). This does not necessarily agree with Susa et al. (2014), who do not find this correlation in their 3D RHD simulations. We discuss this issue in Section 5.2 separately.

Zoom In Zoom Out Reset image size

As Figure 19(b) shows, primordial clouds taken for cases A–C offer rather typical initial conditions for primordial star formation. Our results suggest that the very massive stars, even above the mass range of the pair-instability supernovae (PISNe) progenitors 120 ${M}_{\odot }$ ≲ M_* ≲ 240 ${M}_{\odot }$ (Yoon et al. 2012), can form for such typical cases. A significant number of ∼10³ ${M}_{\odot }$ black holes (BHs) might be produced as remnants of these very massive primordial stars. We shall discuss possible observational signatures of such massive primordial stars in Section 5.4.

Although our 3D simulations show a new 3D effect of the intermittent UV feedback, the ansatz described in Section 2.1.2 can bring some resolution-dependence which can affect the evolution quantitatively. It is thus necessary to explore how resolution affects the results of our simulations. First, we shall focus on case B described in Section 4.2.2, for which UV feedback becomes important when the stellar mass exceeds 200 ${M}_{\odot }$ at t ≳ 2 × 10⁴ years (Figure 14). For this case, we study effects of increasing the grid resolution during early times before UV feedback begins to operate.

We first recalculate the evolution of case B with our standard resolution, turning off UV feedback (case B-NF; see Table 1). We use four model configurations of case B-NF at different epochs for M_* ≃ 0, 40, 70, and 120 ${M}_{\odot }$ as the starting configurations for calculations at higher spatial resolution. Next, we interpolate each of the four starting model configurations onto a grid with double the resolution for each spatial coordinate and follow the evolution for 3000 years. These test cases are named as B-NF-HR2-mX, where the suffix "mX" represents the restarting point at M_* ≃ X ${M}_{\odot }$ (see Table 1). The duration of 3000 years is comparable to the Kepler orbital period P_Kepler ≃ 1600 year for typical values of stellar mass M_* ≃ 50 ${M}_{\odot }$ and disk radius ∼500 au (c.f. Equation (19)). The accretion envelope surrounding the disk typically has a density ≲10⁹ ${{\rm{cm}}}^{-3}$, for which the free-fall timescale is ≳3 × 10³ years. Mass infall from the envelope onto the disk at the higher resolution hardly affects the evolution, whereas cooling, fragmentation, and the associated mass and angular momentum transport within the disk occur on much shorter timescales and therefore should be affected by varying spatial resolution.

Figure 20 compares the stellar mass growth and accretion histories for the examined cases. We see that mass accretion becomes more strongly variable with higher resolution; i.e., the accretion rates vary over shorter timescales and with much higher variations: ${10}^{-4}\;{M}_{\odot }\;{{\rm{yr}}}^{-1}\lesssim {\dot{M}}_{\ast }\lesssim 0.1\;{M}_{\odot }\;{{\rm{yr}}}^{-1}$ (Figure 20(b)). Furthermore, the mass that the star accretes during the 3 × 10³ years is somewhat larger with higher resolution, independent of the different starting points of the doubly resolved cases (Figure 20(a)). This can be understood with the help of Figure 3, which displays the face-on disk structure at t ≃ 1.2 × 10³ years for the different resolutions. The disk clearly has sharper and finer density structure at higher resolution. This is because radiative cooling is less restricted by our ansatz and the disk is thus more gravitationally unstable. Since gravitational torques are caused by the non-axisymmetric disk structure such as spiral arms, increasing the resolution can lead to more efficient mass and angular momentum transport through the disk, which results in a higher stellar mass.

Zoom In Zoom Out Reset image size

We once again double the resolution using case B-NF-HR2-m0 at t ≃ 600 years as a starting model, a time at which M_* ≃ 20 ${M}_{\odot }$, and recalculate the evolution (case B-NF-HR4-m20). We repeat these higher resolution simulations with a larger cooling limit f_limit = 24 (case B-NF-HR4-m20-fc; see Section 2.1.2) and a larger sink cell radius r_min = 52 au (case B-NF-HR4-m20-lsk). The resolution is now 4 times higher than the original case B-NF. Because of the high computational cost, the highest resolution simulations are evolved for only 10³ years. Figure 21 compares the stellar mass growth and accretion histories of these quadruple resolution cases with the results of the lower resolution simulations. The resolution dependences discussed above appear again more prominently.

Zoom In Zoom Out Reset image size

Case B-NF-HR4-m20's accretion history displays lots of spiky features, representing multiple accretion burst events (Figure 21(b)). The stellar mass growth over 10³ years is correspondingly greater at higher resolution (Figure 21(a)). Figures 3(c) and 22 show that, at this high resolution, one of the spiral arm breaks up and creates two fragments. These and Figure 23(a) show the distributions of the Toomre Q-parameter on the equator, defined by

$$\begin{eqnarray}&&Q=\displaystyle \frac{{c}_{s}{\rm{\Omega }}}{\pi G{\rm{\Sigma }}},\end{eqnarray} \tag{ 27 }$$

where c_s is the local sound speed and Σ is the vertical column density. Because the spiral arms have Q ≲ 1, gravitational instability makes them prone to fragmentation. This feature is not seen at lower resolution for the same epoch (Figures 3(a), (b)). The evolution for this case is quite similar to that shown in Vorobyov et al. (2013), who study disk fragmentation using 2D face-on (thin disk) simulations at even higher resolution than ours.

Our results suggest that the stellar mass growth is not hindered, but rather accelerated due to disk fragmentation. This is because newly created fragments rapidly migrate inward through the disk and ultimately accrete onto the star. This behavior is shown in Figure 22, tracking the orbital motion of the fragments seen in Figure 3(c). Both fragments rapidly migrate inward and reach the star in less than 100 years. Recent 3D simulations by Greif et al. (2012) also report a similar rapid migration in highly gravitationally unstable disks. Greif et al. (2012) conclude that the typical timescale for such migrations is roughly the local free-fall timescale, the shortest timescale for any physical process driven by gravity. For our case B-NF-HR4-m20, the density near the disk midplane is n ≳ 10¹¹ ${{\rm{cm}}}^{-3}$, for which the free-fall timescale is t_ff ≲ 2.5 × 10² years. The migration timescale is therefore also of the order of the local free-fall timescale in our calculations. The rapid inward migration of the fragments in gravitationally unstable disks is also commonly seen in simulations of present-day star and planet formation (e.g., Vorobyov & Basu 2006; Baruteau et al. 2011; Machida et al. 2011; Zhu et al. 2012). These studies show that rapid migration occurs roughly over the so-called type I migration timescale

$$\begin{eqnarray}&&{t}_{{\rm{mig,}}{\rm{I}}}=\displaystyle \frac{{h}^{2}}{4C\mu q}\displaystyle \frac{2\pi }{{\rm{\Omega }}}\end{eqnarray} \tag{ 28 }$$

(Tanaka et al. 2002), where $C=2.7+1.1\xi$, ξ is the radial gradient of the column density Σ ∝ r^−ξ, h is the disk aspect ratio, μ = πΣr²/M_*, and q = M_f/M_*, where M_f is the fragment mass. Figure 23(b) shows that the column density profile in our disk is well approximated by Σ ∝ r⁻¹, i.e., ξ = 1. The fragment mass M_f is approximated by the mass enclosed within a spatial scale of the most unstable wavelength $\lambda =2{c}_{s}^{2}/G{\rm{\Sigma }}$,

$$\begin{eqnarray}&&{M}_{f}=\pi {(\lambda /2)}^{2}{\rm{\Sigma }}=\displaystyle \frac{\pi {c}_{s}^{4}}{{G}^{2}{\rm{\Sigma }}}.\end{eqnarray} \tag{ 29 }$$

Figure 23(c) shows that the ϕ-averaged fragment masses estimated by Equation (29) span the range from ∼1 ${M}_{\odot }$ to several 10 ${M}_{\odot }$. Because there are large density fluctuations in the disk, the locally estimated fragment mass varies considerably. Figure 24 shows M_f ∼ 0.1–1 ${M}_{\odot }$ along the spiral arms, where the fragmentation often occurs. We also directly calculate the masses of the two fragments seen in Figure 22(a) (marked as 1 and 2), by summing up the masses in the cells with the column density higher than 5000 g cm⁻² around each fragment. The resulting values, M_f1 ≃ 0.4 ${M}_{\odot }$ and M_f2 ≃ 0.9 ${M}_{\odot }$, are in rough agreement with the estimated values in Figure 24 (see also Figure 23(c)).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 23(d), we plot the type I migration timescale given by Equation (28). Considering the fragment masses discussed above, we estimate that the corresponding migration timescales should be around ∼100 years, which are shorter than or comparable to the Kepler orbital periods. The rapid inward migration seen in our simulations can be well explained by the above evaluation. In addition to the two infalling fragments shown in Figures 3(c) and 22, several similar events were observed during the simulation duration of 10³ years. Some of the accretion bursts seen in Figure 21, although not all, are caused by such events.

As mentioned earlier, we have also performed quadruple resolution simulations with different numerical parameters. For case B-NF-HR4-m20-fc, for instance, we adopt f_limit = 24 instead of the fiducial value f_limit = 12, where f_limit sets the cooling limit threshold (Equation (10)). As described in Section 2.1.2, a larger value for f_limit sets a more stringent limit on cooling. Figure 21 shows that the evolution for this case is quite similar to that for the double resolution case B-NF-HR2-m0, verifying that the resolution dependences examined above really come from the more efficient cooling allowed with higher resolution. For case B-NF-HR4-m20-lsk, on the other hand, we adopt a larger sink r_min = 52 au instead of the fiducial value r_min = 30 au. The stellar mass growth history for this case roughly traces that for case B-NF-HR4-m20, although with less variable accretion rates. The shortest timescale of accretion variability is roughly given by the Kepler orbital timescale near the radial inner boundary at r = r_min.

Next we investigate the effects of varying the grid resolution for evolutionary stages when UV feedback becomes important. For this purpose, we take case C (see Section 4.2.2) as the basis for comparison. Unlike case B considered above, UV feedback operates earlier and at lower stellar masses for case C (e.g., Figures 4 and 5). This enables us to follow the evolution with UV feedback for a longer period of time at a reasonable computational cost. We double the resolution for the initial configuration, M_* = 0 ${M}_{\odot }$, and follow the evolution for 5.5 × 10⁴ years (case C-HR2).

Figure 25 displays the stellar mass (a) accretion rate (b), stellar radius (c), and UV emissivity (d) as a function of time for both case C and C-HR2. In agreement with what we have shown in Section 4.3.1, the accretion rate is more highly variable at the higher resolution (Figure 25(b)). After the different accretion histories over 5.5 × 10⁴ years, however, the final masses differ only by ∼10%. This is because the current evolutionary timescale is much longer than the local Kepler timescale within the disk (Equation (19)), so that we are no longer considering the short term evolution of a nearly isolated disk; the mass supply from the envelope ultimately controls stellar growth via mass accretion through the disk. Since the stellar radius abruptly changes with the variable mass accretion, the UV emissivity also displays very large fluctuations (Figures 25(c), (d)). With the doubled resolution, in particular, the protostellar bloating and subsequent KH contraction are repeated several times because of multiple accretion bursts. The extinction and re-formation of bipolar H ii regions also occurs accordingly. In comparison with the original case C, UV feedback begins to limit stellar mass growth earlier at M_* ≃ 100 ${M}_{\odot }$, but the star nevertheless accretes a larger amount of gas (≃150 ${M}_{\odot }$) after the onset of significant UV feedback. The overall result is a slightly higher stellar mass at the end of the simulation than for the original case C.

Zoom In Zoom Out Reset image size