Gravitational Fragmentation of Extremely Metal-poor Circumstellar Discs
MMNRAS , 1–14 (2021) Preprint 16 February 2021 Compiled using MNRAS L A TEX style file v3.0
Gravitational Fragmentation of Extremely Metal-poor Circumstellar Discs
Kazuhiro Shima, ★ and Takashi Hosokawa † Department of Physics, Graduate School of Science, Kyoto University, Sakyo, Kyoto 606-8502, Japan
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
We study the gravitational fragmentation of circumstellar discs accreting extremely metal-poor ( 𝑍 ≤ − Z (cid:12) ) gas, performinga suite of three-dimensional hydrodynamic simulations using the adaptive mesh refinement code Enzo . We systematically followthe long-term evolution for 2 × years after the first protostar’s birth, for the cases of 𝑍 =
0, 10 − , 10 − , and 10 − Z (cid:12) . Weshow that evolution of number of self-gravitating clumps qualitatively changes with 𝑍 . Vigorous fragmentation induced by dustcooling occurs in the metal-poor cases, temporarily providing ∼
10 self-gravitating clumps at 𝑍 = − and 10 − Z (cid:12) . However,we also show that the fragmentation is a very sporadic process; after an early episode of the fragmentation, the number of clumpscontinuously decreases as they merge away in these cases. The vigorous fragmentation tends to occur later with the higher 𝑍 ,reflecting that the dust-induced fragmentation is most efficient at the lower density. At 𝑍 = − Z (cid:12) , as a result, the clumpnumber stays smallest until the disc fragmentation starts in a late stage. We also show that the clump mass distribution alsodepends on the metallicity. A single or binary clump substantially more massive than the others appear only at 𝑍 = − Z (cid:12) ,whereas they are more evenly distributed in mass at the lower metallicities. We suggest that the disc fragmentation should providethe stellar multiple systems, but their properties drastically change with a tiny amount of metals. Key words: accretion, accretion discs – hydrodynamics – methods: numerical – binaries: general – stars: formation – stars:protostars – early Universe.
The gravitational fragmentation of a circumstellar disc is a possibleprocess that provides multiple stellar systems in the present-day uni-verse (e.g. Kratter & Lodato 2016). In fact, observations are revealingdirect images of this process operating in some nearby star-formingregions (Tobin et al. 2016; Ilee et al. 2018). Theoretical studies in-vestigate necessary conditions for the disc fragmentation to yieldself-gravitating clumps (e.g. Gammie 2001; Takahashi et al. 2016).Numerical simulations demonstrate that the vigorous fragmentationoccurs particularly in an early evolutionary stage when the disc ac-cretes the gas infalling from a surrounding envelope (e.g. Vorobyov& Basu 2010; Machida et al. 2011; Tsukamoto et al. 2013; Oliva &Kuiper 2020).Numerical simulations are a powerful tool to even investigate thestar formation in the early universe, where only the pristine ( 𝑍 =
0) or extremely metal-poor (EMP, 𝑍 ≤ − Z (cid:12) ) gas exists (seeGreif 2015, for a review). Despite significant differences from thepresent-day star formation, particularly in the gas thermal evolution,a circumstellar disc embedded in the accretion envelope commonlyappears for this case (Tan & McKee 2004; Yoshida et al. 2008; Hiranoet al. 2014). Although there is some diversity, three-dimensional(3D) simulations broadly show that the disc becomes gravitationallyunstable and easily fragments for the primordial cases (Saigo et al.2004; Machida et al. 2008; Stacy et al. 2010; Clark et al. 2011; ★ E-mail: [email protected] (KS) † E-mail: [email protected] (TH)
Smith et al. 2011; Greif et al. 2012; Vorobyov et al. 2013; Susa 2013;Hosokawa et al. 2016; Stacy et al. 2016; Regan & Downes 2018;Sharda et al. 2019; Sugimura et al. 2020; Kimura et al. 2020).Despite the consensus that the disc fragmentation occurs in the pri-mordial star formation, it is still challenging to predict the statisticalproperties of multiple stellar systems that finally appear (e.g. Stacy &Bromm 2013). This is due to limitations in the numerics used in theliterature. In the primordial star formation, the effective adiabatic in-dex 𝛾 eff exceeds the critical value 4 / 𝑛 (cid:38) cm − ,with which the hydrostatic protostellar structure appears (Omukai& Nishi 1998; Yoshida et al. 2008). In most of previous simula-tions studying the disc fragmentation, however, such dense gas isnot spatially resolved to prevent the timestep from becoming verysmall. An often-used method is assuming an artificial stiff equationof state (EOS) above a threshold density (e.g. Machida & Nakamura2015; Hirano & Bromm 2017). Some authors employ the so-calledsink method, i.e., introduce a point mass which absorbs the nearbydense gas, representing an accreting protostar (e.g. Bate et al. 1995;Krumholz et al. 2004). These methods effectively mask the densegas with 𝑛 (cid:38) 𝑛 th , and the value of 𝑛 th differs in different studies.There is a trend that simulations assuming the higher 𝑛 th find thefragmentation in the earlier stage of the protostellar accretion (e.g.see appendix in Machida & Doi 2013). Since increasing 𝑛 th resultsin the short timestep, a higher-resolution simulation tends to onlyfollow the shorter-term evolution to save the computational cost. Theabove partly explains why the number of fragments reported in theliterature ranges from a few to ∼ © a r X i v : . [ a s t r o - ph . GA ] F e b K. Shima & T. Hosokawa such a disputed situation. Considering the almost scale-free natureof the governing equations of the fluid dynamics with self-gravity,he suggests that the apparently divergent results in the literature mayindicate similar time evolution described as 𝑁 c , b (cid:39) (cid:18) 𝑡 ff , ad 𝑡 ff , th 𝛥𝑡 (cid:19) . ≡ 𝛥 ˜ 𝑡 . , (1)where 𝑁 c , b is the number of the gravitationally-bound clumps ina given snapshot, 𝛥𝑡 is the elapsed time since the first appearanceof a protostar in the unit of year, and 𝑡 ff , th and 𝑡 ff , ad are the free-fall timescales defined with the threshold number density 𝑛 th and 𝑛 ad = cm − , and 𝛥 ˜ 𝑡 ≡ √︁ 𝑛 th / 𝑛 ad 𝛥𝑡 . Susa (2019) shows that thesimulation results by different authors roughly follow equation (1)once scaled with assumed values of 𝑛 th , though associated with anorder-of-magnitude scatter.Although it is still unknown why the disc fragmentation for theprimordial cases is well described by the simple relation such as equa-tion (1), an important fact is that a barotropic EOS with 𝛾 eff (cid:39) . 𝑛 (cid:46) cm − (Omukai & Nishi 1998). In fact, sev-eral simulations study the disc fragmentation assuming the samebarotropic EOS with 𝛾 eff = . 𝑛 ≤ 𝑛 th , resulting in the evo-lution described by equation (1) (Susa 2019). Since 𝑛 th is only thecharacteristic quantity for this case, the simple scaling of equation(1) may be convincing. This suggests that the disc fragmentationwith a different EOS should provide different evolution of 𝑁 c , b . Forinstance, it is well known that adding a tiny amount of heavy el-ements and dust grains alters the EOS of a collapsing cloud (e.g.Omukai 2000; Bromm et al. 2001; Omukai et al. 2005; Schneideret al. 2006; Smith et al. 2008; Jappsen et al. 2009; Schneider et al.2012; Safranek-Shrader et al. 2014; Chiaki et al. 2015, 2016). Tanaka& Omukai (2014) investigate the evolution of the circumstellar disc inmetal-poor environments developing one-dimensional semi-analyticmodels. They predict that the discs with 𝑍 ∼ − − − Z (cid:12) are sub-ject to the efficient dust cooling and are more unstable than those forthe primordial cases. Machida & Nakamura (2015) consider the discfragmentation with various metallicities 0 ≤ 𝑍 ≤ (cid:12) , performinga suite of 3D numerical simulations. They find qualitative differencesbetween the cases with 𝑍 (cid:46) − Z (cid:12) and 𝑍 (cid:38) − Z (cid:12) ; the vigor-ous disc fragmentation only occurs for the former metal-poor cases.Whereas Machida & Nakamura (2015) use the metallicity-dependentbarotropic EOS, Chiaki & Yoshida (2020) recently report 3D sim-ulations of the disc fragmentation solving the energy equation withrelevant thermal processes coupled with a non-equilibrium chemicalnetwork. They find that for the metal-poor cases with 𝑍 ≤ − Z (cid:12) ,the disc fragmentation does not necessarily prevent the mass growthof the most massive protostar as many clumps are short-lived owingto the frequent merger or tidal disruption events. Although the abovestudies suggest the metallicity-dependence of the disc fragmentation,they both only follow the short-term evolution for ∼
100 yr since thefirst emergence of a protostar.In this paper, we investigate the long-term evolution of thedisc fragmentation with various metallicities in the range of 𝑍 ≤ − Z (cid:12) : 𝑍 =
0, 10 − , 10 − , and 10 − Z (cid:12) . We systematically per-form 3D hydrodynamic simulations that follow the evolution in theearly collapse stage and subsequent accretion stage at the differ-ent metallicities. We examine the metallicity-dependence of the discfragmentation that occurs during the first 2 × yr of the protostellaraccretion. We track the evolution of the number of self-gravitatingclumps for each case, as compiled for the primordial cases by Susa(2019). We consider how the simple scaling given by equation (1)may be applicable or modified for the low-metallicity cases. We also n [cm ]10 T [ K ] Z = 0 Z = 10 Z Z = 10 Z Z = 10 Z Figure 1.
Barotropic equation of state (EOS) assumed for our simulations atdifferent metallicities. The different colors represent different metallicities of 𝑍 = − (orange), 10 − (green), and 10 − Z (cid:12) (brown). Variation ofthe temperature as functions of the number density is evaluated using the one-zone modeling of a collapsing cloud (Omukai et al. 2005; Susa et al. 2015).The steep parts for 𝑛 (cid:38) cm − correspond to the stiff EOS assumed as thepressure floor (also see text). The dashed lines represent the original one-zonemodeling results for which no pressure floor is imposed. investigate how other clump properties, such as their mass distribu-tion, change with increasing the metallicity.The rest of the paper is organized as follows. We describe thenumerical simulation methods in Section 2. We show our simulationresults in Section 3. We finally provide discussion and concludingremarks in Sections 4 and 5. In this study, we model the gas thermal evolution with the pre-calculated barotropic EOS as in Machida & Nakamura (2015), usingthe results of the so-called one-zone modeling of a collapsing cloud(e.g. Omukai 2000). Fig. 1 shows the EOS we use, as variationsof the temperature against the density at different metallicities of 𝑍 =
0, 10 − , 10 − , and 10 − Z (cid:12) . The curve for 𝑍 = < 𝑍 ≤ − Z (cid:12) because the resulting EOS is almost identicalto that for the primordial case. The present-day temperature of thecosmic microwave background is assumed for these cases.In Fig. 1, the local minima of the curves correspond to the pointsof 𝛾 eff =
1, where the fragmentation is expected to be enhanced (e.g.Larson 1985, 2005; Jappsen et al. 2005). We see that the curves exceptfor 𝑍 = 𝑛 ∼ cm − .The former is caused by molecular cooling and the latter by coolingvia dust thermal emission. Since 𝑛 (cid:38) cm − in circumstellar discswe consider, only the dust-induced fragmentation is relevant to ourcases.Although the barotropic EOS is obtained for the evolution duringthe early collapse stage, it has been applied to study the disc fragmen-tation in the literature. Clark et al. (2011) show that, at least for theprimordial cases, applying the barotropic EOS tends to result in thelower disc temperature and thus more fragmentation than solving thethermal and chemical processes in a time-dependent hydrodynamiccode (see also Matsukoba et al. 2021). Nonetheless, Susa (2019) MNRAS , 1–14 (2021) etal-poor disc fragmentation Table 1.
Cloud initial parameters: radius, mass, temperature, and averagedensity. Profile 𝑟 𝑀 𝑇 ¯ 𝑛 (pc) (M (cid:12) ) (K) (cm − )cloud A 1.3 1.1 ×
195 3.8 × cloud B 2.5 × − ×
35 1.9 × Table 2.
Cases considered.Z initial profile 𝛽 (Z (cid:12) )0 cloud A 0.03, 0.06, 0.0910 − cloud A 0.03, 0.06, 0.0910 − cloud A 0.03, 0.06, 0.0910 − cloud B 0.03, 0.06, 0.09 has found that the previous simulations with both approaches showsimilar evolution of the clump number as described by equation (1).Our current method is thus valid for the aim of the current work,whereas quantitative effects of relying on the barotropic EOS is tobe examined in future studies. Our simulations use the adaptive mesh refinement (AMR) hydro-dynamics code,
Enzo (Bryan et al. 2014). The gas is evolved withself-gravity solving the Poisson equation using a multi-grid solver.Regarding the solver of the hydrodynamics, we adopt a three-dimensional implementation of the
Zeus hydro-code (Stone & Nor-man 1992a,b).As for the initial states, we assume two different cloud properties:cloud A for cases with 𝑍 ≤ − Z (cid:12) and cloud B with 10 − Z (cid:12) (seeTables 1 and 2), both of which rigidly rotate at the angular frequency Ω rot . The cloud takes the density profile of a Bonnor-Ebert (BE)sphere (Ebert 1955; Bonnor 1956); a hydrostatic isothermal self-gravitating sphere of gas that is confined by its external pressure.While such a profile is derived analytically, cosmological simulationsshow that it is realized in a primordial star-forming cloud in the so-called "loitering" phase (Bromm et al. 2002; Hirano et al. 2014). TheBE profile is parametrized by the central density 𝑛 , c and temperature 𝑇 , for which we use the values at the loitering point or the localminimum at 𝑛 ∼ cm − in Fig. 1. We use 𝑛 , c = . × cm − and 𝑇 =
195 K for the primordial case (cloud A). We assume thatthe critical BE profile continues until the cloud radius 𝑟 , where theenclosed mass is (cid:39) M (cid:12) , a typical value for the primordial clouds.We also assume the cloud is embedded in the homogeneous mediumthat provides the constant external pressure for 𝑟 > 𝑟 . We enhancethe density by 10 % within the cloud to cause the collapse, and wefurther add a perturbation of 𝑚 = 𝜌 ( 𝑟 )( + 𝛿 cos 2 𝜙 ) ,where 𝜌 is the unperturbed density, 𝛿 the perturbation amplitude,and 𝜙 the azimuthal angle around the rotation (or 𝑧 -) axis. We onlyconsider the cases with 𝛿 = . 𝑍 ≤ − Z (cid:12) are almost identicalto the primordial case for 𝑛 (cid:46) cm − . We thus use the initialcondition of the cloud A also for the low-metallicity cases with 𝑍 ≤ − Z (cid:12) . We only consider a different initial cloud configurationfor 𝑍 = − Z (cid:12) , with which the evolution departs from the othercases at 𝑛 ∼
100 cm − . To characterize the BE sphere for this case, we use 𝑛 , c = . × cm − and 𝑇 =
35 K, the values at the localminimum (cloud B). The following procedure to construct the initialstate is the same as for cloud A.We assume idealized clouds with artificial density perturbationsas the initial conditions, and we also ignore the magnetic fields andturbulence for simplicity. We remark that our simulations followthe evolution of the protostellar accretion for 2 × yr (see Sec-tion 2.3), and the duration corresponds to the free-fall timescale at 𝑛 ∼ cm − . This is much higher than the initial central value,meaning that the disc only accretes the gas coming from a denseaccretion envelope that is set during the cloud collapse. We also notethat the disc fragmentation is not only a process that provides numer-ous self-gravitating clumps. Filamentary structure of star-formingclouds is known to cause fragmentation in general, particularly for thelow-metallicity cases where the dust cooling operates (e.g. Tsuribe& Omukai 2006; Clark et al. 2008; Dopcke et al. 2011, 2013; Chiakiet al. 2016; Sugimura et al. 2017). Our simple setup only allows thefilamentary structure to develop in the circumstellar discs. We isolatethe effects of the disc fragmentation with different metallicities in thecurrent work.The simulation box size is 3 pc for cloud A and 0.75 pc for cloudB on a side, which are (cid:39) 𝑟 . Thebox is covered by 128 root grids initially. Higher grid levels are addedduring the evolution by the adaptive mesh refinement technique, re-ducing the cell size by a factor of two. Spatial cells are refined basedon the requirement that the Jeans length must not fall below 32 cells(e.g. Federrath et al. 2011). We also confirm the numerical conver-gence of our results by performing additional simulations resolvingthe Jeans length by 16 and 64 cells (see Section 4.2). At the maximumrefinement level 16 for cloud A and 14 for cloud B (the minimum cellsize is the same in both clouds), where the Jeans criteria inevitablymust break, we introduce a pressure floor to halt the collapse at a fi-nite density, preventing individual cells from becoming unphysicallymassive. We realize the pressure floor by assuming the stiff EOS with 𝛾 ad = 𝑛 (cid:38) cm − (Takahira et al. 2014), as illustrated inFig. 1. The floor density 10 cm − is similar to that used in Susa(2019), who has only considered the primordial case and shown theevolution well described by equation (1). We calculate the evolution of 12 models in total with 4 differentmetallicities and 3 different degrees of the initial rotation. Thesemodels are summarized in Table 2. We represent the cloud’s rota-tional degree with 𝛽 , or the ratio of the rotation energy 𝐸 rot to thegravitational energy 𝐸 grav 𝛽 ≡ 𝐸 rot (cid:12)(cid:12) 𝐸 grav (cid:12)(cid:12) = Ω 𝑟 𝐺 𝑀 . (2)We take 𝛽 = .
06 as a standard value, which corresponds to theinitial angular frequency Ω rot (cid:39) . × − s − for cloud A and Ω rot (cid:39) . × − s − for cloud B. We also consider the cases with 𝛽 = .
03 and 0.09 for comparisons. Whereas our examined range of 𝛽 is almost the same as in Susa (2019), it is much higher than that inMachida & Nakamura (2015) by more than an order of magnitude.Cosmological simulations show that primordial clouds typically has 𝛽 ∼ . MNRAS , 1–14 (2021)
K. Shima & T. Hosokawa reaches the maximum level. After that, we follow the evolution ofthe protostellar accretion for 2 × yr, during which the disc frag-mentation occurs. The duration is comparable to that in Susa (2019),but it is about ten times longer than in Machida & Nakamura (2015)and Chiaki & Yoshida (2020). Whereas Susa (2019) only considersthe primordial cases, we study the metallicity-dependence of the discfragmentation during the similar long-term evolution in 3D. As described in Section 2.2, we artificially halt the cloud collapseby using the stiff EOS for 𝑛 (cid:38) cm − . Otherwise, the collapsefurther continues, and the timestep becomes smaller and smaller.Following the long-term evolution of the disc fragmentation becomescomputationally infeasible for such a case. An alternative method toachieve the same purpose is using the Lagrangian sub-grid modelsuch as the sink cell/particle method. However, we do not resort tothis in the current work.A disadvantage of the sink method is that the results may dependon the details of the implementation. For instance, we have to relyon the unknown criteria for the sink creation, accretion onto thesink, and mergers between them, though some possible solutionshave been developed (e.g. Federrath et al. 2010; Hubber et al. 2013).Since we also implemented the sink method proposed by Federrathet al. (2010) to Enzo (Shima et al. 2018), we actually performedpreliminary simulations of the disc fragmentation using our modifiedversion of the code. However, it turned out that the simulation withthe sink method is more computationally expensive than that withthe stiff EOS. We thus adopt the current method of the stiff EOS,which is simpler and more efficient than the sink method for ourspecific cases. Susa (2019) has investigated the effects of using thestiff EOS and sink methods in comparisons in smoothed particlehydrodynamics (SPH) simulations. Fortunately, equation (1) welldescribes the evolution of the disc fragmentation observed in bothcases.Since our method defines neither self-gravitating clumps nor pro-tostars, we need a method to detect them. We save the simulationdata every 10 years to follow the evolution of self-gravitating clumpsafter the pressure reaches the floor value. To identify the clumps ateach snapshot, we make use of the finder implemented in yt (Smithet al. 2009; Turk et al. 2011), which employs a contouring algorithmto identify topologically disconnected structures. We begin with thedensity threshold of 𝑛 iso = cm − to find iso-density contoursand identify clump candidates. We continually multiply 𝑛 iso by afactor of 10 and apply the clump finding algorithm again. If a clumpidentified with the lower-density contour turns out to contain twoclumps with the higher-density contour, we regard the number ofclumps as two. We repeat the whole procedure until 𝑛 iso exceedsthe maximum density in the given snapshot. We estimate the massof a clump by summing up the gas contained in a disconnected iso-density contour at the lowest level, above which there are no furthersub-clumps inside.For each clump, we check whether it is gravitationally bound ornot as follows. We consider 𝑊 + 𝐾 + 𝑈 < 𝑊 ( < ) is the total gravitationalenergy, 𝐾 the total kinetic energy, and 𝑈 the total thermal energywithin the clump. We evaluate 𝑊 by summing up the gravitationalbinding energy between all two point cells inside the clump, 𝑊 = − ∑︁ 𝑖, 𝑗,𝑖 ≠ 𝑗 𝐺𝑚 𝑖 𝑚 𝑗 | r 𝑖 − r 𝑗 | , (3)where 𝑚 𝑖 is each cell mass and the 𝑟 𝑖 is the position of the cell. The total kinetic energy is 𝐾 = ∑︁ 𝑖 𝑚 𝑖 (( 𝑢 𝑖 − 𝑢 𝑐 ) + ( 𝑣 𝑖 − 𝑣 𝑐 ) + ( 𝑤 𝑖 − 𝑤 𝑐 ) ) , (4)where ( 𝑢 𝑖 , 𝑣 𝑖 , 𝑤 𝑖 ) is the velocity in each cell inside the clump,and ( 𝑢 𝑐 , 𝑣 𝑐 , 𝑤 𝑐 ) is the velocity of the clump’s center-of-mass. Wenormally only count the gravitationally bound objects, but we alsoinvestigate how much the clump number increases if we skip thebinding check in Section 4.2. Hereafter, we define the origin of theelapsed time 𝛥𝑡 = 𝑛 iso = cm − roughly corresponds to the values where the EOS curves for 𝑍 = − and 10 − Z (cid:12) take the local minima owing to the dust cooling (seeFig. 1). We have tested different initial choices of 𝑛 iso = cm − and 10 cm − . The iso-density contours with 𝑛 iso = cm − cover a large part of the disc rather than the individual clumps. Start-ing the analysis with 𝑛 iso = cm − , on the other hand, returnsalmost the same result as with 𝑛 iso = cm − . We have also testeda different 𝑛 iso -multiplying factor 2 instead of 10, with which theresults hardly change. In what follows we present our simulation results. We first describethe metallicity-dependence of the evolution, considering the caseswith the rotation parameter fixed at 𝛽 = .
06 in Sections 3.1 and 3.2.We next study the effects of varying the parameter 𝛽 in Section 3.3,where the other cases with 𝛽 = .
03 and 0 .
09 are presented. Wefinally investigate how the mass distribution of the self-gravitatingclumps varies with the different metallicities in Section 3.4. ( 𝑍 = ) Fig. 2 shows the face-on images of the disc fragmentation occurringfor the primordial case with 𝛽 = .
06. The snapshots are taken atepochs 𝛥𝑡 = yr after the first clump appears. Thecenter of each panel corresponds to the most massive cell. The crossesin the density maps represent the mass centers of the clumps identifiedby the finder starting with 𝑛 iso = cm − (see Section 2.4). Thenumber of gravitationally bound clumps 𝑁 c , b and maximum clumpmass 𝑚 c , max are also presented in the upper left corner of eachtop panel. The earliest 𝛥𝑡 =
100 yr snapshot clearly shows the discstructure with the two spiral arms, which reflect the initial densityperturbation of 𝑚 = 𝛥𝑡 =
300 yr snapshot, there is still the samebinary found in 𝛥𝑡 =
100 yr near the center, but two more clumps alsoappear in the outer part of the disc after the additional fragmentation.The latter two clumps migrate inward after that, and they merge withthe former clumps. Meanwhile, another clump appears because ofthe fragmentation that occurs a few ×
100 AU away from the center.The last 𝛥𝑡 = yr snapshot consequently shows one clump in theouter large orbit and the other two in a central tight binary system.The binary separation is (cid:39)
24 AU in this epoch. After all, there arealmost always a few clumps during the evolution followed by oursimulation. We also consider the evolution of the clump number inSection 3.2, with reference to Susa (2019).
MNRAS , 1–14 (2021) etal-poor disc fragmentation y [ A U ] t = 100yr N c,b = 2 M c,max = 2.13M t = 300yr N c,b = 4 M c,max = 5.16M t = 1000yr N c,b = 3 M c,max = 19.04M400 200 0 200 400 x [AU]4002000200400 y [ A U ]
400 200 0 200 400 x [AU] 400 200 0 200 400 x [AU] 10 nu m b e r d e n s i t y [ c m ] t e m p e r a t u r e [ K ] Figure 2.
Images of the density-weighted projection maps (through the z-axis) of the number density (top) and temperature (bottom) for the case of 𝑍 = 𝛽 = .
06. The left, middle, and right panels show the snapshots at 𝛥 𝑡 =
100 yr, 300 yr, and 10 yr after the first appearance of a gravitationally bound clump. Weuse the color scales covering the maximum and minimum values that appear through the snapshots for both the top and bottom panels. The center of each panelcorresponds to the position of the densest cell. The white crosses in the upper panels indicate the mass-center positions of the gravitationally bound clumpsidentified by the finder with a threshold of 10 cm − (see Section 2.4). The red cross denotes the most massive clump at each snapshot. The disc radially spreadsfrom left to right as it accretes the gas infalling from the envelope with the higher angular momentum. t [kyr]10203040506070 M t o t a l [ M ] Z = 0 Z = 10 Z Z = 10 Z Z = 10 Z Figure 3.
Time evolution of the total mass of self-gravitating clumps detectedby the finder with the initial iso-contour density 𝑛 iso = cm − . The linecolors represent the same cases with different metallicities as in Fig. 7. As shown in Fig. 3, the total mass of the clumps exceeds 50 M (cid:12) after the first 2 × yr, indicating the mean accretion rate of (cid:39) . (cid:12) yr − , a typical value in the primordial star formation(e.g. Hirano et al. 2014). The most massive clump accretes ∼
20 M (cid:12) of the gas by the epoch of 𝛥𝑡 = yr. Since the clumps repre-sent accreting protostars, radiation emitted from such massive ones may affect the evolution of the disc. The primary feedback is causedby stellar ultra-violet radiation, creating the photoionised and pho-todissociation regions around a protostar (e.g. McKee & Tan 2008;Hosokawa et al. 2011, 2016; Fukushima et al. 2020). Our simulations,where these effects are assumed to be negligible, only consider thedisc fragmentation before the radiative feedback begins to operate.We note that our assumption is more justified for the low-metallicitycases shown below, where the maximum clump mass is much lowerthan the primordial case at the end of the simulations. 𝑍 = − Z (cid:12) Fig. 4 presents how the disc fragmentation proceeds for the case of 𝑍 = − Z (cid:12) and 𝛽 = .
06. We see that the growing disc easily frag-ments, forming many self-gravitating clumps. Such basic evolutionmay appear to be similar to the primordial case, but they are quantita-tively very different. The snapshot of 𝛥𝑡 =
100 yr shows that there aremore than ten fragments along the outstanding spiral structure. Wehave confirmed, by means of the method described in Section 2.4,that these objects are all gravitationally bound. The clump number of 𝑁 c , b =
13 is much more than that for the primordial case at the sameepoch, 𝑁 c , b = 𝑍 = 𝑍 = − Z (cid:12) shows the remarkable temperature MNRAS , 1–14 (2021)
K. Shima & T. Hosokawa y [ A U ] t = 100yr N c,b = 13 M c,max = 1.26M t = 300yr N c,b = 8 M c,max = 3.22M t = 1000yr N c,b = 3 M c,max = 14.31M400 200 0 200 400 x [AU]4002000200400 y [ A U ]
400 200 0 200 400 x [AU] 400 200 0 200 400 x [AU] 10 nu m b e r d e n s i t y [ c m ] t e m p e r a t u r e [ K ] Figure 4.
Same as Fig. 2 but for the case of 𝑍 = − Z (cid:12) . There are more than 10 clumps at the first snapshot 𝛥 𝑡 =
100 yr, but the clump number monotonicallydecreases until the last snapshot at 𝛥 𝑡 = yr. decline for 𝑛 (cid:38) cm − and the local minimum at 𝑛 ∼ cm − .The density and temperature along the spiral arm take these mini-mum values, suggesting the dust-induced fragmentation occurs. Weseparately focus on this feature later in Section 4.1. The subsequent 𝛥𝑡 =
300 yr and 10 yr snapshots show that the number of clumpscontinuously decreases with time, which differs from the primor-dial case. The clumps undergo the complex orbital evolution throughgravitational interaction with each other, and most of them mergeaway during that. As a result, only three clumps survive at the lastsnapshot of 𝛥𝑡 = yr. This number is coincidentally the same asthat for the primordial case at the same epoch.Figs. 3 and 4 indicate that the clumps accrete the gas of (cid:39)
20 M (cid:12) in total, of which the primary one dominates (cid:39)
14 M (cid:12) at the epochof 𝛥𝑡 = yr. The total mass is lower than the primordial case by afactor of a few. The disc-star system accretes the gas from the centralpart of the envelope during 𝛥𝑡 = 𝑡 ff at 𝑛 ∼ cm − . Fig. 1 presents that the temperature at 𝑛 ∼ cm − differs by a factor of two among 𝑍 = 𝑍 = − Z (cid:12) cases. Giventhat the total accretion rate depends on the envelope temperature as (cid:164) 𝑀 ∝ 𝑇 . , the difference in the mass growth histories agrees withour assumed EOS. 𝑍 = − Z (cid:12) Fig. 5 represents the case of 𝑍 = − Z (cid:12) and 𝛽 = .
06. Theevolution for this case quantitatively differs from both of 𝑍 = 𝑍 = − Z (cid:12) cases described above. The first 𝛥𝑡 =
100 yr snapshotshows that there are three fragments along in a straight line due to theinitial perturbation. Whereas such density structure is more or less similar to the primordial case (Fig. 2), the temperature distributionlooks very different owing to the different EOS (Fig. 1). We see a coldpart at 𝑟 (cid:46)
30 AU, which corresponds to the temperature decline for 𝑛 (cid:38) cm − . The outer part shows almost the same temperatureat (cid:39)
400 K, reflecting the plateau for 10 cm − (cid:46) 𝑛 (cid:46) cm − .Subsequently, the evolution similar to 𝑍 = − Z (cid:12) case continues.The snapshot of 𝛥𝑡 =
300 yr displays that there are thirteen clumpsnear the spiral arms. The density and temperature along the spiralarms correspond to the values at the local minimum of the EOScurve at 𝑛 ∼ . cm − . This fact suggests that the dust-inducedfragmentation yields the clumps as in 𝑍 = − Z (cid:12) case, but it occurslater. The clump number slightly decreases by the last snapshot of 𝛥𝑡 = yr because of a few merger events. They are more sparselydistributed than in the previous snapshot after complex gravitationalinteractions. The clump number further continues to decrease untilthe end of the simulation, 𝛥𝑡 = × yr (also see Section 3.2 below).We remark that in the bottom panels the clumps are more out-standing than 𝑍 = − Z (cid:12) cases describe above, representedby the bright (or hot) spots surrounded by the dark (or cold) regions.We interpret this trend as follows. As seen in Fig. 1, the EOS curvesfor 𝑍 ≥ − Z (cid:12) all converge to the same line for 𝑛 (cid:38) cm − ,where 𝛾 eff (cid:39) / 𝑍 . Since the Jeans length is in propor-tion to 𝑛 − / with 𝛾 eff = /
5, the lower-density core has the largersize. Moreover, the disc size at a given epoch 𝛥𝑡 is systematicallysmaller at the higher 𝑍 . Recall that the disc only accretes the gasfrom a central part of the envelope where 𝑡 ff < 𝛥𝑡 . The correspond-ing part is more compact at the higher metallicity because the Jeanslength is smaller with the lower temperature. Therefore, the higher- 𝑍 MNRAS000
5, the lower-density core has the largersize. Moreover, the disc size at a given epoch 𝛥𝑡 is systematicallysmaller at the higher 𝑍 . Recall that the disc only accretes the gasfrom a central part of the envelope where 𝑡 ff < 𝛥𝑡 . The correspond-ing part is more compact at the higher metallicity because the Jeanslength is smaller with the lower temperature. Therefore, the higher- 𝑍 MNRAS000 , 1–14 (2021) etal-poor disc fragmentation y [ A U ] t = 100yr N c,b = 3 M c,max = 0.66M t = 300yr N c,b = 13 M c,max = 1.16M t = 1000yr N c,b = 10 M c,max = 3.86M400 200 0 200 400 x [AU]4003002001000100200300400 y [ A U ]
400 200 0 200 400 x [AU] 400 200 0 200 400 x [AU] 10 nu m b e r d e n s i t y [ c m ] t e m p e r a t u r e [ K ] Figure 5.
Same as Fig. 2 but for the case of 𝑍 = − Z (cid:12) . The clump number takes the maximum at the second snapshot 𝛥 𝑡 =
300 yr (middle panels), which islater than the case of 𝑍 = − Z (cid:12) (cf. Fig. 4). y [ A U ] t = 100yr N c,b = 1 M c,max = 0.21M t = 300yr N c,b = 1 M c,max = 0.22M t = 1000yr N c,b = 3 M c,max = 1.68M 0.44 0.01100 50 0 50 100 x [AU]10050050100 y [ A U ]
100 50 0 50 100 x [AU] 100 50 0 50 100 x [AU] 10 nu m b e r d e n s i t y [ c m ] t e m p e r a t u r e [ K ] Figure 6.
Same as Fig. 2 but for the case of 𝑍 = − Z (cid:12) . In the right top panel, the masses of the companion clumps are also described in the unit of M (cid:12) .MNRAS , 1–14 (2021) K. Shima & T. Hosokawa t [yr]10 N c , b = 0.06 Z = 0 Z = 10 Z Z = 10 Z Z = 10 Z10 t [yr] Figure 7.
Time evolution of the number of gravitationally bound clumpswith different metallicities. The line colors represent the cases with differentmetallicities, whose snapshots are presented in Figs. 2 - 6. The same cloud’srotation parameter 𝛽 = .
06 is assumed for these cases. The bottom and tophorizontal axes represent the time elapsed since the first detection of a boundclump 𝛥 𝑡 and 𝛥 ˜ 𝑡 ≡ √︁ 𝑛 th / 𝑛 ad 𝛥 𝑡 (also see equation 1), where 𝑛 th = cm − and 𝑛 ad = cm − . The black dashed line represents equation (5). Theblue-shaded background denotes the area within × × / disc should accrete the gas with the lower angular momentum, whichexplains its smaller size. The above trend becomes more prominentin the case of 𝑍 = − Z (cid:12) described below. 𝑍 = − Z (cid:12) Fig. 6 shows the evolution for the case of 𝑍 = − Z (cid:12) and 𝛽 = . ∼
300 AU on aside, which is much smaller than in Fig. 5. In this case, the snapshotsfor 𝛥𝑡 =
100 yr and 300 yr both present only one self-gravitatingclump. Although we see the spiral arms develop in the disc, it doesnot cause the vigorous fragmentation for 𝛥𝑡 ≤
300 yr, in contrastto the lower-metallicity cases described above. Such a difference iswell understood with Fig. 1, where the EOS curve for 𝑍 = − Z (cid:12) monotonically increases for 𝑛 (cid:38) cm − . Since the density withinthe disc takes 𝑛 ∼ − cm − , the disc temperature is nearlyconstant at a few ×
10 K. The temperature substantially increasesonly in the interior of an adiabatic core, i.e., for 𝑛 (cid:38) cm − .Accordingly, the central clump is relatively very hot and large againstthe surrounding disc in Fig. 6. The disc fragmentation finally occursafter the epoch of 𝛥𝑡 = , yr. As a result, there are three clumpsin the 𝛥𝑡 = yr snapshot. The bottom panel for this epoch showsthat the spiral arm becomes relatively colder than the surroundinggas. Note that the density just outside the disc gradually drops as itaccretes the gas from the envelope. The EOS curve for 𝑍 = − Z (cid:12) has a shallow minimum at 𝑛 ∼ cm − (Fig. 1), whose featureappears near the outer edge of the disc in the 𝛥𝑡 = yr snapshot.Fig. 3 shows that the mean total accretion rate onto the clumpsfor 𝛥𝑡 = × yr is as low as ∼ − M (cid:12) yr − for 𝑍 = − Z (cid:12) .This is lower than that for the primordial case by a factor of (cid:39) . 𝑍 = − Z (cid:12) case only differs from the other cases. However,its effects on the evolution of the protostellar accretion we consider should be limited because we only focus on the central dense partwith 𝑛 (cid:38) cm − . Such dense gas appears in a late stage of the run-away cloud collapse, which converges to the same similarity solutionregardless of different low-density initial states (e.g. Larson 1969;Yahil 1983; Omukai & Nishi 1998). Fig. 7 summarizes the metallicity-dependent evolution of the numberof clumps for the cases described in Section 3.1. In Fig. 7, we alsooverlay the scaling relation obtained by Susa (2019) for the primordialcases. Since 𝑛 th = cm − for our cases, we rewrite equation (1)as 𝑁 c , b (cid:39) . (cid:18) 𝛥𝑡
100 yr (cid:19) . . (5)We first compare our primordial case to this relation. As seen in thefigure, the corresponding blue line goes slightly below the dashedline of equation (5). Our case indicates that the clump number staysalmost constant at 𝑁 c , b (cid:39) 𝛥𝑡 >
100 yr, whereas equation (5)predicts the monotonic increase of 𝑁 c , b . Our case remains withinthe blue-shaded area, typical scatter of simulation results previouslyreported by different authors, until 𝛥𝑡 ∼ yr. This fact suggeststhat our simulation predicts relatively smaller 𝑁 c , b than other studies.We note that our setup is almost the same as in Susa (2019), exceptthat the basic simulation methods are different; AMR in our studyand SPH in Susa (2019). In fact, no previous 3D AMR simulationresults have been tested against equation (1), and many of previousstudies compiled in Susa (2019) are SPH simulations. Since our aimis investigating the metallicity-dependence of the disc fragmentation,we do not further study what causes the difference in detail. In whatfollows, we suggest that varying the initial cloud rotation (or 𝛽 ) orspatial resolution does not resolve the discrepancy (see Sections 3.3and 4.2).Next, we consider the evolution in the EMP cases described inSections 3.1.2-3.1.4. Fig. 7 clearly shows that the evolution of 𝑁 c , b varies with different metallicities. None of them show the evolutionsimilar to the primordial case. Among them, however, the cases of 𝑍 = − Z (cid:12) and 10 − Z (cid:12) show the common feature; the clumpnumber takes the local maximum 𝑁 c , b , max in an early stage, and itcontinues to decrease afterward. The epoch of 𝑁 c , b , max is somewhatdelayed with increasing the metallicity, as already mentioned in Sec-tion 3.1.3. The evolution of 𝑁 c , b for these cases qualitatively differsfrom the primordial case. The enhancement of 𝑁 c , b is caused by thedust-induced disc fragmentation, which agrees with the semi-analyticmodels developed by Tanaka & Omukai (2014). However, our sim-ulations suggest that fragmentation is a very sporadic process. Afterthe vigorous fragmentation, when 𝑁 c , b takes the maximum, the fur-ther fragmentation ceases for a long time, during which many clumpsmove around and undergo mutual gravitational interactions. The or-dered spiral arm do not grow until 𝑁 c , b substantially drops. As aresult, 𝑁 c , b for these cases becomes comparable to that for 𝑍 = years.In contrast, the case of 𝑍 = − Z (cid:12) displays the opposite trend; 𝑁 c , b remains remarkably small until the efficient fragmentation startsat 𝛥𝑡 (cid:39)
500 yr. The similar trend has been reported by Machida &Nakamura (2015), who study the initial 200 years of the protostellaraccretion with various metallicities. They show that the vigorous discfragmentation only occurs for the cases with 𝑍 (cid:46) − Z (cid:12) , whichagrees with our results if only paying attention to the early evolution.Our long-term simulation shows that 𝑁 c , b starts to increase later, MNRAS , 1–14 (2021) etal-poor disc fragmentation N c , b Z = 0= 0.03= 0.06= 0.09 10 Z = 10 Z= 0.03= 0.06= 0.0910 t [yr]10 N c , b Z = 10 Z= 0.03= 0.06= 0.09 10 t [yr] Z = 10 Z= 0.03= 0.06= 0.0910 t [yr] 10 t [yr] Figure 8.
Same as Fig. 7 but for varying the cloud rotation parameter 𝛽 . Each panel shows the evolution at a given metallicity as indicated at the upper rightcorner. The thin, normal, and thick lines represent the cases with 𝛽 = .
03, 0.06, and 0.09, respectively. because of the delayed dust-induced fragmentation at 𝑍 = − Z (cid:12) .Although we only follow the evolution for 𝛥𝑡 < × yr, the clumpnumber may further increase in the later stage. While we have fixed the cloud rotation parameter at 𝛽 = .
06 forthe models described in Sections 3.1 and 3.2, we here consider themetallicity-dependent evolution with 𝛽 = .
03 and 0 .
09. The nu-merical setup for these additional cases is the same as before butfor varying 𝛽 . Overall, the basic trend we find with 𝛽 = .
06 donot change even in such cases. Fig. 8 summarizes the evolution of 𝑁 c , b with different 𝛽 at 𝑍 =
0, 10 − , 10 − , and 10 − Z (cid:12) . Varying 𝛽 causes some scatter of the lines at each metallicity, but it looks minorcompared to the effects of varying 𝑍 . We only notice that in thecases of 𝑍 = − and 10 − Z (cid:12) the maximum of the clump number 𝑁 c , b , max is highest for 𝛽 = .
06, not for 𝛽 = .
09. We do not see asystematic trend that the more rapid initial cloud rotation leads to themore vigorous disc fragmentation.Fig. 9 displays the metallicity-dependence of the clump numberevolution with the rotation parameter fixed at 𝛽 = .
03 (upper panel)and 0 .
09 (lower panel). The figure shows the similar trend to the caseswith 𝛽 = .
06 (Fig. 7) as expected from Fig. 8. We only see that thescatter of the lines with 𝛽 = .
03 is smaller than the other cases. Thepeak values 𝑁 c , b , max at 𝑍 = − and 10 − Z (cid:12) are relatively smallerthan the counterparts with the higher 𝛽 . Nonetheless, the epochs of 𝑁 c , b , max at a given metallicity do not shift much even if varying 𝛽 . The line scatter with 𝛽 = .
09 looks more similar to that with 𝛽 = .
06 (Fig. 7).We finally highlight the particular case of 𝑍 = − Z (cid:12) and 𝛽 = .
09. In this case, the clump number evolution is similar to the othercases of the same metallicity for 𝛥𝑡 (cid:46) yr. However, the clumpnumber abruptly rises slightly after 𝛥𝑡 = yr and then declinesonly for this case. The corresponding line hence has the doublepeaks, at which 𝑁 c , b , max (cid:39)
10. Note that the horizontal axis is inlogarithmic scale in Figs. 8 and 9, and the declining timescales afterthe peaks are always ∼ yr. What happens here is the same asthat we have described in Section 3.1.2; the dust-induced vigorousdisc fragmentation followed by multiple merger events. This fact alsosuggests that disc fragmentation is a sporadic process, particularlyfor EMP cases. The evolution of the clump number should not bemonotonic but very variable in time. We also feature the second eventof the disc fragmentation in more detail in Section 4.1 later. We have mostly focused on the metallicity-dependence of the clumpnumber evolution caused by the disc fragmentation, inspired by re-cent studies on the primordial star formation (e.g., Susa 2019). Wehere further consider variations of the clump mass distribution withdifferent metallicities. However, recall that the total mass accreted byclumps until a given epoch significantly differs with different metal-licities (Fig. 3). Since this is not caused by the metallicity-dependentnature of the disc fragmentation, we instead consider the relative massdistribution of clumps against the most massive, primary object.We derive the relative clump mass distribution for each run as fol-lows. We examine all the snapshots taken every 10 years after the firstappearance of a self-gravitating clump, except those where there isonly a single object. For a given snapshot, we normalize each clump’smass 𝑚 c by the primary’s mass 𝑚 c , max to build up a probability dis- MNRAS , 1–14 (2021) K. Shima & T. Hosokawa N c , b = 0.03 Z = 0 Z = 10 ZZ = 10 ZZ = 10 Z t [yr]10 N c , b = 0.09 Z = 0 Z = 10 Z Z = 10 Z Z = 10 Z10 t [yr] Figure 9.
Same as Fig. 7 but for the cases with 𝛽 = .
03 (upper panel) and 𝛽 = .
09 (lower panel). tribution function 𝑝 as a function of the mass ratio 𝑥 ≡ 𝑚 c / 𝑚 c , max .We divide unity into 10 equal bins and evaluate 𝑝 ( 𝑥 ) per bin todraw a histogram. We derive such histograms for all the availablesnapshots and then take their average to derive ¯ 𝑝 ( 𝑥 ) . Note that wedo not include the primary clump for each 𝑝 ( 𝑥 ) histogram, becauseotherwise a prominent peak at 𝑥 = 𝑝 ( 𝑥 ) . In this sense, obtained 𝑝 ( 𝑥 ) represents the relative mass distri-bution of companion clumps associated with the primary one. Sincewe have shown that the effects of varying the rotation parameter 𝛽 is relatively minor than the metallicity-dependence (Section 3.3),we further average the cases of 𝛽 = .
03, 0 .
06 and 0 .
09 for eachmetallicity.Fig. 10 shows the results of our analyses. For instance, the panelfor 𝑍 = 𝑥 (cid:39) .
2. This indi-cates that a companion clump with 20 % mass of the primary clumpis the most typical during 𝛥𝑡 = × rm years, the duration of oursimulations. Nonetheless, the histogram has long tails, suggestingthat the companion clumps are somewhat evenly distributed in mass.The panels for 𝑍 = − and 10 − Z (cid:12) also show the wide distribu-tions with the peaks shifted to the lowest-mass bin. In these cases,the companion clumps tend to be less massive than the primary onethan the primordial case. Although not very clear, the histogram for 𝑍 = − Z (cid:12) is more skewed to the lower 𝑥 than for 𝑍 = − Z (cid:12) ,probably reflecting the dust-induced fragmentation is delayed untilthe primary clump grows to become relatively massive (Section 3.1).The panel for 𝑍 = − Z (cid:12) only shows the totally different featuresfrom the others; there are the double peaks at both ends. This isdue to the variation with different 𝛽 , which is remarkable only at 𝑍 = − Z (cid:12) . The histogram only for the case with 𝛽 = .
06 (black dashed line in the bottom right panel) shows a very high peak atthe lowest bin, the trend expected from the cases of 𝑍 = − and10 − Z (cid:12) . As described in Section 3.1.4, the single clump grows untilthe fragmentation starts in a late stage in this case. The other peaknear 𝑥 = 𝛽 = .
03 and 0 .
09, wherea binary system appears in an early stage at 𝛥𝑡 (cid:39) a few ×
100 yrafter the fragmentation. This early fragmentation is only weak andcreates a few self-gravitating clumps at maximum. The binary thencontinues to grow in mass, steadily accreting the gas through a cir-cumbinary disc. Further fragmentation does not occur for a while.The two clumps equally grow, and their mass ratio approaches tounity (Fig. 11, see also Chon & Hosokawa 2019). The more vigorousfragmentation starts slightly after 𝛥𝑡 = yr with 𝛽 = .
09 (Fig. 8),similar to the case with 𝛽 = . 𝑍 = − Z (cid:12) . Inthe lower-metallicity cases, more clumps temporarily appear becauseof the early vigorous fragmentation. However, they do not grow intothe binary system as formed at 𝑍 = − Z (cid:12) . Many clumps undergocomplex gravitational interactions that often cause mergers. A nextfragmentation episode occurs when the clump number settles down toa few, as seen for the case of 𝑍 = − Z (cid:12) and 𝛽 = .
09 (Section 3.3).Such a harsh environment prevents the steady mass growth of twinclumps in a binary. Since the efficient fragmentation begins beforethe end of simulations for 𝑍 = − Z (cid:12) , the binary system may beeventually disrupted or destroyed by gravitational interactions withother clumps. To consider the survival of the binary, we need furtherlong-term simulations that follow the later evolution, a task for futurestudies. As described in Section 3.1.2 and 3.1.3, the vigorous disc fragmenta-tion caused by efficient dust cooling occurs particularly for the casesof 𝑍 = − Z (cid:12) and 10 − Z (cid:12) . We here further investigate how thisprocess develops in more detail. Previous studies have already pro-vided a key concept that the gravitational instability of the spiral armsessentially represents the disc fragmentation. Takahashi et al. (2016)show that the linear stability analysis of a rotating ring or filamentwell describes the spiral-arm instability, with thorough comparisonsto their 2D simulations. Whereas Takahashi et al. (2016) specificallyconsider the fragmentation of a non-accreting massive protoplan-etary disc, Inoue & Yoshida (2018) show that the similar conceptis applicable for the galaxy-formation simulations in 3D. Inoue &Yoshida (2020) further apply the analyses to the cosmological sim-ulation of the primordial star formation performed by Greif et al.(2012). They show that the fragmentation of a rapidly accreting cir-cumstellar disc demonstrated by Greif et al. (2012) is essentially thesame process.Figs. 4 and 5 have already suggested that the temperature anddensity along the spiral arms roughly correspond to the values atthe local minima of the EOS curves at 𝑛 ∼ cm − (Fig. 1),where 𝛾 eff =
1. Fig. 12 now further displays featured images justafter the disc fragmentation for the same cases as in Figs. 4 and 5 inthe left and right columns. These snapshots are taken at 𝛥𝑡 =
80 yrand 340 yr for the former and latter cases, which correspond to theepochs of 𝑁 c , b , max in Fig. 7. The middle column shows the snapshotsat the much later stage of 𝛥𝑡 = 𝑍 = − Z (cid:12) and 𝛽 = .
09. This corresponds to the second peak of 𝑁 c , b , the MNRAS , 1–14 (2021) etal-poor disc fragmentation p ( x ) Z = 0 Z = 10 Z0.0 0.2 0.4 0.6 0.8 1.0 x = m c / m c,max p ( x ) Z = 10 Z 0.0 0.2 0.4 0.6 0.8 1.0 x = m c / m c,max Z = 10 Z Figure 10.
Normalized time-averaged mass distributions of self-gravitating companion clumps at different metallicities. The horizontal axis 𝑥 represents theratio of the companion mass 𝑚 c to the maximum clump mass 𝑚 c , max , and vertical axis represents the probability ¯ 𝑝 ( 𝑥 ) . In each panel, the filled histogramrepresents the probability distribution averaged over the cases with different initial cloud spin parameter 𝛽 = .
03, 0 .
06, and 0 .
09 at each metallicity. The blackdashed line also shows the distribution only with the case of 𝛽 = .
06. In the right bottom panel for 𝑍 = − Z (cid:12) , plotted is 0 . × ¯ 𝑝 ( 𝑥 ) so that it does not gobeyond the panel. second episode of the disc fragmentation that happens after the clumpnumber decreases down to a few (Fig. 8).We find the common feature for all the cases presented in Fig. 12.The geometrically thin spiral arms develop, and it is stretched tobecome the filamentary structure. At first glance, we recognize thatthe self-gravitating clumps are distributed well along the spiral armsand that the iso-density contours delineate them at 𝑛 = cm − .These facts strongly suggest that the spiral-arm instability is also theunderlying physics behind the dust-induced vigorous fragmentation.The dispersion relation derived by the linear stability analysis pre-dicts that the most unstable wavelength is (cid:39) ∼
10 clumps at once, and the violentmotion of these clumps prevents the growth of the spiral arms thatcover the whole disc. It is after the clump number settles down to afew that the large-scale spiral arms start to grow, leading to the nextfragmentation episode.
We have imposed the condition that the Jeans length 𝜆 J must beresolved by at least 32 cells, i.e., 𝜆 J / 𝛥𝑥 =
32, for the simulationspresented above. As described in Section 2.2, we also perform addi-tional simulations with varying the "Jeans criterion" as 𝜆 J / 𝛥𝑥 = 𝜆 J / 𝛥𝑥 =
64. We consider the case of 𝑍 = − Z (cid:12) and 𝛽 = .
06, where the vigorous fragmentation occurs in an early stage (Sec-tion 3.1.2), for such experimental runs.Fig. 13 shows the time evolution of the number of clumps withdifferent criteria. The upper panel shows the evolution of the self-gravitating clump number 𝑁 c , b . We see that the basic evolution doesnot change even if varying the Jeans criteria, though the peak value 𝑁 c , b , max is slightly reduced with 𝜆 J / 𝛥𝑥 =
16. There is a commontrend that 𝑁 c , b gradually decreases after taking peak values at 𝛥𝑡 ∼
100 yr because of many merger events. The number of survivingclumps eventually converges to almost the same value at 𝛥𝑡 ∼ yr.The figure shows that the evolution is particularly similar if 𝜆 J / 𝛥𝑥 >
32, and it well demonstrates the numerical convergence of our results.We also consider the role of checking whether a clump is grav-itationally bound or not (Section 2.4). The lower panel shows theevolution of the clump number without the binding check, 𝑁 c . Wesee the larger variations of the lines than in the upper panel. We detectalmost always more clump candidates with the more stringent cri-terion. This is not surprising because the smaller transient structureof the disc is resolved with the higher-resolution simulation realizedwith the stringent Jeans criterion. In particular, the number of suchtransient structure is more numerous by an order of magnitude thanthe self-gravitating clumps with 𝜆 J / 𝛥𝑥 =
64. This suggests how crit-ical the binding check is for counting the clump number. One maysignificantly overestimate the number if misidentifying the transientstructure.
MNRAS , 1–14 (2021) K. Shima & T. Hosokawa y [ A U ] Z = 10 Z= 0.03 t = 1000 yr100 50 0 50 100 x [AU]10050050100 y [ A U ] Z = 10 Z= 0.09 t = 1000 yr Figure 11.
Comparing images of the number density maps with the differentinitial cloud rotation (top: 𝛽 = .
03, and bottom: 𝛽 = .
09) for 𝑍 = − Z (cid:12) .The snapshots are taken at the epoch of 𝛥 𝑡 = yr for both cases. The massratio of the clumps is almost unity for these cases, whereas it is small with 𝛽 = .
06 at the same epoch (cf. the top right panel in Fig. 6). The color scaleis the same as in Fig. 6. The mass of each clump is also described in the unitof M (cid:12) . We have studied the gravitational fragmentation of accreting circum-stellar discs with various metallicities, by performing a suite of 3Dhydrodynamic simulations using the adaptive mesh refinement code
Enzo . We model the metallicity-dependent EOS of the gas using pre-calculated barotropic EOS at 𝑍 =
0, 10 − , 10 − , and 10 − Z (cid:12) . Asimulation run begins with an idealized rotating cloud characterizedby the spin parameter 𝛽 . We have followed the evolution from theearly collapse to the subsequent accretion stage. In particular, wehave investigated the long-term evolution in the late accretion stagefor 𝛥𝑡 = × yr, which is longer by an order of magnitude than inprevious relevant studies. We have further studied the effects of vary-ing the cloud rotation parameter 𝛽 and the so-called Jeans criteriaperforming the additional simulations. Our findings are summarizedas follows. Our simulations show that the disc fragmentation occurs for all theexamined cases, regardless of the metallicity. However, the resultingevolution shows the clear metallicity dependence. We have paid spe-cial attention to the evolution of the number of self-gravitating clumpsformed by the fragmentation, 𝑁 c , b . The primordial case shows thatthe fragmentation, often followed by clump mergers, steadily contin-ues until the end of the simulation. The clump number stays almostconstant at a few during that, while it does not monotonically in-crease as suggested by Susa (2019). In contrast, the evolution of 𝑁 c , b becomes variable in time with a tiny amount of metals. The vigor-ous fragmentation caused by efficient dust cooling occurs in an earlystage at 𝑍 = − and 10 − Z (cid:12) , as predicted by Tanaka & Omukai(2014) using 1D semi-analytic modeling. The clump number tem-porarily rises to 𝑁 c , b ∼
10, but it continuously decreases as manyclumps merge away after that. The vigorous fragmentation tends tooccur later with the higher 𝑍 , reflecting that the dust-induced frag-mentation is most efficient at the lower density. At 𝑍 = − Z (cid:12) , theclump number is smallest until the efficient fragmentation eventu-ally starts at 𝛥𝑡 ∼ yr. In all the cases, the clump number settlesdown to a few by the specific epoch of 𝛥𝑡 = × years after suchvery different evolution. The above picture does not change even ifvarying the cloud’s initial rotation parameter 𝛽 .We have also analyzed the simulation data to investigate the time-averaged relative mass distribution of the clumps for all the examinedcases. The mass distribution also shows a systematic trend; the com-panion clumps become relatively less massive than the primary ormost massive one with increasing 𝑍 . This reflects the metallicity-dependent evolution described above, i.e., the primary clump hasmore time to accrete the gas until the vigorous fragmentation startsat the higher 𝑍 . On top of this trend, there is another striking differ-ence in the mass distribution at 𝑍 = − Z (cid:12) ; the other peak nearthe high-mass end, representing a binary system with twin clumpssubstantially more massive than the others. We have shown that sucha characteristic system grows through steady accretion from a cir-cumbinary disc, during which the system is not disturbed by otherclumps. It seems that this hardly occurs at 𝑍 = − and 10 − Z (cid:12) because many clumps produced by the early dust-induced fragmen-tation continue to interact with each other violently. Although ourcurrent simulations only follow the initial 2 × years of the pro-tostellar accretion stage, the result suggests how massive and equal-mass binaries form in low-metallicity environments. ACKNOWLEDGEMENTS
We thank Kazuyuki Sugimura, Ryoki Matsukoba, Gen Chiaki,Shigeki Inoue, Naoki Yoshida, Kazuyuki Omukai, and Hajime Susafor useful discussion and comments. This work is financially sup-ported by the Grants-in-Aid for Basic Research by the Ministry ofEducation, Science and Culture of Japan (17H06360, 19H01934:T.H.). The numerical simulations were performed on the Cray XC50(Aterui II) at the Center for Computational Astrophysics (CfCA) ofNational Astronomical Observatory of Japan. The simulation resultsare analyzed using the visualization toolkit for astrophysical data YT(Turk et al. 2011).
DATA AVAILABILITY
The data underlying this article will be shared on reasonable requestto the corresponding author.
MNRAS000
MNRAS000 , 1–14 (2021) etal-poor disc fragmentation y [ A U ] Z = 10 Z= 0.06 t = 80yr Z = 10 Z= 0.09 t = 1310yr Z = 10 Z= 0.06 t = 340yr100 50 0 50 100 x [AU]10050050100 y [ A U ]
100 50 0 50 100 x [AU] 100 50 0 50 100 x [AU] 10 nu m b e r d e n s i t y [ c m ] t e m p e r a t u r e [ K ] Figure 12.
Featured images of the density-weighted projection maps of the number density (top) and temperature (bottom) at the epochs when the number ofclumps remarkably increases due to efficient fragmentation. The left, middle, and right columns of panels represent different runs for ( 𝑍 , 𝛽 ) = ( − Z (cid:12) , . ) , ( − Z (cid:12) , . ) , and ( − Z (cid:12) , . ) , respectively. The elapsed time 𝛥 𝑡 for each snapshot is also presented in each top panel. The white lines in the densitymaps represent the iso-density contours at 𝑛 = cm − . The positions of self-gravitating clumps are denoted with white crosses in the temperature maps. REFERENCES
Bate M. R., Bonnell I. A., Price N. M., 1995, MNRAS, 277, 362Bonnor W. B., 1956, MNRAS, 116, 351Bromm V., Ferrara A., Coppi P. S., Larson R. B., 2001, MNRAS, 328, 969Bromm V., Coppi P. S., Larson R. B., 2002, ApJ, 564, 23Bryan G. L., et al., 2014, ApJS, 211, 19Chiaki G., Yoshida N., 2020, arXiv e-prints, p. arXiv:2008.06107Chiaki G., Marassi S., Nozawa T., Yoshida N., Schneider R., Omukai K.,Limongi M., Chieffi A., 2015, MNRAS, 446, 2659Chiaki G., Yoshida N., Hirano S., 2016, MNRAS, 463, 2781Chon S., Hosokawa T., 2019, MNRAS, 488, 2658Clark P. C., Glover S. C. O., Klessen R. S., 2008, ApJ, 672, 757Clark P. C., Glover S. C. O., Smith R. J., Greif T. H., Klessen R. S., BrommV., 2011, Science, 331, 1040Dopcke G., Glover S. C. O., Clark P. C., Klessen R. S., 2011, ApJ, 729, L3Dopcke G., Glover S. C. O., Clark P. C., Klessen R. S., 2013, ApJ, 766, 103Ebert R., 1955, Z. Astrophys., 36, 222Federrath C., Banerjee R., Clark P. C., Klessen R. S., 2010, ApJ, 713, 269Federrath C., Sur S., Schleicher D. R. G., Banerjee R., Klessen R. S., 2011,ApJ, 731, 62Fukushima H., Hosokawa T., Chiaki G., Omukai K., Yoshida N., Kuiper R.,2020, MNRAS, 497, 829Gammie C. F., 2001, ApJ, 553, 174Greif T. H., 2015, Computational Astrophysics and Cosmology, 2, 3Greif T. H., Bromm V., Clark P. C., Glover S. C. O., Smith R. J., KlessenR. S., Yoshida N., Springel V., 2012, MNRAS, 424, 399Hirano S., Bromm V., 2017, MNRAS, 470, 898Hirano S., Hosokawa T., Yoshida N., Umeda H., Omukai K., Chiaki G., YorkeH. W., 2014, ApJ, 781, 60Hosokawa T., Omukai K., Yoshida N., Yorke H. W., 2011, Science, 334, 1250 Hosokawa T., Hirano S., Kuiper R., Yorke H. W., Omukai K., Yoshida N.,2016, ApJ, 824, 119Hubber D. A., Walch S., Whitworth A. P., 2013, MNRAS, 430, 3261Ilee J. D., Cyganowski C. J., Brogan C. L., Hunter T. R., Forgan D. H.,Haworth T. J., Clarke C. J., Harries T. J., 2018, ApJ, 869, L24Inoue S., Yoshida N., 2018, MNRAS, 474, 3466Inoue S., Yoshida N., 2020, MNRAS, 491, L24Jappsen A. K., Klessen R. S., Larson R. B., Li Y., Mac Low M. M., 2005,A&A, 435, 611Jappsen A.-K., Klessen R. S., Glover S. C. O., Mac Low M.-M., 2009, ApJ,696, 1065Kimura K., Hosokawa T., Sugimura K., 2020, arXiv e-prints, p.arXiv:2012.01452Kratter K., Lodato G., 2016, ARA&A, 54, 271Krumholz M. R., McKee C. F., Klein R. I., 2004, ApJ, 611, 399Larson R. B., 1969, MNRAS, 145, 271Larson R. B., 1985, MNRAS, 214, 379Larson R. B., 2005, MNRAS, 359, 211Machida M. N., Doi K., 2013, MNRAS, 435, 3283Machida M. N., Nakamura T., 2015, MNRAS, 448, 1405Machida M. N., Omukai K., Matsumoto T., Inutsuka S.-i., 2008, ApJ, 677,813Machida M. N., Inutsuka S.-i., Matsumoto T., 2011, ApJ, 729, 42Matsukoba R., Vorobyov E. I., Sugimura K., Chon S., Hosokawa T., OmukaiK., 2021, MNRAS, 500, 4126McKee C. F., Tan J. C., 2008, ApJ, 681, 771Oliva G. A., Kuiper R., 2020, A&A, 644, A41Omukai K., 2000, ApJ, 534, 809Omukai K., Nishi R., 1998, ApJ, 508, 141Omukai K., Tsuribe T., Schneider R., Ferrara A., 2005, ApJ, 626, 627MNRAS , 1–14 (2021) K. Shima & T. Hosokawa N c , b w/ binding check J / x = 16 J / x = 32 J / x = 6410 t [yr]10 N c w/o binding check10 t [yr] Figure 13.
Time evolution of the number of clumps found by the finderwithout (upper) and with (lower) the gravitational binding energy check inthe case of 𝑍 = − Z (cid:12) and 𝛽 = .
06. In both panels, the sold line representsour standard case with the Jeans criterion 𝜆 J / 𝛥 𝑥 =
32, i.e., the Jeans length isalways resolved by at least 32 cells. The dot-dashed and dashed lines representthe less and more stringent criteria, 𝜆 J / 𝛥 𝑥 =
16 and 64, respectively. Thegreen line terminates at 𝛥 𝑡 (cid:39)
800 yr since we stop the corresponding high-resolution simulation at this point because of the heavy computational cost.Regan J. A., Downes T. P., 2018, MNRAS, 475, 4636Safranek-Shrader C., Milosavljević M., Bromm V., 2014, MNRAS, 438, 1669Saigo K., Matsumoto T., Umemura M., 2004, ApJ, 615, L65Schneider R., Omukai K., Inoue A. K., Ferrara A., 2006, MNRAS, 369, 1437Schneider R., Omukai K., Bianchi S., Valiante R., 2012, MNRAS, 419, 1566Sharda P., Krumholz M. R., Federrath C., 2019, MNRAS, 490, 513Shima K., Tasker E. J., Federrath C., Habe A., 2018, PASJ, 70, S54Smith B., Sigurdsson S., Abel T., 2008, MNRAS, 385, 1443Smith B. D., Turk M. J., Sigurdsson S., O’Shea B. W., Norman M. L., 2009,ApJ, 691, 441Smith R. J., Glover S. C. O., Clark P. C., Greif T., Klessen R. S., 2011,MNRAS, 414, 3633Stacy A., Bromm V., 2013, MNRAS, 433, 1094Stacy A., Greif T. H., Bromm V., 2010, MNRAS, 403, 45Stacy A., Bromm V., Lee A. T., 2016, MNRAS, 462, 1307Stone J. M., Norman M. L., 1992a, ApJS, 80, 753Stone J. M., Norman M. L., 1992b, ApJS, 80, 791Sugimura K., Mizuno Y., Matsumoto T., Omukai K., 2017, MNRAS, 469,4022Sugimura K., Matsumoto T., Hosokawa T., Hirano S., Omukai K., 2020, ApJ,892, L14Susa H., 2013, ApJ, 773, 185Susa H., 2019, ApJ, 877, 99Susa H., Doi K., Omukai K., 2015, ApJ, 801, 13Takahashi S. Z., Tsukamoto Y., Inutsuka S., 2016, MNRAS, 458, 3597Takahira K., Tasker E. J., Habe A., 2014, ApJ, 792, 63 Tan J. C., McKee C. F., 2004, ApJ, 603, 383Tanaka K. E. I., Omukai K., 2014, MNRAS, 439, 1884Tobin J. J., et al., 2016, Nature, 538, 483Tsukamoto Y., Machida M. N., Inutsuka S.-i., 2013, MNRAS, 436, 1667Tsuribe T., Omukai K., 2006, ApJ, 642, L61Turk M. J., Smith B. D., Oishi J. S., Skory S., Skillman S. W., Abel T., NormanM. L., 2011, ApJS, 192, 9Vorobyov E. I., Basu S., 2010, ApJ, 719, 1896Vorobyov E. I., DeSouza A. L., Basu S., 2013, ApJ, 768, 131Yahil A., 1983, ApJ, 265, 1047Yoshida N., Omukai K., Hernquist L., 2008, Science, 321, 669This paper has been typeset from a TEX/L A TEX file prepared by the author.MNRAS000