Formation of Massive Primordial Stars: Intermittent UV Feedback with Episodic Mass Accretion
Takashi Hosokawa, Shingo Hirano, Rolf Kuiper, Harold W. Yorke, Kazuyuki Omukai, Naoki Yoshida
cc (cid:13) Draft version September 6, 2018
Preprint typeset using L A TEX style emulateapj v. 5/2/11
FORMATION OF MASSIVE PRIMORDIAL STARS:INTERMITTENT UV FEEDBACK WITH EPISODIC MASS ACCRETION
Takashi Hosokawa , Shingo Hirano , Rolf Kuiper , Harold W. Yorke ,Kazuyuki Omukai , Naoki Yoshida Draft version September 6, 2018
ABSTRACTWe present coupled stellar evolution (SE) and 3D radiation-hydrodynamic (RHD) simulations of theevolution of primordial protostars, their immediate environment, and the dynamic accretion historyunder the influence of stellar ionizing and dissociating UV feedback. Our coupled SE-RHD calculationsresult in a wide diversity of final stellar masses covering 10 M (cid:12) (cid:46) M ∗ (cid:46) M (cid:12) . The formationof very massive ( (cid:38) M (cid:12) ) stars is possible under weak UV feedback, whereas ordinary massive(a few × M (cid:12) ) stars form when UV feedback can efficiently halt the accretion. This may explainthe peculiar abundance pattern of a Galactic metal-poor star recently reported by Aoki et al. (2014),possibly the observational signature of very massive precursor primordial stars. Weak UV feedbackoccurs in cases of variable accretion, in particular when repeated short accretion bursts temporarilyexceed 0 . M (cid:12) yr − , causing the protostar to inflate. In the bloated state, the protostar has lowsurface temperature and UV feedback is suppressed until the star eventually contracts, on a thermaladjustment timescale, to create an H ii region. If the delay time between successive accretion burstsis sufficiently short, the protostar remains bloated for extended periods, initiating at most only shortperiods of UV feedback. Disk fragmentation does not necessarily reduce the final stellar mass. Quitethe contrary, we find that disk fragmentation enhances episodic accretion as many fragments migrateinward and are accreted onto the star, thus allowing continued stellar mass growth under conditionsof intermittent UV feedback. This trend becomes more prominent as we improve the resolution of oursimulations. We argue that simulations with significantly higher resolution than reported previouslyare needed to derive accurate gas mass accretion rates onto primordial protostars. Subject headings: cosmology: theory – early universe – galaxies: formation – stars: formation –accretion INTRODUCTION
Modern cosmology predicts that the early universe isdominated by massive ( M ∗ (cid:38) M (cid:12) ) primordial stars(e.g., Bromm et al. 2009; Greif 2015). The formationof such massive primordial stars is thought to beginin mini-haloes at z (cid:39) −
30 (e.g., Abel et al. 2002;Yoshida et al. 2003), forming tiny embryo protostars( M ∗ ∼ − M (cid:12) ) after the gravitational collapse of natalclouds (e.g., Omukai & Nishi 1998; Yoshida et al. 2008).These protostars grow in mass by accreting surround-ing gas (e.g., Omukai & Palla 2003; Tan & McKee 2004).Their evolution in the accretion stage and how they inter-act with their environment determine the characteristicsof the emerging massive stars.Recent theoretical studies suggest that the typical fi-nal stellar mass will be somewhat smaller than that ofthe cloud mass, ∼ M (cid:12) (e.g., Bromm et al. 2002). Research Center for the Early Universe, the Universityof Tokyo, Tokyo 113-0033, Japan; [email protected], [email protected] Department of Physics, School of Science, the University ofTokyo, Tokyo 113-0033, Japan University of T¨ubingen, Institute of Astronomy and Astro-physics, Auf der Morgenstelle 10, D-72076 T¨ubingen, Germany Jet Propulsion Laboratory, California Institute of Technol-ogy, Pasadena CA 91109, USA Astronomical Institute, Tohoku University, Sendai 980-8578,Japan Kavli Institute for the Physics and Mathematics of the Uni-verse (WPI), Todai Institutes for Advanced Study, the Univer-sity of Tokyo, Chiba 277-8583, Japan
Stellar radiative feedback is one of the crucial processeswhich regulate the mass accretion (e.g., Omukai & In-utsuka 2002; McKee & Tan 2008). Hosokawa et al.(2011, 2012b, hereafter HS11 and HS12) show that,with coupled stellar evolution (SE) and 2D axisymmetricradiation-hydrodynamic (RHD) simulations, UV feed-back eventually halts the accretion of material onto theprotostar. As the stellar UV emissivity rises with increas-ing stellar mass, a bipolar H ii region forms and dynam-ically expands, first slowing, then reversing the inwardflow of gas in the accretion envelope. The circumstellaraccretion disk also loses gas via photoevaporation, beingexposed to the stellar UV radiation. Hirano et al. (2014,hereafter HR14) have greatly extended the above workto follow the evolution in more than one hundred pri-mordial clouds taken from 3D cosmological simulations.Although UV feedback eventually shuts off mass accre-tion in general, the resulting final stellar masses havegreat diversity, spanning 10 M (cid:12) (cid:46) M ∗ (cid:46) M (cid:12) (seealso Susa et al. 2014; Hirano et al. 2015, hereafter HR15).Whereas the 2D simulations have revealed the dy-namic nature of stellar UV feedback, 3D simulationsshow an additional aspect of mass accretion through self-gravitating circumstellar disks (e.g., Stacy et al. 2010;Clark et al. 2011; Greif et al. 2011, 2012; Vorobyov et al.2013; Machida & Doi 2013). In particular, these studiesshow that disks readily fragment via gravitational insta-bility. The fate of the fragments seems to have a great di-versity; some will survive, being expelled from the systemvia gravitational ejection, and others will migrate inward a r X i v : . [ a s t r o - ph . GA ] A p r T. Hosokawa et al.through the disk to merge with each other or fall ontothe central star causing accretion bursts. Through suchrecurrent merger and ejection events, binary or multiplestellar systems will finally appear (e.g., Stacy & Bromm2013).The final stellar mass should ultimately be determinedby the complex interplay between UV feedback and thevarious 3D effects described above. Performing SE-RHDsimulations including all of these effects is still challeng-ing, but nevertheless remains the ultimate goal for com-putational studies. Whereas the first 3D simulationsmostly focused on the early evolution for 10 − yearsafter the birth of the protostar, Stacy et al. (2012) followfor longer periods ( (cid:39) × years), namely until UVfeedback begins to operate. They show that a bipolarH ii region grows from a binary system formed via diskfragmentation. Susa (2013) and Susa et al. (2014) fol-low for even longer periods ( (cid:39) years) of evolution,including the effect of photodissociating (FUV) feedbackwith their 3D simulations, and show that mass accretiononto the star is ultimately halted by the feedback effects.The above studies suggest that both UV feedback anddisk fragmentation might co-operatively reduce the finalstellar masses. This is because the accreting gas is di-vided among multiple protostars and the accretion ratesonto each protostar is correspondingly reduced (e.g., Pe-ters et al. 2010). In fact, the mass spectrum derived bySusa et al. (2014) covers a somewhat lower mass rangethan the 2D results in HR14.In this paper, we cast a new light on a different mode ofinterplay between UV feedback and mass accretion in 3D.As described above, highly variable time-dependent massaccretion histories are commonly seen in 3D numericalsimulations. Even if the disk undergoes fragmentation,a large fraction of the fragments do not actually survive.For instance, Greif et al. (2012) show that about 2 / M ∗ (cid:38) − M (cid:12) yr − dramatically changes the proto-stellar structure, causing the abrupt expansion of the star(e.g., Hosokawa et al. 2012a, 2013; Schleicher et al. 2013).Since a star’s UV emissivity strongly varies with its ef-fective temperature or radius, the resulting UV feedbackexpected for accretion bursts potentially show a new fea-ture during the accretion phase.Here, we present the results of coupled SE-RHD nu-merical simulations (SE in 1D; RHD in 3D) to follow thestellar growth and evolution with variable mass accre-tion, including the effects of stellar UV feedback. As inour previous studies HS11 and HR14, we employ a stel-lar evolution code to numerically solve the interior struc-ture under the assumption of spherical symmetry andmass accretion. We follow the long-term ( ∼ years)evolution of the surrounding gas including the effectsof both the hydrogen-ionizing (EUV: hν ≥ . . ≤ hν ≤ . M ∗ (cid:46) M (cid:12) ) stars. On the otherhand, it is also possible to form very massive ( (cid:38) M (cid:12) )stars, whose observational signature might have been dis-covered in the Galactic metal-poor star SDSS J001820.5-093939. (Aoki et al. 2014). In direct contrast to popularbelief, disk fragmentation does not necessarily reduce thefinal stellar mass; episodic accretion and/or successivemergers can actually enhance stellar mass growth. Forthese cases, the rapid and sometimes extreme evolutionof the stellar radius with variable mass accretion causesthe recurrent extinction and formation of H ii regions.Such intermittent UV feedback does not efficiently haltgas accretion onto the protostar until very massive starsfinally form. We also stress that performing the stellarevolution calculations together with 3D RHD simulationsin realtime enables us to follow the evolution consistently.This is an essential feature of our SE-RHD code.The rest of the paper is organized as follows. We firstdescribe the numerical methods and the cases consideredin Sections 2 and 3. The main results are presented inSection 4. We discuss the implications of our results inSection 5 and summarize our conclusions in Section 6. NUMERICAL METHOD
The evolution in the protostellar accretion stage is fol-lowed with a hybrid SE-RHD code by solving the 3DRHD of the accretion flow in the vicinity of a protostarsimultaneously with the evolution of the protostar insidea sink cell of the RHD grid. This approach is basicallythe same as that developed by HS11. The dynamics ofthe accretion flow are followed by performing a 3D RHDsimulation, whereby for each time step the characteristicsof the protostar (mass, radius, intrinsic luminosity, accre-tion luminosity, FUV flux, EUV flux) provide boundaryconditions for the RHD code at the inner edge of the sinkcell, necessary for the radiation transfer, heating/coolingand chemistry modules of the RHD code. The RHDsimulation in turn provides the time-dependent mass ac-cretion rate into the sink cell. The protostar’s evolutionis calculated by numerically solving the interior struc-ture with the additional feature of mass accretion ontothe stellar surface (e.g., Kuiper & Yorke 2013; Vorobyov& Basu 2015). We describe the methodologies for theRHD simulations and for stellar evolution calculationsin Sections 2.1 and 2.2 respectively.
Radiation Hydrodynamics of Accretion Flow
Simulation Overview and Numerical Configuration
We use a modified version of the public multi-dimensional magneto-hydrodynamics code
PLUTO
Table 1. Cosmological Cases Considered
Cases N R · N θ · N φ r max (10 AU) M g , tot ( M (cid:12) ) ∆ t (10 yr) M ∗ , ( M (cid:12) ) M ∗ , ( M (cid:12) )A 96 · ·
64 150 1 . × ∗ · ·
64 150 1 . × ∗ † · ·
64 150 7B-NF-HR2 ‡ -m0 (cid:63) · ·
128 4 1 . × · ·
128 4 0.3B-NF-HR2-m70 128 · ·
128 4 0.3B-NF-HR2-m120 128 · ·
128 4 0.3B-NF-HR4 ‡ -m20 192 · ·
256 0.67 0.1B-NF-HR4-m20-fc (cid:5) · ·
256 0.67 0.1B-NF-HR4-m20-lsk • · ·
256 0.67 0.1C 64 · ·
64 4 1 . ×
12 286.0 160.3C-HR2-m0 128 · ·
128 4 5.5D 64 · ·
64 4 490 8 67.4 55.5D-NF 64 · ·
64 4 8D-DF † · ·
64 4 8E 96 · ·
64 150 2 . × Col. 2: cell numbers, Col. 3: position of the radial outer boundary, Col. 4: total mass contained in the computational domain, Col. 5: timeduration followed Col. 6: final stellar masses in the current work, Col. 7: final stellar masses obtained in 2D SE-RHD simulations (Hiranoet al. 2014) ∗ : Stellar masses at the end of the simulations, t = 7 × years. † : Suffix “NF” represent casess with no UV feedback and “DF” with only dissociation (FUV) feedback. ‡ : Suffix “HR2” represents cases with doubled and “HR4” with quadrupled spatial resolutions. (cid:63) : Suffix “mX” represents the points for the restart after doubling the resolution when M ∗ (cid:39) X M (cid:12) . (cid:5) : Stringent limit for the cooling with f limit = 24 (see Section 2.1.2). • : Larger sink cell with r min (cid:39)
52 AU ( r min = 30 AU for all other cases). The governing equations of hydrodynamics are ∂ρ∂t + ∇ · ( ρ v ) = 0 , (1) ∂ ( ρ v ) ∂t + ∇ · ( ρ v ⊗ v ) = − ρ ∇ Φ − ∇ p, (2) ∂e∂t + ∇ · ( e v ) = − p ∇ · v + Γ − Λ , (3) p = ( γ − e, (4)where ρ and p are the gas mass density and pressure, v the 3D velocity vector, Φ the gravitational potential, e the gas thermal energy density, Γ and Λ the heatingand cooling rates per volume respectively. We allow theadiabatic exponent γ to vary depending on the chem-ical abundances and gas temperature (e.g., Omukai &Nishi 1998). We calculate the energy source term Γ − Λwith the same thermal processes as considered in HS11(Table S1 in HS11, but see Section 2.1.3 below for somedifferences). Note that equation (2) does not include theviscosity terms which enabled angular momentum trans-port under 2D axial symmetry (e.g., α -viscosity, Shakura& Sunyaev 1973). Instead, gravitational torques inducedby the non-axisymmetric disk structure (e.g., spiral arms,clumps) operate in 3D.As in Kuiper et al. (2011), we solve the above equationsin spherical coordinates. The spatial grids for the polar( θ ) and azimuthal ( φ ) directions are homogeneously dis-tributed over 0 ≤ θ ≤ π and 0 ≤ φ ≤ π . The angularresolutions are thus ∆ θ = π/N θ and ∆ φ = 2 π/N φ , where N θ and N φ are the grid numbers. We set ∆ θ = ∆ φ bytaking N φ = 2 N θ (Table 1). The radial grid size ∆ r logarithmically increases with increasing r as∆ r ( r ) = r (10 f − , (5) where f ≡ log( r max /r min ) /N r , N r is the number of ra-dial grid cells, and r min ( r max ) is the radial inner (outer)boundary of the computational domain.In order to avoid very short timesteps required for thedense gas near and within the star, we do not resolvethe central structure. We instead assume a fixed cen-tral spherical sink cell of radius r min , within which asingle star grows. The sink cell is assumed to be semi-permeable, meaning accretion flow into the sink is al-lowed, but there is no outflow. We evaluate mass accre-tion rates onto the star by measuring the gas inflow ratesinto the sink (also see Section 2.2). We set the sink radius r min = 30 AU for our fiducial cases. Although the radiusof an accreting protostar is generally much smaller thanthis sink size when ˙ M ∗ (cid:46) − M (cid:12) yr − (e.g., Omukai& Palla 2003), an accreting star expands rapidly when˙ M ∗ (cid:38) − M (cid:12) yr − , and inflates up to red supergiantradii. Hosokawa et al. (2012a) derive the mass-radiusrelationship for such a “supergiant protostar” accretingmaterial above the critical rate, R ∗ ( M ∗ ) (cid:39) .
14 AU (cid:18) M ∗ M (cid:12) (cid:19) / , (6)which is comparable to our sink cell size.In the following, we will see that accreting protostarsenter the supergiant stage for many of the examinedcases. Even with the sink size r min = 30 AU, our spheri-cal coordinate system offers relatively high spatial res-olution in the vicinity of the sink, (∆ r × r ∆ θ ) min (cid:39) .
57 AU × .
12 AU for fiducial cases. This is even betterthan the finest resolution used in HS11 ∼
10 AU, so thatthe disk scale height is always resolved with several gridcells even in the innermost regions.With this sink-cell technique, the gravitational poten-tial Φ has two components: one arising from the central T. Hosokawa et al.star Φ ∗ and one from the gas’ self-gravity Φ g ,Φ = Φ ∗ + Φ g . (7)The specific gravitational force due to the central star is − ∇ Φ ∗ = − GM ∗ r e r , (8)where e r is the radial unit vector. The potential forthe gas’ self-gravity Φ g is obtained by solving Poisson’sequation ∆Φ g = 4 πGρ. (9)We make use of the Poisson solver developed in Kuiperet al. (2010a), which has been well tested and used withthe PLUTO code.
Spatially Resolving the Jeans Length
Truelove et al. (1997) have suggested that, in order toavoid artificial fragmentation in numerical simulations,the Jeans length has to be spatially resolved by at leastfour cells. Recent studies argue that even higher resolu-tions are needed to correctly capture the gravity-driventurbulence that generally appears in the star formationprocess (e.g., Federrath et al. 2011; Turk et al. 2012;Meece et al. 2014), suggesting that the Jeans lengthshould be resolved with 32-64 cells. Unfortunately, itis computationally prohibitive to consider the long-term( ∼ years) evolution while fully satisfying this condi-tion. We adopt the following ansatz instead:For each cell and each timestep, we calculate the localJeans length λ J and monitor the ratio of λ J to the greaterof ∆ r and r ∆ θ . When this ratio falls below a certainfixed number f limit , we artificially squelch the coolingfunction in the energy equation, i.e.,Λ → λ J max(∆ r, r ∆ θ ) < f limit . (10)We normally take f limit = 12 as our fiducial value. Inpractice, for a smooth transition, we evaluate the ratio ξ ≡ f limit max(∆ r, r ∆ θ ) /λ J and multiply the cooling rateby a suppressing factor when ξ exceeds unity, C limit = exp (cid:34) − (cid:18) ξ − . (cid:19) (cid:35) , (11)which falls below 10 − for ξ (cid:38) . r, r ∆ θ ) overallmeans that the condition in equation (10) is satisfied forfewer cells and cooling operates without suppression fora larger number of cells. Because the disk becomes moreunstable with more efficient cooling, disk fragmentationoccurs more readily in simulations at higher grid resolu-tion, which is also seen in other studies (e.g., Machida& Doi 2013). We study effects of varying the spatialresolution in Section 4.3 and show that, contrary to thepredictions of previous studies, more fragmentation doesnot necessarily reduce the final stellar masses, but ratherfacilitates the formation of massive stars by inducing in-termittent UV feedback more often. Radiation TransferStellar UV radiation — We consider different kinds of ra-diation which play different roles in our simulations. Thestellar UV flux is the primary source of UV feedbackwhich potentially halts the mass accretion onto stars andfixes the final stellar masses. We solve the transfer ofionizing (EUV) photons emitted by a protostar with afrequency-averaged approximation,1 r ∂ ( r F EUV ) ∂r = − n (1 − x ) σ EUV F EUV , (12)where F EUV is the photon number flux above the Lymanlimit, n the particle density, x the degree of ionization,and σ EUV the photoionization cross section per parti-cle, which is a function of the stellar effective tempera-ture T eff . Although the frequency-dependent UV trans-port 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 aeffect largely modifies the dynamics especially with densemedia in the vicinity of the star. Equation (12) is easilyintegrated F EUV = S EUV πr exp( − τ EUV ) , (13)where S EUV is the stellar EUV emissivity, and τ EUV theoptical depth given by τ EUV = (cid:90) rr min n (1 − x ) σ EUV dr (cid:48) (14)(We implicitly assume no UV absorption within the sinkcell).We solve equations (13) and (14) making use of theadopted spherical grid (Section 2.1.1), i.e., along the ra-dial cells having the same angular coordinate ( θ j , φ k ).The obtained EUV flux F EUV is used to provide thephotoionization rate in the chemical rate equations andthe associated heating rate in the energy equation. Weconsider the ionization of hydrogen only, treating heliumas “hydrogenic”, i.e., singly ionized within an H ii region.For the purposes of our investigation, this approximationis not critical, but certainly could influence the precisenumerical values of our results.Within an H ii region there is a component of diffuseEUV radiation emitted by the recombining gas, namelyrecombination of hydrogen directly into the ground stateand recombinations of helium into the lowest states.HS11 have separately solved the transfer of diffuse hy-drogen recombination photons using the flux-limited dif-fusion (FLD) approximation (Levermore & Pomraning1981; Richling & Yorke 1997). However, UV feedbackon the accreting material is mostly caused by the directEUV stellar flux (also see Tanaka et al. 2013). Thus, weonly solve for the transfer of the direct stellar EUV com-ponent and use the “on-the-spot” approximation (eachrecombination of hydrogen directly into the ground statecauses the immediate ionization of a hydrogen atom inthe vicinity). In praxis this is done by using a hydrogenrecombination coefficient that excludes recombinationsdirectly into the ground state.As for the FUV radiation transfer, we only consider thestellar direct component, by exactly the same method us-ing the self-shielding function (Draine & Bertoldi 1996)assive Primordial Stars 5as in HS11. We recognize that by excluding the diffusecomponents of EUV and FUV, one will encounter ”shad-owing” effects caused by UV-opaque clumps and the ac-cretion disk itself. Our previous comparison studies withthe 2D SE-RHD code have shown that this simplificationdoes not significantly change the accretion history or fi-nal masses of the stars formed, which is the primary goalof this investigation. Molecular hydrogen line emission — We follow the methodused in HS11 for calculating the cooling rate via H lineemission; we sum up the contributions from all “allowed”quadrupole transitions between the rotational and vibra-tional levels of J = 0 −
20 and v = 0 −
2. To estimate thetrapping effect for each transition, we approximate its op-tical depth τ by the local Jeans length λ J , i.e., τ = αλ J ,where α is the transition’s absorption coefficient. This isa crude approximation compared to recently developednovel (and much more computationally expensive) meth-ods (e.g., Clark et al. 2012; Hirano & Yoshida 2013; Greif2014; Hartwig et al. 2015). In our current simulations,however, it is not reasonable to model the detailed ther-mal disk structure, because of our limitations in followingthe exact thermal state of the densest gas (Section 2.1.2).For a given spatial resolution, our treatment underesti-mates the cooling and thus leads to less gravitationallyunstable disks. It therefore affects the quantitative be-havior of disk fragmentation, such as the exact numberand sizes of the emerging fragments. To assess the impactof the details of disk fragmentation on our final results,we vary our spatial resolution (see Section 4.3). Continuum emission — As described in HS11, continuumcooling via H − free-bound emission is important for thegas infalling onto the disk surface. HS11 have imple-mented this cooling process by solving continuum radi-ation transfer using the FLD approximation. In their2D simulations, however, the H − continuum emissionis always optically thin throughout the evolution. Thisknowledge allows us to use a simpler approximation forthe current 3D simulations; we simply subtract the H − bounding energy 0.755 eV for each formation reaction H+ e → H − + γ (R10 in Table 2) from the gas’ inter-nal energy. We have also implemented this method inour 2D SE-RHD code and confirmed the validity of thisapproximation.We do not include the cooling via H collision-inducedcontinuum emission (CIE, Omukai & Nishi 1998). Thisis because CIE only becomes important for the densemolecular gas with n (cid:38) cm − , which only occa-sionally occurs in our simulations. Such dense gas nor-mally lies close to the central star, primarily within thesink cell. Only with the enhanced spatial resolution con-sidered in Section 4.3 below, however, does the densityattain ∼ cm − in the densest parts of disks. Notethat we do not pursue the exact thermal state of suchtiny dense parts anyway, because of our inability to fol-low the long-term evolution at sufficiently high spatialresolution (see discussion in Section 2.1.2). Chemical Network for Primordial Gas
We have implemented an updated version of the pri-mordial chemistry network used in HS11 and HR14 tothe
PLUTO code. Here, we solve for the non-equilibrium
Table 2. Chemical Reactions IncludedNo. Reactions ReferencesR1 H + e → H + + 2 e 1R2 H + + e → H + γ † H − + H → H + e 3R4 H + H + → H + H 4R5 H + e → † H + H → † → H + H 6R8 2 H + H → → → H − + γ γ → H + + e 2R12 H + γ → → H + + e + H 7R14 H − + e → H + 2 e 1R15 H − + H + → ∗ H − + H + → H +2 + e 4R17 ∗ H − + γ → H + e 10R18 ∗ H + H + → H + γ ∗ H + H → H + H + ∗ H + e → ∗ H + H − → H + H 11 † reaction rates updated with respect to HR14, ∗ additional reactions added with respect to HR14REFERENCES: (1) Abel et al. 1997, (2) Osterbrock & Fer-land 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) John1988, (11) Millar 1991 chemical abundances of H, H , H + , e, H − , and H with-out the approximation of using the equilibrium value forH − (as employed by e.g., Abel et al. 1997, HS11). Thereactions included are summarized in Table 2.As described in Section 3.1, HD molecules will formduring the early collapse in some cases. The resulting HDmolecular cooling reduces the temperature and density inthe accretion envelope, which in turn leads to a some-what slower initial protostellar growth (e.g., Yoshidaet al. 2007, HS12). However, this effect ends during theearly collapse by the time n ∼ cm − and is unimpor-tant thereafter. While we include deuterium chemistryfor cosmological simulations following the early collapsestage (see Section 3.1 below), we do not in the currentnetwork for the accretion stage. We have checked withour previous 2D SE-RHD code with the same setting,starting after the density exceeds ∼ cm − as in HS12,that omitting the deuterium chemistry hardly changesthe results. Stellar Evolution Calculations
We follow the evolution of the central protostar by solv-ing the four stellar structure equations taking account ofmass accretion (e.g., Stahler et al. 1980, 1986): (cid:18) ∂r∂M (cid:19) t = 14 πρr , (15) (cid:18) ∂P∂M (cid:19) t = − GM πr , (16) (cid:18) ∂L∂M (cid:19) t = (cid:15) − T (cid:18) ∂s∂t (cid:19) M , (17) T. Hosokawa et al. (cid:18) ∂s∂M (cid:19) t = GM πr (cid:18) ∂s∂p (cid:19) T (cid:18) LL s − (cid:19) C, (18)where M is the Lagrangian mass coordinate, (cid:15) the spe-cific energy production rate by nuclear fusion, s the spe-cific entropy, and L s the radiative luminosity for an adi-abatic temperature gradient. The coefficient C in equa-tion (18) is given by the mixing-length theory for convec-tive zones and is unity otherwise. Hydrogen burning viapp1-, pp2- and pp3-chains and via the CNO-cycle as wellas helium burning via triple- α reactions are included.In order to construct stellar models by solving theabove equations, we adopt the numerical code STELLAR ,originally developed by Yorke & Bodenheimer (2008),which uses the Henyey method. This differs from codeswe used in HS11 and HR14, which employ shooting meth-ods to solve the same equations (e.g., Omukai & Palla2003; Hosokawa & Omukai 2009). We have confirmedthat, despite a number of technical differences, both nu-merical codes provide essentially the same results for var-ious different accretion histories (e.g., Hosokawa et al.2010, 2013). Here, we adopt
STELLAR , because the codeoffers relatively stable convergence properties for highlyvariable accretion rates. Kuiper & Yorke (2013) im-plemented
STELLAR as their stellar evolution module in
PLUTO to study present-day high-mass star formation as-suming 2D axial symmetry. We make use of their imple-mentation by suitably modifying the material functions(e.g., EOS, chemistry module, and opacity tables) forprimordial star formation.We start the stellar evolution calculations after the sinkmass exceeds 1 M (cid:12) via mass accretion, assuming an ini-tial 1 M (cid:12) polytropic star model obtained by solving theLane-Emden equation with the index n = 1 .
5. A timesequence of stellar models are constructed as the stellarmass subsequently increases at the rate provided by the3D RHD simulations. Varying the characteristics of theinitial model (mass, radius, entropy, etc.) only affects theearly evolution long before stellar UV feedback begins tooperate (e.g., Hosokawa et al. 2010).For the SE calculations, we use accretion rates aver-aged over a short duration: 300 years for the defaultcases and 30 years for cases with enhanced spatial reso-lution. This “smearing” of the instantaneous accretionrate measured at r = r min simulates the transport ofmaterial through the innermost part of the disk enclosedwithin the sink. Such a smearing and delaying timescaleis given by the local viscous timescale, presumably sev-eral times longer than the Kepler orbital period, P Kepler (cid:39) . (cid:16) r
30 AU (cid:17) / (cid:18) M ∗ M (cid:12) (cid:19) − / . (19)Note that accretion variability considered below occurson much longer timescales than our assumed averagingtimescale.Mass accretion onto a protostar generally occursthrough a circumstellar disk with finite angular momen-tum. Effects of disk accretion for the stellar evolutioncalculations can be modeled by considering proper outerboundary conditions (e.g., Palla & Stahler 1992). A sim-ple way is adding the accreted gas to the stellar outer-most zone with the same thermal state as in the atmo-sphere. This supposes an extreme “cold” disk accretion, for which the gas lands on the stellar surface gently, hav-ing enough time to adjust its thermal state to that in thestellar atmosphere. The flow from the disk only coversa small part of the surface, and the other part freely ra-diates into a vacuum. In reality, however, the gas newlyjoining the star will carry some additional heat into theopaque interior. We approximate this “warm” accretionby allowing a fraction of the accretion luminosity to beadvected into the star (e.g., Siess et al. 1997), L ∗ , acc ≡ ηL acc = η GM ∗ ˙ M ∗ R ∗ , (20)where η is a parameter, which we normally set to η =0 .
01. As shown in Hosokawa et al. (2013), however, theresulting stellar evolution is almost independent of vary-ing η , except at early stages. Increasing η from 0 to 1simply leads to earlier bloating of the star. Even for ex-treme “cold” disk accretion ( η = 0), bloating occurs longbefore radiative feedback becomes significant.The stellar EUV and FUV emissivities are calculatedassuming a black-body spectrum S EUV = 4 πR ∗ (cid:90) ∞ . πB ( T eff ) hν dν, (21) S FUV = 4 πR ∗ (cid:90) . . πB ( T eff ) hν dν, (22)(McKee & Tan 2008, HS11), where the effective temper-ature is given by T eff = (cid:18) L ∗ πR ∗ σ (cid:19) / , (23)where σ is Stefan-Boltzmann constant. We do not in-clude the accretion luminosity (1 − η ) L acc (see equation23), which presumably has spectral characteristics differ-ing from a black-body at T eff . This omission has littleeffect on our results, because the accretion luminosity L acc is generally much smaller than the stellar luminos-ity L ∗ at the stage when UV feedback becomes effective.The contribution of the accretion luminosity may be-come large when the accretion rate suddenly rises as in aburst event. For these cases, however, the protostar alsorapidly inflates for peak rates exceeding ∼ − M (cid:12) yr − .The stellar UV emissivity accordingly drops, which di-minishes the UV feedback (also see Section 4.2.2 below).For accretion rate spikes below ∼ − M (cid:12) yr − , the ac-cretion luminosity does not significantly affect the totalluminosity. CASES CONSIDERED
Primordial Cloud Collapse in CosmologicalSimulations
We study the formation of primordial stars in a fullcosmological context. To this end, we first follow theformation and gravitational collapse of a number of pri-mordial clouds utilizing a modified version of 3D SPH/N-body code
GADGET -2 (e.g., Springel 2005; Yoshida et al.2006; Hirano & Yoshida 2013). The procedures for thisstep are well described in our previous investigationsHR14 and HR15, who in total have investigated the evo-lution of more than 1600 clouds in a standard ΛCDMassive Primordial Stars 7
Fig. 1.—
Physical properties of the cosmological cases considered here (colored rhombuses denoting cases A–E; see also Table 1) comparedto the cosmological samples analyzed in Hirano et al. (2015) (cyan dots).
Left column: physical properties of gas clouds measured at theself-gravitating cloud (Jeans) scale, (a) infall rates, (b) masses, and (c) spin parameters as a function of the formation redshifts.
Rightcolumn: same as the left column but for the dark matter halo (virial) scale. The masses and spin parameters presented in (e) and (f) areonly for the DM components. universe. The hierarchical zoom-in and particle-splittingtechniques have been used to follow the cloud evolutionwith sufficiently high spatial resolution in cosmologicalvolumes.We focus on five representative cases A–E taken fromthe cosmological samples of HR14 and HR15. For eachcase, we followed the run-away collapse until the centraldensity reaches (cid:39) − cm − (HR14). Here, we studythe subsequent evolution in 3D, using the method de-scribed in Section 2. For initial conditions, we remap theparticle-based cosmological simulation data onto our 3Dspherical grid. The coordinate origin is set at the den-sity maximum and we allow a single star to grow withinthe assumed sink, normally r <
30 AU. The directionof the polar axis is taken to be that of the angular mo-mentum vector for the gas enclosed in the inner region r ≤ . − . ∼ years after the formation of the pro-tostars. Since accretion histories are predestinated bythe envelope structure set in the early collapse stage,higher stellar masses result from a higher global infallrate < πr ρv r > during the collapse, where v r is the ra-dial component of the velocity. In particular, a good es-timator of the final stellar mass can be constructed usingthe infall rate ˙ M J measured at the self-gravitating cloud(Jeans) scale, defined as the radius r J where the ratio ofthe enclosed gas mass M enc ( r ) to the the local Bonner-Ebert mass M BE reaches its maximum as a function of ra-dius r , i.e., M enc ( r J ) /M BE = max( M enc ( r ) /M BE ). HR15have derived a fitting formula relating the final stellarmasses and cloud-scale infall rates, M ∗ = 250 M (cid:12) (cid:32) ˙ M J . × − M (cid:12) yr − (cid:33) . , (24)where the infall rates are measured when the cloud cen-tral density reaches 10 cm − . Figure 1 shows the diver-sity of the cloud-scale infall rates and relevant quantities:the masses and spin parameters at the cloud (Jeans) andhalo (virial) scales. We see that the five cases A–E con-sidered here well cover the distribution of ∼ ∼ M ∗ (cid:39) M (cid:12) . This suggests the presence of T. Hosokawa et al. s t e ll a r m a ss : M ∗ [ M fl ] (a) AB CDE years ]10 -7 -6 -5 -4 -3 -2 -1 a cc . r a t e s : ˙ M ∗ [ M fl y r − ] (b) Fig. 2.—
Time evolution of the stellar masses (panel a) andaccretion rates (panel b) for cases A (green), B (blue), C (black),D (red), and E (magenta). The time t = 0 corresponds to theformation of embryo protostars. a fiducial initial condition for primordial star formation.Among our five cases, cases A–C are more representativeof such a fiducial condition than cases D and E.Case E provides the lowest final stellar mass M ∗ (cid:39) . M (cid:12) in 2D. For this case, HD line cooling also oper-ates in addition to H cooling in the collapse stage. Thetemperature consequently drops to a few ×
10 K, lim-ited by the cosmic microwave background floor. HR14show that the HD-cooling mode can be triggered whenthe collapse is somewhat slower than for purely sphericalfree-fall. Such a slow collapse can occur, for example,when rotational support decelerates the collapse. In theprotostellar accretion stage, UV feedback typically lim-its the stellar masses to a few tens of M (cid:12) (e.g., HS12;HR14).We assume that the primordial clouds for cases A–E do not experience any external feedback. In reality,however, FUV radiation from nearby stars can alter thethermal evolution during collapse and consequently theresulting final stellar masses (e.g., Oh & Haiman 2002;O’Shea & Norman 2008). HR15 show that, by consid-ering the spatial distribution of FUV fields created byPop III stars at different redshifts, these so-called PopIII.2 stars (e.g., Bromm et al. 2009; O’Shea et al. 2008)become comparable in number to the so-called Pop III.1primordial stars, for which the external feedback is neg-ligible. Interestingly, HR15 show that the final stellarmasses and cloud-scale infall rates still have the samecorrelation among the different sub-populations of III.1and III.2 stars. We examine how this correlation candiffer in 3D, focusing on the interplay between variablemass accretion and UV radiative feedback.In addition to the cosmological cases described above,we perform simulations starting with idealized initialconditions: a homogeneous, zero-metallicity gas cloud insolid body rotation with constant density, temperature,and chemical composition. The results of these idealized !" (cid:3) '"( )%*%+,+-$./%01 %*%/4,5%2 ! (cid:3) '6(%)%*%+,+-$./%01 %*%/7,8%2 ! (cid:3) +999 % : ; (cid:3) <%4%1$ (cid:3) <%8%1$ (cid:3) 'B(%)%*%+,+-$./%01 %*%8C,8%2 ! (cid:3) B=>?DA%%E$A (cid:3) G (cid:3) G (cid:3) G (cid:3) G (cid:3) +,$.8% % +,$./% % +99% % +9%%HI%BD G4 J (cid:3) Fig. 3.—
Gas column density along the polar direction of theself-gravitating circumstellar disks obtained with different spatialresolution (cases B-NF, B-NF-HR2-m0, B-NF-HR4-m20 for panelsa, b, and c, respectively; see text and Table 1). The snapshots areshown at the same epoch t = 1 . × years, but the accreted stel-lar masses differ due to different prior accretion histories. In eachpanel, the black contours mark the positions where the Toomre-Qparameter takes the value of Q = 0 . simulations are discussed in the Appendix. These moresimplified initial conditions allow us to test the effects ofvarying numerical settings in a controlled manner. RESULTS
Overall Results
We first summarize the simulation results for our fidu-cial cases A–E. As seen in Figure 2-(a), the stellar massassive Primordial Stars 9 !" (cid:3) ' (cid:3) ! (cid:3) ( (cid:3) ) (cid:3) *$+,$-"*.-$%/01%2%3$456"37%8 (cid:3) (cid:3) @$3A>;*<%/B+= (cid:3) *%D%EFGH$IJ%<- L %D%MNJFN%K ! (cid:3) *%D%EFGM$IJ%<- L %D%MEEFO%K ! (cid:3) *%D%EFME$IJ%<- L %D%EJJFM%K ! (cid:3) *%DMFJH$IJ%<- L %D%EHFJ%K ! (cid:3) *%D%EFEP$IJ%<- L %D%HHFE%K ! (cid:3) !!""" $ % (cid:3) Q%ERR%/B+= (cid:3) @$3A>;*<%%26;?6$-%@"3.$ (cid:3)
STT%-$?;A: (cid:3)
HFR$IJ%%UHHU%%OJGFP%%ENGFP%%HRFR (cid:3)
EFR$IER%%HFM$IRG%%EFR$IRP%%HFM$IRU%%EFR$IRJ (cid:3)
ERFR%%UFNM%%HFEN%%EFPG%%EFR (cid:3)
Fig. 4.—
Temperature, density, and velocity structure of cases A–E in a plane containing the polar axis at the time when the polarextension of H ii regions first reaches 10 AU. In each panel, the gas temperature is shown to the left and the density distribution tothe right. Velocity vectors v <
10 km s − are depicted by small blue and orange arrows, whereby the color gradation represents themagnitude. Velocities v >
10 km s − are depicted by white arrows with lengths varying in proportion to the magnitude (see the legend inthe lower-right panel). The small 3D magenta contour at the panel center delineates the iso-density surface for 10 cm − , which roughlytraces an accretion disk around the central star. Each panel also shows the elapsed time and stellar mass. !" (cid:3) ' (cid:3) ! (cid:3) ( (cid:3) ) (cid:3) * + %,-./$0%10"234,%% ? %6%+@<8@%> ! (cid:3) ? %6%+778A%> ! (cid:3) ? %6%7<<8+%> ! (cid:3) ? %6%7:8<%> ! (cid:3) ? %6%::87%> ! (cid:3) !!""" $ % (cid:3) CD104,5 (cid:3) )D104,5 (cid:3)
E87% % % % % % (cid:3) Fig. 5.—
Same as Figure 4 but for the chemical structure in the accretion envelope. In each panel, the color image shows the distributionof the hydrogen molecular number fraction in the same vertical plane as in Figure 4. The 3D red contour delineates the structure ofionization fronts, within which the hydrogen atomic fractions is below 0.5. R ∗ [ R fl ] (a) stellar radius A BCDE stellar mass: M ∗ [ M fl ]10 S E UV , S F UV [ s ec − ] (b) UV emissivity BD Fig. 6.—
Evolution of the stellar radius (panel a) and UV emissiv-ities (panel b) as a function of accumulated stellar mass. The differ-ent colors represent the same cases as in Figure 2. In panel (a), thethick gray dotted and long-dashed lines represent the mass-radiusrelationships for zero-age main-sequence stars and for supergiantprotostars (equation 6). In panel (b), the dashed line representsthe stellar ionizing (EUV; hν ≥ . . ≤ hν ≤ . growth histories show a great variety in the ∼ yearsafter the formation of embryo protostars. Ordinary mas-sive stars with M ∗ < M (cid:12) result in cases D and E.In particular, the final stellar mass is only ∼ M (cid:12) forcase E. For cases A–C, on the other hand, the stellarmasses reach several hundred solar masses at the end ofthe simulations.Figure 2-(b) displays the highly variable accretion his-tories for all cases. Such a high degree of variabil-ity is a well-known feature of mass and angular mo-mentum transport through self-gravitating circumstellardisks. Figure 3 presents examples of the density struc-ture of self-gravitating disks seen in our simulations. Thedisk has non-axisymmetric spiral arms, which producethe gravitational torques leading to variable mass accre-tion. In Figure 2-(b), we see that the mass accretionfinally ceases for cases C-E; the mean accretion rates forthe last 10 years having fallen down to ∼ − M (cid:12) yr − due to UV feedback. For cases A and B, however, theaccretion rates are still higher than 10 − M (cid:12) yr − at t = 7 × years. Although the simulations end at thispoint due to high computational costs, the stellar masseswill continue to increase via accretion, because UV feed-back has not yet become fully effective.Figures 4 and 5 shows that, for all cases, bipolar H ii regions form and grow within the larger H photodisso-ciation regions (PDRs) within a few × years afterstellar birth. Although the snapshots look quite similar,the stellar masses at the evolutionary stage of H ii regionformation differ significantly among the cases; the stellarmass is lower in alphabetical order for cases A through E. For cases A–C, it is not until the stellar masses ex-ceed 100 M (cid:12) that the H ii regions appear. For cases Dand E, however, only a few × M (cid:12) stars can create H ii regions.The diversity of stellar masses which first form anH ii region comes from the different protostellar evo-lution with different accretion histories (Fig. 6). Forcases D and E, the protostar contracts from its bloatedstate when M ∗ (cid:38) M (cid:12) ; this is the so-called Kelvin-Helmholtz (KH) contraction phase (e.g., Omukai & Palla2003). The stellar UV emissivities dramatically increaseduring this stage as the effective temperature rises as T eff ∝ R − / ∗ . Stellar contraction continues until the stararrives at the zero-age main-sequence (ZAMS). As a re-sult, the H ii region emerges for M ∗ (cid:39) M (cid:12) . For casesB and C, however, KH contraction is interrupted by theabrupt stellar expansion occurring around M ∗ (cid:39) M (cid:12) .For case A, there is no contraction stage until the stellarmass exceeds 100 M (cid:12) . It is only after the stellar massreaches a few hundred solar masses that the protostarsubstantially contracts. Since the stellar UV emissiv-ity sensitively depends on the effective temperature, theepoch of stellar contraction actually controls when theH ii region first appears.For cases A–C, the protostar enters the supergiantstage with rapid mass accretion ˙ M ∗ (cid:38) . M (cid:12) yr − be-fore KH contraction begins (Section 2). Several studiespresuppose that supergiant protostars only appear in ex-treme cases of primordial star formation, for which veryhigh accretion rates are available (e.g., direct-collapsemodel, see Bromm & Loeb 2003). However, HR14 andHR15 show with 2D SE-RHD simulations that the super-giant stage can appear even for more typical Pop III.1cases when accretion rates are relatively high but notnecessarily extremely high. By contrast, our 3D SE-RHDsimulations now demonstrate that the supergiant stageis a common feature of primordial star formation. This isdue to the highly variable mass accretion onto primordialprotostars caused by self-gravitating disks, for which theaccretion rates temporarily and recurrently exceed thecritical rate for expansion to the supergiant stage. Asa result of such variable mass accretion, in fact, stellarcontraction and re-expansion are repeated several timesfor cases A and B. The volatile evolution of the stellarradius results in strong fluctuations of the stellar UVemissivity. This variability causes recurrent extinctionand re-formation of H ii regions (see Section 4.2.2 be-low). UV feedback only intermittently operates and can-not efficiently stop mass accretion. Very massive starscan potentially form by this mechanism.Although our results are resolution dependent becauseof the ansatz explained in Section 2.1.2, the basic featuresdescribed become more prominent at higher spatial res-olution (Section 4.3). This is demonstrated in Figure 3,showing that higher resolution results in sharper andfiner density substructure within the disk. The resultingmass accretion is more variable with increased spatialresolution, bringing the central protostar to the super-giant stage more often. Disk fragmentation (e.g., Fig. 3-c) does not necessarily reduce the final stellar masses.The emerging fragments cause repeated accretion burstswhen falling onto the star, which can extinguish the stel-lar UV emission and thus ultimately enhances stellarassive Primordial Stars 11 !" (cid:3) ! " ! $ % & ' ( (cid:3) '"(%)%*%+,-$./%01 %*%45,6%2 ! (cid:3) '7(%)%*%8,4$.+%01 %*%//,+%2 ! (cid:3) '9(%)%*%-,5$.+%01 %*%6-,+%2 ! (cid:3) ':(%)%*%;,8$.+%01 %*%6<,+%2 ! (cid:3) =&> (cid:3) ?@@%1$ABCD (cid:3) Fig. 7.—
Time evolution of the temperature distribution in a plane containing the polar axis for case D. The stellar mass and evolutionarytime are given in the insets of panels (a)–(d). The color scale for temperature is the same as in Figure 4. The 3D white and red contoursdelineate the structure of dissociation and ionization fronts, within which the hydrogen atomic and molecular fractions are below 0.5 and10 − , respectively. mass growth. Mass Accretion under the Stellar UV Feedback
Massive Stars with Tens of Solar Masses:UV Feedback Halts the Mass Accretion
We first describe in detail the evolution of cases D andE, which produced ordinary massive stars with tens ofsolar masses. As shown below, the evolution for thesecases is similar to that found in our previous 2D SE-RHDsimulations (HS11; HS12).
Case D — Figure 7 presents the evolution of the tem-perature structure for case D as an H PDR and an H ii region expand into the accretion envelope. After the pro-tostar enters the KH contraction stage (see Fig. 6), thestellar FUV emissivity rises and forms a bipolar H PDR(Fig. 7-a). At the same time, the stellar EUV emissivityincreases and slightly later creates a bipolar H ii regionwithin the larger PDR (Fig. 7-b). The PDR and H ii region continue to grow and expand as shown for theepoch t = 5 × years (Fig. 7-c), but these regionsare still centrally “pinched” because a circumstellar diskblocks the stellar UV radiation close to the disk plane. years ]020406080100120140 s t e ll a r m a ss : M ∗ [ M fl ] Case D
I+D FeedbackD FeedbackNO Feedback
Fig. 8.—
Example of the impact of differing stellar UV feedbackon stellar mass growth. The black solid (case D), magenta solid(case D-DF), and blue bashed lines (case D-NF) represent the casesincluding both the ionizing and dissociating feedback, with thedissociating feedback only, and without UV feedback, respectively.
However, even the gas protected from the UV photonsis disturbed by shocks that propagate into the envelopeahead of the ionization fronts. Comparing Figure 7-(b)and (c), we see that even the gas outside the PDR is2 T. Hosokawa et al. !" (cid:3) !* (cid:3) !+ (cid:3) (/ $ (cid:3) < $6$=0>$; ! (cid:3) < $6$=?@$; ! (cid:3) < $6$@712$; ! (cid:3) Fig. 9.—
Impact of differing stellar UV feedback on the accretionenvelope: (a) no feedback (case D-NF), (b) only photo-dissociatingfeedback (case D-DF), and (c) both photo-ionizing and dissociatingfeedback (case D). The temperature, density and velocity distribu-tions are shown at the same epochs, 5 × years after the birthof the central protostar. Note that the stellar masses differ amongthe cases because of the different mass accretion histories (also seeFig. 8). The illustration is in the same style as in Figure 4. heated up by shock compression. HS11 have shown thatthe mass influx from the envelope onto the disk is hin-dered by this effect. Moreover, the disk loses mass viaphotoevaporation, being exposed to the stellar EUV ra-diation. The H ii region is filled by the photoevaporatingoutflow launched from the disk surface. Ultimately, thecircumstellar disk is destroyed due to photoevaporation,allowing the H ii region to expand in equatorial directions(Fig. 7-d), and accretion onto the star is effectively shut off for M ∗ (cid:39) . M (cid:12) at t = 8 × years.In Figure 8, we plot the evolution of the stellar mass forcase D (solid black line). For comparison, we show thestellar mass growth for case D-NF (blue dashed line), forwhich all UV feedback effects had been turned off. Wesee that mass accretion onto the protostar is significantlyreduced by UV feedback effects. This figure also displaysthe stellar mass growth for a test case (case D-DF) thatonly included photodissociation (FUV) feedback (solidmagenta line). Comparing these mass growth historieswe conclude that FUV feedback does hinder mass accre-tion somewhat, but its impact is much weaker than EUVfeedback.Figure 9 shows the variation of the envelope tempera-ture - density - velocity structure with differing UV feed-back effects. For case D-NF with neither EUV nor FUVfeedback (Fig. 9-a), the disk is surrounded by a warm( T (cid:39) −
500 K) envelope, where compressional heat-ing due to gas infall is balanced by H molecular linecooling. The gas temperature is higher ( T (cid:39) r (cid:46) AU, where H − free-bound emis-sion is the main cooling process. For case D-DF with onlyFUV feedback (Fig. 9-b), the warm envelope is larger, be-cause FUV radiation reduces the abundances of efficientcoolants via H photodissociation and H − photodetach-ment. However, this only occurs in polar regions. Be-cause the bulk of stellar FUV radiation is blocked by thedisk, the envelope structure is hardly affected in equato-rial directions, where the disk mostly accretes gas fromthe envelope. As seen in Fig. 9-c EUV feedback dra-matically changes the envelope structure. In addition tothe bipolar photoevaporation flow within the H ii region,portions of the envelope outside the H ii region are ex-pelled from the computational grid as shocks propagatethrough the envelope and spread the pressure excess ofionized gas.Our results demonstrate that FUV feedback does re-duce mass accretion rates, which qualitatively agreeswith 3D RHD simulations performed by Susa (2013).However, the FUV feedback seen in our simulationsseems to be weaker than in Susa (2013), who claims thatFUV feedback is largely responsible for terminating massaccretion onto primordial protostars. We discuss possiblecauses of this difference in Section 4.2.3 below. Case E — The evolution of case E is qualitatively similarto that for case D, except that HD molecular line emis-sion also contributes to the cooling in the early collapsestage. During collapse, in general, the thermal evolutionof the central homogeneous core strongly affects the ra-dial structure of the accretion envelope. Since the radiusof the core is well approximated by the Jeans length,the envelope temperature and density at a given radius r obey r = λ J ∝ T / ρ − / , i.e., a lower temperatureimplies a lower density. In case E the accretion enve-lope thus has relatively low temperatures and densities,which can been seen in Figure 4. KH contraction beginsat lower stellar masses with correspondingly lower accre-tion rates (Fig. 6). The final stellar mass is also lowerthan for case D, M ∗ (cid:39) . M (cid:12) . Massive Stars with Hundreds of Solar Masses:Intermittent UV Feedback with Variable Mass Accretion assive Primordial Stars 13 !" (cid:3) '"(%)%*%+,-+$.+%/0 %*%34+%1 ! (cid:3) '5(%677%/0 %*%349%1 ! (cid:3) ':(%+77%/0 %*%3;3%1 ! (cid:3) '$(%6,3$.+%/0 %*%+<;%1 ! (cid:3) '=(%6,-$.+%/0 %*%-+4%1 ! (cid:3) + % > : %% ' ; , $. - % ? @ ( (cid:3) 'A(%3,9$.B%/0 %*%B6B%1 ! (cid:3) B,7$.+% % B,7$.B% % B77% % B7 (cid:3) Fig. 10.—
The extinction and re-formation of H ii regions for case B. In each panel, the 3D yellow contour delineates the structure of theionization front (as in Figures 5 and 7 with the red contours), and the color image shows a 2D slice of the temperature distribution (colorscale as in Figure 4). The panels (a)–(f) show a time sequence for t > . × years after the birth of the protostar; the elapsed timeafter the epoch of panel (a) is given for (b)–(f). Note that the box size is more than 50 times larger than in Figure 7. R ∗ [ R fl ] -5 -4 -3 -2 -1 ˙ M ∗ [ M fl y r − ] (a) stellar radius & acc. rates Case B years ]10 S E UV , S F UV [ s ec − ] (b) EUV/FUV emissivity Fig. 11.—
Effects of the variable mass accretion on the evolutionof the stellar radius (panel a) and UV emissivities (panel b) forcase B. In panel (a), the black and magenta lines show the timeevolution of the stellar radius and accretion rates respectively. Inpanel (b), the dashed and dotted lines represent the stellar EUVand FUV emissivities as in Figure 6-(b). The vertical blue dottedline marks the epoch of the stellar inflation at t = 4 . × years(also see Fig. 12). HW Yorke, Jet Propulsion Laboratory, Pasadena using HWYplt V2.4 package
Mass [ M ]100 101 102 R ad i u s [ R ] Halo B D H Case B
Fig. 12.—
Evolution of the stellar interior structure as a functionof increasing mass for case B. The black solid and dashed lines rep-resent the radial positions of the photosphere and mass coordinates80%, 60%, 40%, and 20% of the total stellar mass. The yellow re-gions denote convective zones and white regions radiative zones,each without active nuclear burning. The pink and brown stripesdenote convective and radiative deuterium-burning zones, respec-tively and green depicts the central convective hydrogen-burningzone. Layers with active nuclear burning are indicated if the de-pletion time of the “burning” species is shorter than the lifetime ofa main-sequence star with the equivalent mass. The vertical bluedotted line marks the same epoch as indicated in Figure 11. !" (cid:3) ! " ! $ % & ' ( (cid:3) '"(%)%*%+,-$./%01 %*%+4+%2 ! (cid:3) '5(%)%*%/,-$./%01 %*%464%2 ! (cid:3) '7(%)%*%8,9$./%01 %*%/:-%2 ! (cid:3) ';(%)%*%<,+$./%01 %*%8/9%2 ! (cid:3) =1$ B+ C%'D$E(%F% (cid:3) @>A5$1%;$@ B4 C%'1GHI)(%F% (cid:3)
Fig. 13.—
The evolution of pressure (left-half of each panel) and density (right-half) under the influence of intermittent UV feedback forcase B. Panels (b)–(d) show the same epochs as in Figure 10 (d)–(f) but on a smaller physical scale. The 2D velocity vectors are presentedin the same style as in Figure 4. In each panel, the short dashed line delineates the location of the hydrogen ionization front and the solidgray line depicts the dissociation front, whereby the fractional atomic and molecular abundances are 0.5 and 10 − , respectively. The thinblue lines in the left panels show the pressure contours for P = 1 and 0 . × − dynes cm − . years ]0100200300400500600 s t e ll a r m a ss : M ∗ [ M fl ] Case B
I+D FeedbackNO Feedback
Fig. 14.—
Same as Figure 8 but for case B, i.e., the black solidand blue dashed lines represent the evolution for cases B and B-NF.
Next we describe the results for cases A–C, for whichvery massive ( M ∗ (cid:38) M (cid:12) ) stars form that passthrough a supergiant evolutionary stage. Case B — Figures 10 and 11 show that for case B, vari-able mass accretion leads to abrupt and radical changesof the stellar radius and the associated recurrent extinc-tion and re-formation of H ii regions. A bipolar H ii re-gion, which have formed prior to t (cid:39) . × years(Fig. 10-a), begins to retreat in a direction away fromthe star a few hundred years later (Fig. 10-b,c), shortlyafter the accretion rates rise above ∼ − M (cid:12) yr − andthe protostar rapidly inflated to a supergiant (Fig. 11-a).Due to the protostar’s expansion the stellar EUV emis-sivity accordingly drops by almost ten orders of mag-nitude (Fig. 11-b), initiating the extinction of the H ii region.Figure 12 depicts the evolution of the stellar inte-rior structure for this case B. Although the star ini-tially begins KH contraction when M ∗ (cid:39) M (cid:12) , thecontraction is interrupted by a phase of rapid accre-tion at M ∗ (cid:39) M (cid:12) . The star quickly and substan-tially bloats up again, because the accretion rate exceedsassive Primordial Stars 15 ∼ − M (cid:12) yr − , and only begins to contract after theaccretion rate falls during a subsequent quiescent phase(also see Fig. 11-a). The abrupt expansion depicted indetail in Figures 10 and 11 occurs for M ∗ (cid:39) M (cid:12) ; itis the third of such events. Figure 12 shows that, evenduring the supergiant stage, an inner core with a largefraction of the total mass continues to contract while re-leasing gravitational energy. Only a surface layer with atiny fraction of the stellar mass participates in the bloat-ing up, trapping a part of the energy transferred from thecontracting core. We note that hydrogen burning occursin the stellar center for M ∗ (cid:38) M (cid:12) and conclude thatrapid mass accretion will even trigger the bloating of anH-burning star and extinguish its stellar UV emissivity.The H ii region begins to recombine each time it losesits exciting UV source. Because the recombinationtimescale varies as t rec ∝ /n HII , hydrogen recombina-tion occurs in the densest region first. Within the bipo-lar H ii region, the radial density gradient is establishedby the photoevporation flow with geometrical dilution,which roughly follows n HII ∼ cm − × (cid:16) r AU (cid:17) − (cid:16) v evp
100 km sec − (cid:17) − (cid:32) ˙ M evp − M (cid:12) yr − (cid:33) , (25)where ˙ M evp and v evp are the photoevaporation rate andflow velocity, normalized with typical values (e.g., Hol-lenbach et al. 1994; Tanaka et al. 2013). The structure ofrecombining ionized gas thus resembles a “recombinationfront”, starting at the dense inner regions and propagat-ing outward through the bipolar H ii region (Fig. 10-b,c).The recombination timescale for the current case is t rec ∼
100 years (cid:16) n HII cm − (cid:17) − , (26)which matches the fast propagation of the recombinationfront shown in Figure 10-(b-c).Even after the ionized gas has recombined, the hotatomic gas still remains as a “remnant” of the H ii region(or fossil H ii region) for a long time (Fig. 10-c,d). Thisrecombined atomic gas still flows outward through iner-tia, but it is no longer replenished via photoevaporationof the disk. The temperature of the expanding recom-bined atomic gas gradually decreases due to adiabaticcooling and Ly- α cooling. About 1 . × years afterthe extinction of the UV source (Fig. 10-d), for instance,the remnant neutral gas has (cid:39) − ∼ pcscales. The central regions r (cid:46) AU have the coldest( T (cid:46)
100 K) gas because of strong adiabatic cooling.Whereas the signature of a fossil H ii region remainsover ∼ pc scales for many thousands of years, influences ofthe stellar UV radiation in the stellar vicinity r (cid:46) AUdisappear much quicker. Figure 13-(a) shows that thepressure structure of the accretion envelope is greatlymodified by a bipolar H ii region if present. The radialpressure gradient created by the photoevaporation out-flow spreads over the envelope even beyond the H ii re-gion (also see Section 4.2.1). However, after the stellarUV emissivity shuts off, these pressure structures quicklydisappear. At the epoch depicted in Figure 13-(c), al-most no signatures of UV feedback are present. The M ∗ (cid:39) M (cid:12) protostar quiescently accretes the gasthrough the circumstellar disk with negligible stellar ra-diative feedback.Figure 14 shows the stellar mass growth history forthis case. We see that UV feedback does hinder the massaccretion for 2 (cid:46) ( t/ years) (cid:46) .
5, when a bipolar H ii region emerges (also see Figs. 11 and 13-a). However,the accretion flow is no longer halted once the H ii regiondisappears at t (cid:39) × years. The gas accumulatedin the envelope thus far begins to fall back toward thestar. As a result, the stellar mass rapidly increases duringthe period 4 . (cid:46) ( t/ years) (cid:46)
6, during which theprotostar remains in the supergiant stage emitting fewionizing photons.After the star has accreted most of the gas accumu-lated in the envelope, the accretion rate gradually de-creases for t (cid:38) × years. The protostar then beginsto contract, and the stellar UV emissivity dramaticallyincreases again (Fig. 11). Figure 10-(f) shows that abipolar H ii region begins to grow again within the hot( T (cid:39) ii region that disappeared after the star had in-flated. Figure 13-(d) shows that the photoevaporationoutflow refills the H ii region. The pressure structure ofthe accretion envelope is accordingly modified by the ex-panding H ii region. At this point UV feedback begins toregulate the mass accretion as before.For this case B, however, UV feedback only intermit-tently operates during the 7 × years considered, sothat the stellar mass growth is not efficiently reduced.Figure 14 shows that, in fact, UV feedback only tem-porarily postpones stellar mass growth. At the end of thesimulation, the stellar mass is slightly higher than for thecomparison case B-NF, for which UV feedback has beenturned off. This slightly accelerated stellar mass growthcan be understood as follows. As described above, thegas accumulated in the envelope outside the disk rapidlyfalls back toward the star-disk system once the H ii re-gion disappears for t (cid:38) × years. With the resultingrapid mass supply from the envelope, the disk becomeshighly gravitationally unstable. This promotes the gravi-tational torque operating in the disk, which consequentlyhelps the stellar mass growth. Case A: — Intermittent UV feedback is also seen in caseA. Figure 15 shows that variable accretion causes theabrupt expansion of the stellar radius and rapid dropof UV emissivities as in case B. This case exemplifies acharacteristic feature of stellar evolution with accretionbursts. For instance, let us focus on the short accretionburst of duration ∼ years occurring at t (cid:39) years.The high rate of accretion for this burst causes the pro-tostar to inflate, terminating KH contraction. The pro-tostar remains swollen in the supergiant stage for ∼ years, long after the accretion rate had dipped below For case B-NF, there is a big jump of the stellar mass at t (cid:39) . × years. Although we do not resolve those additionalindividual infalling clumps that would appear at higher resolution(see Figure 3), variable accretion is a general feature even at thedefault resolution. Interestingly, the inner part of the B-NF diskis greatly tilted at t (cid:39) . × years. It is known that externaltidal perturbations from the surrounding non-axisymmetric enve-lope could excite the “bending wave”, which enhances the diskinclination and angular momentum transport (e.g., Papaloizou &Terquem 1995). R ∗ [ R fl ] -5 -4 -3 -2 -1 ˙ M ∗ [ M fl y r − ] (a) stellar radius & acc. rates Case A years ]10 S E UV , S F UV [ s ec − ] (b) EUV/FUV emissivity Fig. 15.—
Same as Figure 11 but for case A. !" (cid:3) '$( (cid:3) *%7%898$:;%+/ = %7%8>;%< ! (cid:3) = %7%3@@%< ! (cid:3) $: ; % & A (cid:3) Fig. 16.—
Example of a precessing accretion disk with a bipolarH ii region as seen for case A. The panels (a) and (b) depict thedensity and ionization structure as viewed from the same anglefor the evolutionary times 2 . × years and 3 × years afterthe birth of the protostar. The dark blue, cyan, and green-graycontours represent cut-outs of the iso-density surfaces for n H =10 , 10 , and 10 cm − , respectively. The red contours delineatethe H ii region. ˙ M ∗ (cid:46) − M (cid:12) yr − . This is because, whereas the pro-tostar can react almost immediately to a high accretionrate and quickly expands, it can contract only on a KHtimescale even if the accretion drops to zero (e.g., Sakuraiet al. 2015). Thus, even short episodes of accretion burstscan have a long-lasting impact on the star’s evolution. Asimilar feature is also seen for the second burst eventat t (cid:39) × years. This effect can only be correctlysimulated by solving for the stellar interior structure si-multaneously with time-dependent accretion histories, asprovided by SE-RHD simulations.In our simulations, circumstellar disks normally formperpendicular to the polar axis of the spherical coordi-nate system because of our choice of the polar axis ori-entation (see Section 3.1). The orientation of the diskoccasionally fluctuates with time, however, especially forcases A and B, for which the spin parameters measuredat the cloud scale are relatively low (Fig. 1). For in-stance, Figure 16 shows the density contours of the diskand ionization contours of the bipolar H ii region at differ-ent epochs for case A. Although viewing from the sameangle, the orientations of the H ii region as well as thedisk changes with time. Since the primordial cloud formsthrough cosmological structure formation and is not nec-essarily perfectly axisymmetric, material falling onto thedisk has different mean angular momentum vectors atdifferent times. Similar phenomena have been also re-ported for present-day star formation, whereby turbulentmotions in molecular clouds easily twist the distributionof the angular momentum (e.g., Bate et al. 2010; Field-ing et al. 2015). Due to this disk “precession”, the stellarEUV radiation can ionize the envelope gas over ratherwide opening angles. The orientation of the disk at anyparticular time also depends on whether or not UV feed-back is included in the simulations. This is not surprisingbecause UV feedback blows away the gas in the polar re-gions and modifies the pressure structure even near theequator. As a result, UV feedback can modify the accre-tion history of both angular momentum and mass fromthe envelope onto the disk. Case C: — For our case C, which also shows the inter-mittent UV feedback in an early stage, we follow themuch longer evolution than the above. For cases A andB, mass accretion is still on-going under intermittent UVfeedback for 7 × years, because episodes of rapid massaccretion are able to interrupt the KH contraction suf-ficiently often. Although it is uncertain how long thisevolution can continue for these cases, UV feedback willlikely fix the final stellar masses, if KH contraction re-mains uninterrupted for a sufficiently long time. Our caseC exemplifies this evolution. After the protostar entersthe supergiant stage for M ∗ (cid:38) M (cid:12) , the star againbegins KH contraction, after the stellar mass reaches ∼ M (cid:12) at t (cid:39) years (Fig. 6). The bipolar H ii re-gion quickly expands as the star approaches the ZAMS.Mass accretion is not immediately shut off by UV feed-back, however, but continues at a much reduced rate. Asa result, the star’s mass almost doubles by t (cid:39) . × years, when the accretion rate falls below 10 − M (cid:12) yr − . Comparison to 2D SE-RHD Results
As mentioned in Section 3.1, HR14 have previouslystudied the evolution of protostellar accretion for theassive Primordial Stars 17 time [ 10 years ] s t e ll a r m a ss : M ∗ [ M fl ] ABCDE
Fig. 17.—
Stellar mass growth histories obtained in the current3D SE-RHD simulations compared to our previous 2D results (Hi-rano et al. 2014). The solid lines represent the 3D cases with thesame colors as in Figure 2-(a); dashed lines depict the 2D resultsfor the same clouds. R ∗ [ R fl ] (a) stellar radius ABD stellar mass: M ∗ [ M fl ]10 S E UV [ s ec − ] (b) EUV emissivity Fig. 18.—
The protostellar evolution of radius and UV emis-sivities for the 3D cases considered here in comparison to those ofour previous 2D calculations (Hirano et al. 2014). This illustrationis almost the same as Figure 6; the panels (a) and (b) show theevolution of the stellar radius and EUV emissivity with increas-ing mass. For clarity only cases A, B, and D are presented. Thethick dashed lines represent the 2D results for the same clouds asin Figure 17. same primordial clouds considered above, but under theassumption of 2D axial symmetry. Figure 17 presentsthe stellar mass growth histories for both the current 3Dand the previous HR14 2D results in comparison. We seethat the 3D results roughly match the overall trends ofthe 2D results, in spite of a number of technical differ-ences between the adopted numerical codes. Comparingthe 2D and 3D results for each case, the stellar masses ata given epoch differ by a factor of a few at most. Somedifferences come from the fact that, for a 2D simulation,the original 3D structure of a natal cloud is modified tofit the 2D axisymmetric grid by averaging physical quan-tities over the azimuthal direction. Moreover, the 2D s t e ll a r m a ss : M ∗ ( M fl ) (a) PISNAKCF -5 -4 -3 -2 -1 infall rates at cloud scale: ˙ M J ( M fl yr − )10 -2 -1 p r o b a b ili t y (b) Fig. 19.—
Panel (a) : correlation between the final stellar massesand infall rates measured at the cloud (Jeans) scale when the cloudcentral density is (cid:39) cm − . The rhombuses and triangles de-note the current 3D results for cases A (green), B (blue), C (black),D (red), and E (magenta). The triangles are lower limits to themass, indicating that mass accretion continues beyond t = 7 × years, the end of the simulations. The cyan dots represent the 2DSE-RHD results for about 110 cosmological cases studied by Hi-rano et al. (2014). The dashed line represents the fitting formulafor these 110 cases given by equation 24. The horizontal yellow bardenotes the expected mass range of progenitors of pair-instabilitysupernovae, 120 M (cid:12) (cid:46) M ∗ (cid:46) M (cid:12) (Yoon et al. 2012). Thearrows mark the estimated stellar masses of Pop III supernovaprogenitors, which have been estimated from the abundance pat-terns of recently-discovered Galactic metal deficit stars: (A) SDSSJ001820.5-093939.2 (Aoki et al. 2014), (K) SMSS J031300.36-670839.3 (Keller et al. 2014), (C) SDSS J102915+172927 (Caffauet al. 2011), and (F) SD 1313-0019 (Frebel et al. 2015). The verticalbar associated with each arrow represents the possible mass rangediscussed in the literature (e.g., Schneider et al. 2012; Ishigaki et al.2014). Panel (b) : probability distribution of the cloud-scale infallrates in the 1540 cosmological cases studied by Hirano et al. (2015)(also see their Fig. 17). results also depend on the somewhat ad-hoc use of an α -viscosity to model the gravitational torque operatingin self-gravitating disks.Figure 18 shows the protostellar mass and UV emissivi-ties as a function of stellar mass for several representativecases. For case D, the evolution is quite similar between3D and 2D cases. The qualitative differences appear forcases A and B, where both the stellar radius and EUVemissivity abruptly change for M (cid:12) (cid:38) M (cid:12) because ofthe variable accretion encountered in 3D (also see Sec-tion 4.2.2). In our previous 2D cases, by contrast, therewas an intermediate stage between the supergiant andZAMS stage, where the stellar radius gradually increaseswith the mass (named the “P2” path in HR14). Stellarevolution calculations demonstrate that this behavior isgenerally present, when the accretion rates are in therange of 4 × − M (cid:12) yr − (cid:46) ˙ M ∗ (cid:46) − M (cid:12) yr − (e.g.,Omukai & Palla 2003). For 2D cases, this stage appearsbecause the accretion rates smoothly evolved with time.In 3D, however, the highly variable accretion rates donot remain in this narrow range very long. The inter-mediate P2 stage thus does not occur in our 3D cases.8 T. Hosokawa et al.Instead, the protostar evolves back and forth betweenthe supergiant and ZAMS stages depending on detailsof the variable accretion history. This effect reduces thestellar EUV emissivity for M ∗ (cid:38) M (cid:12) in comparisonwith the 2D cases, as episodes of rapid accretion bringthe star back to the supergiant stage several times.Figure 19-(a) displays the stellar masses at the endof the 3D simulations as a function of the infall ratesmeasured at the cloud (Jeans) scale (see Section 3.1),demonstrating that the 3D results still show a similarcorrelation as found in 2D. Considering the fact that themasses for cases A and B are only lower limits, the corre-lation found in 3D appears to have a steeper slope thanfor the 2D results (equation 24), albeit with much poorerstatistical significance (five cases in 3D). This does notnecessarily agree with Susa et al. (2014), who do not findthis correlation in their 3D RHD simulations. We discussthis issue in Section 5.2 separately.As Figure 19-(b) shows, primordial clouds taken forcases A–C offer rather typical initial conditions for pri-mordial star formation. Our results suggest that the verymassive stars, even above the mass range of the pair-instability supernovae (PISN) progenitors 120 M (cid:12) (cid:46) M ∗ (cid:46) M (cid:12) (Yoon et al. 2012), can form for suchtypical cases. A significant number of ∼ M (cid:12) BHsmight be produced as remnants of these very massiveprimordial stars. We shall discuss possible observationalsignatures of such massive primordial stars in Section 5.4.
Varying Spatial Resolution
Disk Fragmentation and Accretion Bursts
Although our 3D simulations show a new 3D effectof the intermittent UV feedback, the ansatz describedin Section 2.1.2 can bring some resolution-dependencewhich can affect the evolution quantitatively. It is thusnecessary to explore how resolution affects the results ofour simulations. First, we shall focus on case B describedin Section 4.2.2, for which UV feedback becomes impor-tant when the stellar mass exceeds 200 M (cid:12) at t (cid:38) × years (Fig. 14). For this case, we study effects of increas-ing the grid resolution during early times before UV feed-back begins to operate.We first recalculate the evolution of case B with ourstandard resolution, turning off UV feedback (case B-NF; see Table 1). We use four model configurations ofcase B-NF at different epochs for M ∗ (cid:39)
0, 40, 70, and120 M (cid:12) as the starting configurations for calculations athigher spatial resolution. Next, we interpolate each ofthe four starting model configurations onto a grid withdouble the resolution for each spatial coordinate and fol-low the evolution for 3000 years. These test cases arenamed as B-NF-HR2-mX, where the suffix “mX” repre-sents the restarting point at M ∗ (cid:39) X M (cid:12) (see Table 1).The duration of 3000 years is comparable to the Keplerorbital period P Kepler (cid:39) M ∗ (cid:39) M (cid:12) and disk radius ∼
500 AU (c.f.equation 19). The accretion envelope surrounding thedisk typically has a density (cid:46) cm − , for which thefree-fall timescale is (cid:38) × years. Mass infall from theenvelope onto the disk at the higher resolution hardly af-fects the evolution, whereas cooling, fragmentation, andthe associated mass and angular momentum transportwithin the disk occur on much shorter timescales and therefore should be affected by varying spatial resolu-tion.Figure 20 compares the stellar mass growth and ac-cretion histories for the examined cases. We see thatmass accretion becomes more strongly variable withhigher resolution; i.e., the accretion rates vary overshorter timescales and with much higher variations:10 − M (cid:12) yr − (cid:46) ˙ M ∗ (cid:46) . M (cid:12) yr − (Fig. 20-b). Fur-thermore, the mass that the star accretes during the3 × years is somewhat larger with higher resolution,independent of the different starting points of the doubly-resolved cases (Fig. 20-a). This can be understood withthe help of Figure 3, which displays the face-on diskstructure at t (cid:39) . × years for the different reso-lutions. The disk clearly has sharper and finer densitystructure at higher resolution. This is because radiativecooling is less restricted by our ansatz and the disk isthus more gravitationally unstable. Since gravitationaltorques are caused by the non-axisymmetric disk struc-ture such as spiral arms, increasing the resolution canlead to more efficient mass and angular momentum trans-port through the disk, which results in a higher stellarmass.We once again double the resolution using case B-NF-HR2-m0 at t (cid:39)
600 years as a starting model, a time atwhich M ∗ (cid:39) M (cid:12) , and recalculate the evolution (caseB-NF-HR4-m20). We repeat these higher resolution sim-ulations with a larger cooling limit f limit = 24 (case B-NF-HR4-m20-fc; see Section 2.1.2) and a larger sink cellradius r min = 52 AU (case B-NF-HR4-m20-lsk). Theresolution is now 4 times higher than the original caseB-NF. Because of the high computational cost, the high-est resolution simulations are evolved for only 10 years.Figure 21 compares the stellar mass growth and accre-tion histories of these quadruple resolution cases withthe results of the lower resolution simulations. The res-olution dependences discussed above appear again moreprominently.Case B-NF-HR4-m20’s accretion history displays lotsof spiky features, representing multiple accretion burstevents (Fig. 21-b). The stellar mass growth over 10 years is correspondingly greater at higher resolution(Fig. 21-a). Figures 3-(c) and 22 show that, at this highresolution, one of the spiral arm breaks up and createstwo fragments. These and Figure 23-(a) show the dis-tributions of the Toomre Q -parameter on the equator,defined by Q = c s Ω πG Σ , (27)where c s is the local sound speed and Σ is the verticalcolumn density. Because the spiral arms have Q (cid:46) s t e ll a r m a ss : M ∗ [ M fl ] Case B (a) years ]10 -4 -3 -2 -1 a cc . r a t e s : ˙ M ∗ [ M fl y r − ] (b) Fig. 20.—
Comparison of the mass growth histories with differentspatial resolutions in halo B cases. The panels (a) and (b) showthe time evolution of the stellar masses and accretion rates ontothe protostar. In both panels, the thick dashed lines representcase B-NF. The solid lines represent the cases with 2 times higherspatial resolution: B-NF-HR2-m0 (red), B-NF-HR2-m40 (blue), B-NF-HR2-m70 (magenta), and B-NF-HR2-m120 (green). The filledcircles mark the starting points of these doubled-resolution cases,i.e., when the stellar mass is (cid:39) M (cid:12) , 40 M (cid:12) , 70 M (cid:12) , and 120 M (cid:12) for case B-NF (also see text). s t e ll a r m a ss : M ∗ [ M fl ] (a) Case B defaultlarger f_limit: 24larger sink: 52AU
600 800 1000 1200 1400 1600time [ years ]10 -3 -2 -1 a cc . r a t e s : ˙ M ∗ [ M fl y r − ] (b) Fig. 21.—
Same as Figure 20 but including cases at 4 times higherspatial resolution than case B-NF (thick dashed lines) and 2 timeshigher resolution than case B-NF-HR2-m0 (dotted lines). The solidlines represent cases (see text and Table 1) B-NF-HR4-m20 (black),B-NF-HR4-m20-fc (blue), and B-NF-HR4-m20-lsk (red). In panel(b), accretion histories for cases B-NF, B-NF-HR2-m0, and B-NF-HR4-m20 only are shown for clarity. !" $&$23(20 ! (cid:3) !4 $&$28(,0 ! (cid:3) )669: (cid:3) ;<7=>?$@*?/A%-$ (cid:3) B$'(*+2$ $ B$'(*+,$ $ B$'66$ $ B$'6$CD$;> BE F$ (cid:3) ' (cid:3) E (cid:3) Fig. 22.—
Inward migration of two fragments (crosses) for caseB-NF-HR4-m20. As in Figure 3, colors represent the gas columndensity but for a smaller region around the central star. The black(white) contours delineate the boundaries where the Toomre-Q pa-rameter takes the value Q = 1 . Q = 0 . Q < Panel (a):
Evolutionary age t = 1 . × years (the same moment as in Fig. 3) Panel (b): (c). Both fragments rapidly migrate inward and reachthe star in less than 100 years. Recent 3D simulationsby Greif et al. (2012) also report a similar rapid migra-tion in highly gravitationally unstable disks. Greif et al.(2012) conclude that the typical timescale for such migra-tions is roughly the local free-fall timescale, the shortesttimescale for any physical process driven by gravity. Forour case B-NF-HR4-m20, the density near the disk mid-plane is n (cid:38) cm − , for which the free-fall timescaleis t ff (cid:46) . × years. The migration timescale is there-fore also of the order of the local free-fall timescale inour calculations. The rapid inward migration of the frag-ments in gravitationally unstable disks is also commonlyseen in simulations of present-day star and planet forma-tion (e.g., Vorobyov & Basu 2006; Baruteau et al. 2011;Machida et al. 2011; Zhu et al. 2012). These studies showthat rapid migration occurs roughly over the so-calledtype I migration timescale t mig , I = h Cµq π Ω (28)(Tanaka et al. 2002), where C = 2 . . ξ , ξ is the0 T. Hosokawa et al. Fig. 23.—
Radial distributions of physical quantities for the cir-cumstellar disks shown in Figure 22: (a) Toomre Q -parameter, (b)vertical column density Σ, (c) fragment mass M f (equation 29),and (d) type I migration timescale (equation 28). In each panel,the red curves pertain to the evolutionary state depicted in Fig-ure 3-(c) and Figure 22-(a), whereas the blue curves correspond to80 years later (see Figure 22-(b)). These two epochs are also dis-played in Figure 24. In each panel, solid curves depict quantitiesaveraged over the azimuthal φ direction. The dashed lines in panel(a) represent the minimum values Q min ( r ) in an annulus of radius r , which roughly trace the Q values through either spiral arms orfragments. In panel (b), the thick dotted line shows the profileΣ ∝ r − for reference. The colored dashed lines in panels (c) and(d) represent the mininum values of the estimated fragment massand corresponding migration timescale in the annulus of radius r .In panel (c), the red filled circles indicate the masses of the twofragments marked in Figure 22-(a). In panel (d), we also plot theKepler orbital period P Kepler = 2 π (cid:112) r /GM ∗ for M ∗ = 45 . M (cid:12) with a black dashed line. !" $&$23(20 ! (cid:3) !4 $&$28(,0 ! (cid:3) )669: (cid:3) ;."<=*>%$="//$ (cid:3) ?$,6$ $ ?$,$ $ ?$6(,$ $ ?$6(6,$@0 ! A$ (cid:3) ' (cid:3) B (cid:3) Fig. 24.—
Same as Figure 22 but for the fragment masses eval-uated locally with equation (29). In panel (a), we only mark thepositions of the two fragments referred to in Figures 22 and 23 forthis epoch. radial gradient of the column density Σ ∝ r − ξ , h is thedisk aspect ratio, µ = π Σ r /M ∗ , and q = M f /M ∗ , where M f is the fragment mass. Figure 23-(b) shows that thecolumn density profile in our disk is well approximatedby Σ ∝ r − , i.e., ξ = 1. The fragment mass M f isapproximated by the mass enclosed within a spatial scaleof the most unstable wavelength λ = 2 c s /G Σ, M f = π ( λ/ Σ = πc s G Σ . (29)Figure 23-(c) shows that the φ -averaged fragment massesestimated by equation (29) span the range from ∼ M (cid:12) to several 10 M (cid:12) . Because there are large density fluctu-ations in the disk, the locally-estimated fragment massvaries considerably. Figure 24 shows M f ∼ . − M (cid:12) along the spiral arms, where the fragmentation often oc-curs. We also directly calculate the masses of the twofragments seen in Figure 22-(a) (marked as 1 and 2),by summing up the masses in the cells with the columndensity higher than 5000 g cm − around each fragment.The resulting values, M f (cid:39) . M (cid:12) and M f (cid:39) . M (cid:12) ,are in rough agreement with the estimated values in Fig-ure 24 (see also Figure 23-(c)).In Figure 23-(d), we plot the type I migration timescalegiven by equation (28). Considering the fragment massesdiscussed above, we estimate that the corresponding mi-assive Primordial Stars 21gration timescales should be around ∼
100 years, whichare shorter than or comparable to the Kepler orbital peri-ods. The rapid inward migration seen in our simulationscan be well explained by the above evaluation. In addi-tion to the two infalling fragments shown in Figures 3-(c)and 22, several similar events were observed during thesimulation duration of 10 years. Some of the accretionbursts seen in Figure 21, although not all, are caused bysuch events.As mentioned earlier, we have also performed quadru-ple resolution simulations with different numerical pa-rameters. For case B-NF-HR4-m20-fc, for instance, weadopt f limit = 24 instead of the fiducial value f limit =12, where f limit sets the cooling limit threshold (equa-tion 10). As described in Section 2.1.2, a larger value for f limit sets a more stringent limit on cooling. Figure 21shows that the evolution for this case is quite similar tothat for the double resolution case B-NF-HR2-m0, ver-ifying that the resolution dependences examined abovereally come from the more efficient cooling allowed withhigher resolution. For case B-NF-HR4-m20-lsk, on theother hand, we adopt a larger sink r min = 52 AU in-stead of the fiducial value r min = 30 AU. The stellarmass growth history for this case roughly traces that forcase B-NF-HR4-m20, although with less variable accre-tion rates. The shortest timescale of accretion variabilityis roughly given by the Kepler orbital timescale near theradial inner boundary at r = r min . Stellar UV Feedback with Variable Accretion
Next we investigate the effects of varying the grid res-olution for evolutionary stages when UV feedback be-comes important. For this purpose, we take case C (seeSection 4.2.2) as the basis for comparison. Unlike caseB considered above, UV feedback operates earlier and atlower stellar masses for case C (e.g., Figs. 4 and 5). Thisenables us to follow the evolution with UV feedback fora longer period of time at a reasonable computationalcost. We double the resolution for the initial configura-tion, M ∗ = 0 M (cid:12) , and follow the evolution for 5 . × years (case C-HR2).Figure 25 displays the stellar mass (a) accretion rate(b), stellar radius (c), and UV emissivity (d) as a functionof time for both case C and C-HR2. In agreement withwhat we have shown in Section 4.3.1, the accretion rateis more highly variable at the higher resolution (Fig. 25-b). After the different accretion histories over 5 . × years, however, the final masses differ only by ∼
10 %.This is because the current evolutionary timescale ismuch longer than the local Kepler timescale within thedisk (equation 19), so that we are no longer consideringthe short term evolution of a nearly isolated disk; themass supply from the envelope ultimately controls stel-lar growth via mass accretion through the disk. Since thestellar radius abruptly changes with the variable mass ac-cretion, the UV emissivity also displays very large fluc-tuations (Fig. 25-c, d). With the doubled resolution, inparticular, the protostellar bloating and subsequent KHcontraction are repeated several times because of multi-ple accretion bursts. The extinction and re-formation ofbipolar H ii regions also occurs accordingly. In compari-son with the original case C, UV feedback begins to limitstellar mass growth earlier at M ∗ (cid:39) M (cid:12) , but the starnevertheless accretes a larger amount of gas ( (cid:39) M (cid:12) ) Fig. 25.—
Effects of doubling the spatial resolution for a long-term evolution under the influence of UV feedback for case C. Theevolution of (a) stellar mass, (b) accretion rate, (c) stellar radius,and (d) emissivity of ionizing photons is presented. In each panel,the black and green lines represent the fiducial case C and a testcase with 2 times higher spatial resolution, C-HR2 (also see Ta-ble 1). after the onset of significant UV feedback. The overallresult is a slightly higher stellar mass at the end of thesimulation than for the original case C. DISCUSSION
Impact of Dissociation Feedback
As described in Section 4.2.1, the accretion limitingeffect of stellar dissociating (FUV) radiation is modestin our simulations; FUV feedback slightly reduces massaccretion onto the star, but does not completely stopit, as opposed to ionizing (EUV) feedback (Figs. 8 and9). This is in contrast to the conclusions of Susa (2013)and Susa et al. (2014), who show that mass accretioncan be halted by FUV feedback alone. EUV feedback isnot included in their SPH/N-body simulations because of2 T. Hosokawa et al.the limited spatial resolution, especially in polar regionswhere an H ii region should first appear. To examine theroles of FUV feedback for our simulations in comparisonwith theirs, we return to our case D-DF, for which EUVfeedback had been artificially turned off.As described in Susa (2013), a circumstellar disk ofdensity n (cid:38) cm − is ultimately dispersed by FUVfeedback in their simulations. This is because H forma-tion heating mainly due to the three-body reaction effi-ciently heats the gas when the dense gas ( n (cid:38) cm − )is exposed to FUV radiation. As the disk disperses, thegas density around stars rapidly falls to n (cid:46) cm − ,resulting in low accretion rates onto the stars. For ourcase D-DF, however, FUV radiation heats up the gas pri-marily in the polar regions throughout the course of 8 × years evolution. The dense disk with n (cid:38) cm − survives in our case D-DF simulation, being protectedagainst the FUV radiation by efficient H self-shieldingin the inner regions of the disk. The pressure excess inthe polar PDR is much less than what an H ii region canprovide when stellar EUV radiation is present. Unlikethe effects of EUV feedback, the pressure structure inthe envelope is hardly modified outside the PDR. There-fore mass accretion onto the disk continues from the en-velope as well as from the disk onto the star with littlemodification by FUV feedback.It seems that the disagreement outlined above stemsfrom different efficiencies of H self-shielding in the vicin-ity of the star. Both simulations suffer from the inabil-ity to exactly model this effect in the innermost regions r (cid:46)
30 AU, currently masked by the sink. Numerical sim-ulations devoted to resolving small-scale structure nearthe stellar surface will be needed in the future.
Formation of Binary and Multiple Stellar Systems
Our current simulations do not show clear indicationsfor the formation of binary or multiple stellar systems.As explained in Section 4.3, disks will readily fragment,especially in simulations with high spatial resolution,but fragmentation actually enhances angular momentumtransport in the disk and mass accretion onto the centralstars via tidal torques and the rapid inward migration ofthe fragments. This is in disagreement with recent simu-lations claiming the formation of multiple stellar systems(see Section 1).First of all, our inability to form binary and multiplestellar systems may be partly due to a selection bias;we did not choose primordial clouds that have stronglyasymmetric or complex structure as starting configura-tions. Among HR14’s cosmological sample, for instance,some clouds have multiple density peaks emerging beforethe formation of embryo protostars (see also Turk et al.2009). With typical separations of 0 . − (cid:39) (cid:46) AU) bi-nary or multiple stellar systems could form via disk frag-mentation. Whereas a number of the fragments migrateinward to fall onto the star, some could potentially sur-vive orbiting around the star or even as ejecta via gravita-tional interactions with other fragments. Our numericalmethod using the spherical coordinates and a stationarycentral sink does allow the formation of such systems. In fact, Vorobyov & Basu (2015) report finding such eventsin their simulations using the similar numerical settings.We cannot discount the potential claim that, had we fol-lowed the long-term evolution at sufficiently high spatialresolution, we might also observe such events in our cases(e.g., Stacy & Bromm 2013). However, Vorobyov & Basu(2015) also point out that the ejection events are muchrarer than the burst accretion events initiated by the mi-gration of fragments. If the fragments grow to becomemassive stars emitting a copious amount of UV photons,then our current framework will have the limitation ofnot dealing with the multiple light sources.It is also possible that binary or multiple stellar sys-tems could form within our 30 AU sink. For present-daymassive star formation, close binaries are the rule ratherthan the exception (e.g., Zinnecker & Yorke 2007). Forexample, the brightest Trapezium star, Ori Θ C has com-ponents “1” and “2” separated by ∼
10 AU. One of the bi-nary components of the second brightest star, Ori Θ A1,is itself a spectroscopic binary with a 1 AU separation.Ori Θ E and Ori Θ B1 are also spectroscopic binarieswith separations ∼ . ∼ .
13 AU. Ori Θ D is abinary with a ∼
10 AU separation. All in all, the Trapez-ium has at least 14 close components in 5 multiple stellarsystems. If such systems are also common among primor-dial stars, close binary BH systems could form after acommon envelope phase and eventually merge emittinggravitational waves. Interestingly, the recent detectionGW150914 reported by Abbott et al. (2016) may reallyrepresent the signature of primordial close binaries (Kin-ugawa et al. 2014).In general, the formation efficiency of binary or mul-tiple systems depends on the production and survivalrates of the fragments and perhaps on fission of rapidlyrotating protostars, of which we know very little. Asdemonstrated in Section 4.3, the fragmentation ratesin our simulations clearly depend on the spatial resolu-tion used to follow the long-term evolution. Machida &Doi (2013) investigate the spatial resolution necessary tocapture disk fragmentation with numerical convergence,meaning that the results no longer depend on the resolu-tion. They conclude that numerical codes need to resolvestructure at least down to the scale ∼ . − . ∼
10 years could be followed. Note that eventhese studies may suffer from limited spatial resolution;the numbers of cells assigned per Jeans length are 8 and32 in Machida & Doi (2013) and Greif et al. (2012), re-spectively, whereas it is known that significantly higherresolution is required to also capture the turbulent mo-tions expected in the early collapse stage (e.g., Turk et al.2012).We conclude that with current computational capabil-ities it is necessary to use sink cells or sink particles tofollow the long term evolution up to the point that UVfeedback becomes important. For such cases, disk frag-mentation and the resulting stellar multiplicity will nec-essarily depend on the adopted spatial resolution. Stel-lar multiplicity should affect the evolution. Susa (2013)allow the creation of multiple sink (star) particles with-assive Primordial Stars 23out further mergers, and show that UV feedback frommultiple protostars jointly halt the mass accretion. Thisfundamentally different implementation of sinks may ex-plain why Susa et al. (2014) find no correlation betweenthe final stellar masses and cloud-scale infall rates, incontrast to our results (see Section 4.2.3).The evolution of the stellar radius is important for de-termining merger rates with other stars or fragments aswell as the associated UV emissivities. Once the proto-star inflates entering the supergiant stage after an accre-tion burst event, the large stellar size will further enhancethe merger rates of nearby fragments and close compan-ions. If a supergiant protostar is in a close binary sys-tem, for example, mass exchange with the companionstar could cause the companion to inflate and initiatecommon envelope evolution of the primary and compan-ion stars. As a protostar inflates, it can consume thegas in the inner parts of the circumstellar accretion disk,enhancing the instantaneous accretion rate even more.These processes can help the formation of very massivestars, in addition to the intermittent UV feedback shownin Section 4.2.2.
Effects of Magnetic Fields
We have not included effects of magnetic fields in ourcurrent simulations. However, theoretical studies predictthat the field strength can be significantly amplified bythe turbulent dynamo process during early collapse (e.g.,Schleicher et al. 2010; Turk et al. 2012). Magnetic fieldscould affect the evolution in the late accretion stage ina number of ways. The turbulent motion in accretiondisks can further increase the field strength and add tothe magnetic pressure support of the disk (e.g., Tan &Blackman 2005). The angular momentum of the accret-ing material can be reduced by magnetic breaking (e.g.,Machida & Doi 2013). Both effects could potentiallysuppress disk fragmentation. Nevertheless, some accre-tion variability should remain as long as the gravitationaltorque controls the mass and angular momentum trans-fer in disks.Magnetically-driven outflows might appear intermit-tently with the variable mass accretion (e.g., Machidaet al. 2011; Machida & Hosokawa 2013). It is interest-ing to speculate on the potential interaction between thestellar UV and outflow feedback mechanisms. Whereasthe UV feedback is turned off by burst accretion (Sec-tion 4.2.2), the outflow power will be enhanced becausea fraction of the accreting gas is generally diverted intothe outflow (e.g., Matzner & McKee 1999). The UV andoutflow feedback mechanisms may actually complementeach other; the outflow is relatively weak as an H ii regionexpands, and the outflow is activated as the H ii regiondisappears. Future numerical simulations such as thoseperformed for the present-day high-mass star formation(e.g., Kuiper et al. 2015) will explore these effects. Signatures of Primordial Stars in GalacticMetal-poor Stars
Although only a handful cases were studied above, itwould be useful to confront our results with current ob-servational constraints on the masses of primordial stars.Figure 19-(a) includes the recently estimated masses ofPop III supernova progenitors derived from abundance patterns in Galactic metal-poor stars. Caffau et al.(2011), Keller et al. (2014), and Frebel et al. (2015) haveshown that the abundance patterns of extremely metal-poor ([Fe/H] < −
4) stars can be explained by chemicalyields of core collapse supernovae (CCSN), whose pro-genitors are several × M (cid:12) stars. Among our simu-lations, such ordinary massive stars ultimately form forcases D and E, after the UV feedback efficiently shuts offmass accretion (Section 4.2.1).Aoki et al. (2014) recently found a peculiar abundancepattern in SDSS J001820.5-093939.2, which does notmatch the yields of CCSN but seems to indicate the ex-istence of more massive populations: progenitors of pair-instability supernovae (PISN; 120 M (cid:12) (cid:46) M ∗ (cid:46) M (cid:12) ,e.g., Yoon et al. 2012) or the core-collapse of very massivestars ( M ∗ (cid:38) M (cid:12) , e.g., Ohkubo et al. 2009). For ourcases A–C, the very massive stars above the PISN rangeform due to the process of intermittent (burst) accre-tion and its associated variability of UV feedback (Sec-tion 4.2.2). Although currently absent among the fivecases considered here, we strongly suspect that we wouldalso have formed stars in the PISN range, if a greaternumber of cases had been considered. Figure 19 sug-gests that such stars will appear for the cases that havethe cloud-scale infall rates between the values for casesA and C. However, Aoki et al. (2014) actually show thatsuch very massive stars may be preferred to the PISNprogenitors to explain the peculiar abundance pattern.Note that this object has a higher metallicity [Fe/H] (cid:39) − . M ∗ (cid:38) M (cid:12) primordial starsconstitute a majority (Fig. 19), so it might be possible tofind similar signatures in other nearby stars. However,this still may be challenging, because such energetic su-pernovae, could potentially expel most of the gas out ofa halo and suppress subsequent star formation (Cooke& Madau 2014). It will be necessary to determine thestar formation efficiency after energetic supernovae in or-der to assess the mass distribution of massive primordialstars using stellar archaeology. Formation of Supermassive Stars: an ExtremeCase
Several authors have considered the formation of pri-mordial supermassive ( M ∗ ∼ M (cid:12) ) stars (SMSs) asa promising pathway to seed supermassive BHs in theearly universe (e.g., Bromm & Loeb 2003). A popularmodel supposes that such an extreme mode of primor-dial star formation should occur in the so-called atomic-cooling haloes exposed to strong FUV radiation. Withonly atomic hydrogen as an available coolant, the cloudcollapse advances almost isothermally at a higher tem-perature ( T (cid:39) ∼ Myr, duringwhich the stellar mass could attain values (cid:38) M (cid:12) (e.g., Regan & Haehnelt 2009; Latif et al. 2013; Inayoshiet al. 2014). Recent numerical simulations show that,also for such cases, circumstellar disks easily fragmentdue to the gravitational instability (e.g., Regan et al.4 T. Hosokawa et al.2014; Becerra et al. 2015). It is also expected that theemerging fragments will rapidly migrate inward throughthe disk (e.g., Inayoshi & Haiman 2014), as for normalPop III cases (see Section 4.3.1). The resulting accretionrates will thus show significant variability.As discussed above, stellar UV feedback will be veryweak whenever the protostar enters the supergiant stagefor accretion rates exceeding 0 . M (cid:12) yr − . However, theprotostar will contract toward the ZAMS under realisticconditions of variable accretion, each time the accretionrate falls below 0 . M (cid:12) yr − . Stellar evolution calcu-lations show that such a slowly accreting star can reachthe ZAMS, if this lower rate of accretion continues for (cid:38) years (Sakurai et al. 2015). As it approaches theZAMS, the stellar UV emissivity rises sharply, initiat-ing very powerful UV feedback (Section 4.2.2). For thenormal primordial cases studied above, a short burst ofrapid mass accretion will cause the star to bloat up againand the star can again accrete more available material.Although case B has reached nearly 600 M (cid:12) through thisrecurrent accretion process and is still accreting materialafter 7 × years, our simulations are still far from themass 10 M (cid:12) and average accretion rates 0 . M (cid:12) yr − necessary to produce such an SMS within a Myr. Ournewly developed SE-RHD numerical code will be usefulto explore the interplay between variable accretion andUV feedback for the formation of SMSs.Needless to say, other physical processes which havenot been fully considered, such as stellar multiplicity(Section 5.2) and magnetic fields (Section 5.3), can alsojointly modify the protostellar accretion histories. Sincethe formation of SMSs is an extreme version of the forma-tion of very massive stars studied in Section 4.2.2, bothcases share many common theoretical issues to be solvedin future studies. CONCLUSIONS
We have studied the long term evolution during theaccretion stage of primordial star formation, focusing onthe interplay between highly variable accretion onto thestar and stellar UV feedback. To follow the evolutionby numerical simulations, we have developed a new 3DSE-RHD hybrid code to solve for the radiation hydro-dynamics of the accretion flow simultaneously with thestellar evolution of the accreting protostar. The evolu-tion of the protostellar accretion flow under the influenceof both ionizing (EUV) and photodissociating (FUV) ra-diative feedback has been followed for over ∼ years.We have examined the evolution of five different primor-dial clouds occurring in representative dark matter halosfound in cosmological structure formation simulations.Many numerical parameters (e.g., spatial resolution, sinkcell size) used in our simulations have been varied to as-sess their effects on the results. Our findings are summa-rized as follows:1. The resulting stellar masses M ∗ show a great di-versity: 10 M (cid:12) (cid:46) M ∗ (cid:46) M (cid:12) , which alsohas been inferred in previous studies (e.g., HR14,HR15; Susa et al. 2014). The stellar mass tendsto be higher for clouds with a higher average in-fall rate, a cloud parameter set during the earlycollapse stage before the formation of an embryoprotostar. This roughly agrees with our previous findings discussed in HR14.2. Ordinary massive ( M ∗ (cid:46) M (cid:12) ) stars form incases with relatively slow mass accretion. For thesecases, UV feedback efficiently shuts off accretiononce a protostar reaches the late KH contractionstage, as shown in our previous 2D studies. On theother hand, very massive ( M ∗ (cid:38) M (cid:12) ) starsalso form when relatively rapid accretion is avail-able. For such cases, the protostar goes througha “supergiant protostar” evolutionary stage, whichrecurrently appears with the highly variable massaccretion associated with self-gravitating circum-stellar disks. As a result of the fluctuating stel-lar radius, characterized by alternating periods ofcontraction and inflation, the formation and extinc-tion of H ii regions also recur repeatedly. Over thecourse of time UV feedback only operates intermit-tently, and does not efficiently halt mass accretiononto the star.3. Very massive stars form in three out of the fiveexamined cases (and two of these three yield onlylower limits to the final stellar masses). Accordingto a statistical study by HR15, who have analyzedphysical properties of more than 1500 primordialclouds, the formation of such very massive starscan be a major mode of primordial star formation.4. According to our systematic assessment of varyingthe numerical parameters, our current simulationsdo not converge with increasing spatial resolution.As far as we can tell with our present suite of nu-merical tools, however, variable accretion and theresulting intermittent UV feedback becomes evenmore prominent at higher spatial resolution.Our results show that the interplay between stellar UVfeedback and variable mass accretion can help the forma-tion of very massive primordial stars, making this a likelypathway to form a significant number of very massivestars in the early universe, circumventing potential barri-ers, such as stellar UV feedback and disk fragmentation.The peculiar abundance pattern recently discovered ina Galactic metal-poor star (Aoki et al. 2014) could bethe observational signature indicative of such very mas-sive stars. With the resolution dependence of our resultsdiscussed in Section 4.3, this basic proof of concept in-vestigation still requires further detailed study to verifythe above idea.As an extension of cases considered here, the forma-tion of even more massive (or supermassive) stars maybe possible under certain favorable conditions; the feasi-bility of this scenario will be examined in future studies.Such studies will be necessary to determine the maximummass of primordial stars, a key quantity when consider-ing the origin of the supermassive black holes observedin the early universe.The authors thank Hajime Susa, Eduard Vorobyov,Masayuki Umemura, Shu-ichiro Inutsuka, and KenNomoto for fruitful discussions and comments. The nu-merical simulations were performed on the Cray XC30at the Center for Computational Astrophysics, CfCA, ofassive Primordial Stars 25the National Astronomical Observatory of Japan. Por-tions of this work were conducted at the Jet PropulsionLaboratory, California Institute of Technology, operat-ing under a contract with the National Aeronautics andSpace Administration (NASA). This work was financiallysupported by the Grants-in-Aid for Basic Research bythe Ministry of Education, Science and Culture of Japan(25800102, 15H00776: TH, 25287040: KO, 25287050:NY) and by Grant-in-Aid for JSPS Fellows (SH). RK ac-knowledges funding within the Emmy Noether ResearchGroup on “Accretion Flows and Feedback in RealisticModels of Massive Star Formation” granted by the Ger-man Research Foundation (DFG) under grant no. KU2849/3-1.6 T. Hosokawa et al. REFERENCESAbbott, B. P., et al. 2016, Physical Review Letters, 116, 061102Abel, T., Anninos, P., Zhang, Y., & Norman, M. L. 1997, NewAstronomy, 2, 181Abel, T., Bryan, G. L., & Norman, M. L. 2002, Science, 295, 93Aoki, W., Tominaga, N., Beers, T. C., Honda, S., & Lee, Y. S.2014, Science, 345, 912Baruteau, C., Meru, F., & Paardekooper, S.-J. 2011, MNRAS,416, 1971Bate, M. R., Lodato, G., & Pringle, J. E. 2010, MNRAS, 401,1505Becerra, F., Greif, T. H., Springel, V., & Hernquist, L. E. 2015,MNRAS, 446, 2380Bromm, V., Coppi, P. S., & Larson, R. B. 2002, ApJ, 564, 23Bromm, V., & Loeb, A. 2003, ApJ, 596, 34Bromm, V., Yoshida, N., Hernquist, L., & McKee, C. F. 2009,Nature, 459, 49Caffau, E., et al. 2011, Nature, 477, 67Clark, P. C., Glover, S. C. O., & Klessen, R. S. 2012, MNRAS,420, 745Clark, P. C., Glover, S. C. O., Smith, R. J., Greif, T. H., Klessen,R. S., & Bromm, V. 2011, Science, 331, 1040Cooke, R. J., & Madau, P. 2014, ApJ, 791, 116DeSouza, A. L., & Basu, S. 2015, MNRAS, 450, 295Draine, B. T., & Bertoldi, F. 1996, ApJ, 468, 269Federrath, C., Sur, S., Schleicher, D. R. G., Banerjee, R., &Klessen, R. S. 2011, ApJ, 731, 62Fielding, D. B., McKee, C. F., Socrates, A., Cunningham, A. J.,& Klein, R. I. 2015, MNRAS, 450, 3306Forrey, R. C. 2013, ApJ, 773, L25Frebel, A., Chiti, A., Ji, A. P., Jacobson, H. R., & Placco, V. M.2015, ApJ, 810, L27Galli, D., & Palla, F. 1998, A&A, 335, 403Greif, T. H. 2014, MNRAS, 444, 1566—. 2015, Computational Astrophysics and Cosmology, 2, 3Greif, T. H., Bromm, V., Clark, P. C., Glover, S. C. O., Smith,R. J., Klessen, R. S., Yoshida, N., & Springel, V. 2012,MNRAS, 424, 399Greif, T. H., Springel, V., White, S. D. M., Glover, S. C. O.,Clark, P. C., Smith, R. J., Klessen, R. S., & Bromm, V. 2011,ApJ, 737, 75Hartwig, T., Clark, P. C., Glover, S. C. O., Klessen, R. S., &Sasaki, M. 2015, ApJ, 799, 114Hirano, S., Hosokawa, T., Yoshida, N., Omukai, K., & Yorke,H. W. 2015, MNRAS, 448, 568Hirano, S., Hosokawa, T., Yoshida, N., Umeda, H., Omukai, K.,Chiaki, G., & Yorke, H. W. 2014, ApJ, 781, 60Hirano, S., & Yoshida, N. 2013, ApJ, 763, 52Hollenbach, D., Johnstone, D., Lizano, S., & Shu, F. 1994, ApJ,428, 654Hosokawa, T., & Omukai, K. 2009, ApJ, 691, 823Hosokawa, T., Omukai, K., & Yorke, H. W. 2012a, ApJ, 756, 93Hosokawa, T., Omukai, K., Yoshida, N., & Yorke, H. W. 2011,Science, 334, 1250Hosokawa, T., Yorke, H. W., Inayoshi, K., Omukai, K., &Yoshida, N. 2013, ApJ, 778, 178Hosokawa, T., Yorke, H. W., & Omukai, K. 2010, ApJ, 721, 478Hosokawa, T., Yoshida, N., Omukai, K., & Yorke, H. W. 2012b,ApJ, 760, L37Inayoshi, K., & Haiman, Z. 2014, MNRAS, 445, 1549Inayoshi, K., Omukai, K., & Tasker, E. 2014, MNRAS, 445, L109Ishigaki, M. N., Tominaga, N., Kobayashi, C., & Nomoto, K.2014, ApJ, 792, L32John, T. L. 1988, A&A, 193, 189Karlsson, T., Johnson, J. L., & Bromm, V. 2008, ApJ, 679, 6Keller, S. C., et al. 2014, Nature, 506, 463Kinugawa, T., Inayoshi, K., Hotokezaka, K., Nakauchi, D., &Nakamura, T. 2014, MNRAS, 442, 2963Kreckel, H., Bruhns, H., ˇC´ıˇzek, M., Glover, S. C. O., Miller,K. A., Urbain, X., & Savin, D. W. 2010, Science, 329, 69Kuiper, R., Klahr, H., Beuther, H., & Henning, T. 2010a, ApJ,722, 1556—. 2011, ApJ, 732, 20Kuiper, R., Klahr, H., Dullemond, C., Kley, W., & Henning, T.2010b, A&A, 511, A81Kuiper, R., & Klessen, R. S. 2013, A&A, 555, A7 Kuiper, R., & Yorke, H. W. 2013, ApJ, 772, 61Kuiper, R., Yorke, H. W., & Turner, N. J. 2015, ApJ, 800, 86Latif, M. A., Schleicher, D. R. G., Schmidt, W., & Niemeyer, J.2013, MNRAS, 433, 1607Levermore, C. D., & Pomraning, G. C. 1981, ApJ, 248, 321Machida, M. N., & Doi, K. 2013, MNRAS, 435, 3283Machida, M. N., & Hosokawa, T. 2013, MNRAS, 431, 1719Machida, M. N., Inutsuka, S.-i., & Matsumoto, T. 2011, ApJ,729, 42Martin, P. G., Schwarz, D. H., & Mandy, M. E. 1996, ApJ, 461,265Matzner, C. D., & McKee, C. F. 1999, ApJ, 510, 379McKee, C. F., & Tan, J. C. 2008, ApJ, 681, 771Meece, G. R., Smith, B. D., & O’Shea, B. W. 2014, ApJ, 783, 75Mignone, A., Bodo, G., Massaglia, S., Matsakos, T., Tesileanu,O., Zanni, C., & Ferrari, A. 2007, ApJS, 170, 228Millar, T. J. 1991, A&A, 242, 241Oh, S. P., & Haiman, Z. 2002, ApJ, 569, 558Ohkubo, T., Nomoto, K., Umeda, H., Yoshida, N., & Tsuruta, S.2009, ApJ, 706, 1184Omukai, K. 2001, ApJ, 546, 635Omukai, K., & Inutsuka, S.-i. 2002, MNRAS, 332, 59Omukai, K., & Nishi, R. 1998, ApJ, 508, 141Omukai, K., & Palla, F. 2003, ApJ, 589, 677O’Shea, B. W., McKee, C. F., Heger, A., & Abel, T. 2008, inAmerican Institute of Physics Conference Series, Vol. 990, FirstStars III, ed. B. W. O’Shea & A. Heger, D13O’Shea, B. W., & Norman, M. L. 2008, ApJ, 673, 14Osterbrock, D. E., & Ferland, G. J. 2006, Astrophysics of gaseousnebulae and active galactic nuclei, ed. Osterbrock, D. E. &Ferland, G. J.Palla, F., Salpeter, E. E., & Stahler, S. W. 1983, ApJ, 271, 632Palla, F., & Stahler, S. W. 1992, ApJ, 392, 667Papaloizou, J. C. B., & Terquem, C. 1995, MNRAS, 274, 987Peters, T., Klessen, R. S., Mac Low, M.-M., & Banerjee, R. 2010,ApJ, 725, 134Regan, J. A., & Haehnelt, M. G. 2009, MNRAS, 396, 343Regan, J. A., Johansson, P. H., & Haehnelt, M. G. 2014,MNRAS, 439, 1160Richling, S., & Yorke, H. W. 1997, A&A, 327, 317Sakurai, Y., Hosokawa, T., Yoshida, N., & Yorke, H. W. 2015,MNRAS, 452, 755Schleicher, D. R. G., Banerjee, R., Sur, S., Arshakian, T. G.,Klessen, R. S., Beck, R., & Spaans, M. 2010, A&A, 522, A115Schleicher, D. R. G., Palla, F., Ferrara, A., Galli, D., & Latif, M.2013, A&A, 558, A59Schneider, R., Omukai, K., Limongi, M., Ferrara, A., Salvaterra,R., Chieffi, A., & Bianchi, S. 2012, MNRAS, 423, L60Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337Siess, L., Forestini, M., & Bertout, C. 1997, A&A, 326, 1001Smith, R. J., Hosokawa, T., Omukai, K., Glover, S. C. O., &Klessen, R. S. 2012, MNRAS, 424, 457Springel, V. 2005, MNRAS, 364, 1105Stacy, A., & Bromm, V. 2013, MNRAS, 433, 1094Stacy, A., Greif, T. H., & Bromm, V. 2010, MNRAS, 403, 45—. 2012, MNRAS, 2508Stahler, S. W., Palla, F., & Salpeter, E. E. 1986, ApJ, 302, 590Stahler, S. W., Shu, F. H., & Taam, R. E. 1980, ApJ, 241, 637Susa, H. 2013, ApJ, 773, 185Susa, H., Hasegawa, K., & Tominaga, N. 2014, ApJ, 792, 32Tan, J. C., & Blackman, E. G. 2005, MNRAS, 362, 983Tan, J. C., & McKee, C. F. 2004, ApJ, 603, 383Tanaka, H., Takeuchi, T., & Ward, W. R. 2002, ApJ, 565, 1257Tanaka, K. E. I., Nakamoto, T., & Omukai, K. 2013, ApJ, 773,155Tielens, A. G. G. M., & Hollenbach, D. 1985, ApJ, 291, 722Truelove, J. K., Klein, R. I., McKee, C. F., Holliman, II, J. H.,Howell, L. H., & Greenough, J. A. 1997, ApJ, 489, L179Turk, M. J., Abel, T., & O’Shea, B. 2009, Science, 325, 601Turk, M. J., Oishi, J. S., Abel, T., & Bryan, G. L. 2012, ApJ,745, 154Vorobyov, E. I., & Basu, S. 2006, ApJ, 650, 956—. 2015, ApJ, 805, 115Vorobyov, E. I., DeSouza, A. L., & Basu, S. 2013, ApJ, 768, 131Whalen, D., & Norman, M. L. 2008, ApJ, 673, 664 assive Primordial Stars 27
Table A1. Ideal Cases for Homogeneous Clouds ConsideredCases N R · N θ · N φ ∆ t (10 yr) M ∗ ( M (cid:12) )HM2D 64 · · † · · · ·
64 3 70.8HM-ssk • · ·
64 2.6 71.1HM-HR2 † · ·
128 2.6 69.0Col. 2: numbers of grid cells, Col. 3: time duration of simulations. Col. 4: final stellar masses. The radial outer boundary is at r max =4 × AU, symmetry is assumed with respect to the equatorial plane ( θ = 0), and the initial cloud mass is 157 M (cid:12) for all cases. † : The suffixes have the same meanings as in Table. 1. • : Smaller sink with r min (cid:39)
10 AU ( r min = 30 AU for the other cases). s t e ll a r m a ss : M ∗ [ M fl ] (a) years ]10 -6 -5 -4 -3 -2 -1 a cc . r a t e s : ˙ M ∗ [ M fl y r − ] (b) Fig. A1.—
Evolution of the stellar mass (panel a) and accretion rates (panel b) for test cases starting with a homogeneous cloud(Table A1). The black solid and dashed lines represent the 2D axisymmetric cases with (HM2D) and without (HM2D-NF) UV feedback.The other lines represent 3D cases with different numerical settings: the fiducial case HM (blue), HM-ssk with a smaller sink (green), andHM-HR2 with doubled spatial resolution (magenta). The origin of time ( t = 0) marks the epoch of stellar birth after the runaway cloudcollapse.Yoon, S.-C., Dierks, A., & Langer, N. 2012, A&A, 542, A113Yorke, H. W., & Bodenheimer, P. 2008, in Astronomical Societyof the Pacific Conference Series, Vol. 387, Massive StarFormation: Observations Confront Theory, ed. H. Beuther,H. Linz, & T. Henning, 189Yoshida, N., Abel, T., Hernquist, L., & Sugiyama, N. 2003, ApJ,592, 645 Yoshida, N., Omukai, K., & Hernquist, L. 2007, ApJ, 667, L117—. 2008, Science, 321, 669Yoshida, N., Omukai, K., Hernquist, L., & Abel, T. 2006, ApJ,652, 6Zhu, Z., Hartmann, L., Nelson, R. P., & Gammie, C. F. 2012,ApJ, 746, 110Zinnecker, H., & Yorke, H. W. 2007, ARA&A, 45, 481APPENDIX A. IDEAL TEST CASES WITH A HOMOGENEOUS CLOUD
In addition to the cosmological cases studied above, we also perform simulations starting with idealized initialconditions: a zero-metallicity homogeneous gas cloud with the initial density ρ = 3 . × − cm − , temperature T = 200 K, and chemical compositions X H = 5 × − , X e = X HII = 10 − , where X i is the number ratio betweenthe species i and hydrogen nuclei. We assume that the cloud is initially rigidly rotating at an angular velocityΩ = 4 × − Hz, which is about 13 % of the Kepler value for the enclosed gas mass (cid:112) GM ( < r ) /r = (cid:112) πGρ / r . Since the initial Jeans length (cid:39) × AU is smaller than the cloud’s radial extension4 × AU, the cloud begins to gravitationally collapse in the well-known self-similar fashion (e.g., Omukai & Nishi1998). For the current test cases, we use our modified version of the mesh-based
PLUTO code from the beginning,unlike the cosmological cases for which the N-body/SPH code
GADGET -2 had been used for the early collapse stage (seeSection 2). The grid configuration is almost the same as described in Section 2, except that the current computationaldomain assumes mirror symmetry with respect to the equator throughout the evolution. We follow the collapsewithout a central sink cell until the central density reaches ∼ cm − and the minimum Jeans length becomes8 T. Hosokawa et al. R ∗ [ R fl ] (a) stellar radius stellar mass: M ∗ [ M fl ]10 S E UV [ s ec − ] (b) UV emissivity Fig. A2.—
Evolution of the stellar radius (panel a) and ionizing (EUV) photon emissivity (panel b) for the test cases starting with ahomogeneous cloud. The different lines represent the same cases as in Figure A1. In panel (a), the dotted line shows the mass-radiusrelationship for zero-age main-sequence stars. poorly resolved. We then introduce the sink cell to follow the subsequent evolution in the protostellar accretion stageusing either 2D or 3D versions of our SE-RHD hybrid code.For all examined 3D cases, the evolution is similar to that obtained for the cosmologcal cases D and E (see Sec-tion 4.2.1); M ∗ (cid:39) − M (cid:12) stars form after UV feedback caused by an expanding bipolar H ii region efficiently haltsmass accretion. For the 2D cases we used a so-called α -viscocity with α = 1 . M ∗ (cid:38) M (cid:12) as the UVemissivity rises in the late KH contraction stage (Fig. A2). The 3D cases show a comparable evolution, except for thelarge fluctuations of accretion rates seen in Figure A1-(b). This is due to the manner of mass and angular momentumtransport in self-gravitating circumstellar disks. We note that, for these idealized 3D cases, the assumed initial axialsymmetry is well preserved until the disk begins to grow during the accretion stage. Thus, even with a much less dis-turbed structure of the accretion envelope than for our cosmological cases, highly variable accretion occurs in general.Since, however, the mean accretion rate is relatively low for these test cases, the protostar continues to contract tothe ZAMS without being brought to the supergiant stage by accretion bursts. This is still the same even with thedoubled spatial resolution (case HM-HR2), though the scatter of the accretion rates is larger than the fiducial caseHM, as analogously seen for cosmological cases B and C (Section 4.3).For the above cases, the stellar radius is always much smaller than the default sink size 30 AU. We thus also examinethe effects of reducing the sink size to 10 AU (case HM-ssk). Figure A1 shows that the resulting accretion history ismore variable with the smaller sink. This trend is also consistent with our findings for cosmological case B, for which alarger sink resulted in less accretion variability (Fig. 21). Although the stellar mass growth is even more rapid than forthe other 3D cases for the first 10 years (Fig. A1-a), UV feedback limits the mass accretion once a bipolar H ii regionbegins to expand. Since the 10 AU sink is still larger than the stellar radius, we still need further studies dedicatedto resolving the flow connecting the disk inner edge and stellar surface. It is also possible that tight binaries withseparations much smaller than ∼∼