Formation of Massive and Wide First-star Binaries in Radiation Hydrodynamics Simulations

We study the formation of Pop III stars by performing radiation hydrodynamics simulations for three different initial clouds extracted from cosmological hydrodynamics simulations. Starting from the cloud collapse stage, we follow the growth of protostars by accretion for $\sim 10^5$ yr until the radiative feedback from the protostars suppresses the accretion and the stellar properties are nearly fixed. We find that the Pop III stars form in massive and wide binaries/small-multiple stellar systems, with masses $>30\,M_\odot$ and separations $>2000$ au. We also find that the properties of the final stellar system correlate with those of the initial clouds: the total mass increases with the cloud-scale accretion rate, and the angular momentum of the binary orbit matches that of the initial cloud. While the total mass of the system in our simulations is consistent with our previous single-star formation simulations, individual masses are lower due to mass sharing, suggesting potential modification in the extent of feedback from Pop III stars in the subsequent evolution of the Universe. We also identify such systems as mini-binaries embedded in a wider outer multiple-star system, which could evolve into progenitors for observed gravitational wave events.

Previous numerical studies that followed the growth of embryonic Pop III stars up to the end of the accretion (e.g., Hosokawa et al. 2011;Susa 2013;Stacy & Bromm 2013;Susa et al. 2014;Hirano et al. 2014Hirano et al. , 2015;;Stacy et al. 2016;H16) demonstrate that the Pop III stars tend to be more massive than predicted from the present-day Salpeter-like initial mass function (IMF).This allows Pop III stars to play crucial roles in driving the subsequent cosmic evolution through various processes, such as the reionization of the intergalactic medium by stellar radiation (e.g., Schaerer 2002;Wise et al. 2012;Ricotti et al. 2016), the metal enrichment of the pristine gas by supernovae (e.g., Woosley et al. 2002;Chiaki et al. 2018;Abe et al. 2021), and the seeding of supermassive black hole (BH) progenitors (e.g., Alvarez et al. 2009;Jeon et al. 2012).
Observations show that massive stars are predominantly in multiple systems in the nearby Universe (e.g., Duchêne & Kraus 2013).If this is also the case in the early Universe, massive Pop III stars would also be in multiples.In this case, mass reservoir in the parental cloud will be shared among forming stars and the individual stellar masses will decrease.Such Pop III binaries are attracting attention as they may later evolve into binary black holes detectable by gravitational waves (e.g., Kinugawa et al. 2014;Hartwig et al. 2016;Abbott et al. 2016;Tanikawa et al. 2021;Liu & Bromm 2021, but see also Belczynski et al. 2017) or the first X-ray binaries, which are important heating sources of the intergalactic medium and should be constrained by future 21cm line observations (e.g., Mirabel et al. 2011;Dewdney et al. 2009).The X-ray background can also change the mode of the Pop III star formation (e.g., Ricotti 2016;Park et al. 2021aPark et al. ,b, 2023)).
On the theoretical side, seminal works found that Pop III stars in fact end up in multiple-star systems, by way of numerical simulations using particle-based codes (Susa 2013;Stacy & Bromm 2013;Susa et al. 2014;Stacy et al. 2016).In those calculations, however, protostellar feedback may have been underestimated because the density-dependent resolution of the particle-based codes was insufficient to directly follow the propagation of the ionizing radiation in the low-density polar regions around each protostar (Susa 2013).
Later, in Sugimura et al. (2020) (hereafter Paper I), we confirmed that Pop III stars do form in a multiple-stellar system by performing a simulation with a newly developed radiation hydrodynamics (RHD) code, SFUMATO-RT, which employs the adaptivemesh-refinement (AMR; Berger & Colella 1989;Berger & Oliger 1984) and the adaptive-ray-tracing (ART; Abel & Wandelt 2002) techniques to accurately follow the ionizing radiation from multiple protostars.In that work, the formation of multiple protostars by gas fragmentation and their long-term evolution are followed until the radiation feedback terminates the accretion, thereby fixing the stellar properties.
Quite recently, other groups also performed similar simulations using other AMR codes (Latif et al. 2022;Park et al. 2023) and a moving-mesh code (Jaura et al. 2022) with different radiation transfer modules.While all groups agree on the formation of multiple stars, there is a discrepancy on the effectiveness of radiation feedback.Jaura et al. (2022) argued that the trapping of ionizing radiation near protostars significantly weakens the radiative feedback unlike in the other simulations where efficient feedback is observed.The cause of this discrepancy is still unclear and further investigation is needed.
In Paper I, our simulation was limited only to a single case, despite the known diversity of the birth environments (e.g., Hirano et al. 2014).We could not discuss such statistics as the relation between the final stellar mass and the natal cloud properties, found in 2D simulations (Hirano et al. 2014).Furthermore, in Paper I, we just focused on presenting the overall evolution and the nature of resulting stellar system, leaving analysis of detailed physical process to future work.
In this paper, we simulate Pop III star formation in three different clouds extracted from cosmological simulations, using the same code as in Paper I.Among the three clouds studied, two are new cases while the one is the same as that already presented in Paper I. Based on the results, we examine the relation between the properties of the final stellar systems, such as the total mass and binary separation, and those of the natal clouds.We also identify and describe in detail some interesting phenomena that can play important roles in Pop III star formation.
The rest of this paper is organized as follows.In Section 2, we describe our numerical methods in detail.In Section 3, we present our simulation results, from which we then investigate relations between the stellar system and natal clouds in Section 4. In Section 5, we discuss the role of radiative feedback and the numerical resolution effects.We conclude the paper in Section 6.The appendices provide a detailed description of simulation methods and a supplementary analysis of circum-stellar disks in our simulations.

SFUMATO-RT
Our radiation hydrodynamics code SFUMATO-RT has been developed to follow the Pop III star formation and was first used in Paper I. The code is an extended version of an AMR self-gravitational magnetohydrodynamics (MHD) code SFUMATO (Matsumoto 2007;Matsumoto et al. 2015), with addition of new modules to solve the non-equilibrium chemistry and thermodynamics of primordial gas under the influence of radiation from multiple protostars (see Sadanari et al. 2021Sadanari et al. , 2023;;Kimura et al. 2023, for other extensions).

Hydrodynamics
We use SFUMATO's modules to solve hydrodynamics with self-gravity and sink particles that represent accreting protostars.We do not use the MHD module since we ignore magnetic fields in this work.
The basic equations of hydrodynamics are where ρ is the mass density, v the velocity, P the pressure, g the gravitational acceleration, Γ the heating rate, Λ the cooling rate, and γ the adiabatic exponent.The mass density and pressure are related to the number density of hydrogen atoms n H as ρ = n H µ m p and P = n H i y(i) k B T , respectively, with µ the mean molecular weight per hydrogen atom, T the temperature, m p the proton mass, k B the Boltzmann constant, and y(i) the chemical abundances defined as the abundance ratio of species i to the hydrogen nuclei.We determine Γ, Λ, γ, and µ from T and y(i) (see, e.g., Omukai & Nishi 1998).The chemical and thermal model is described in Section 2.1.2.In this work, we use the hydrodynamical scheme with second-order accuracy in space and time.
We consider both the self-gravity of the gas and the gravity by the sink particles.The total gravitational acceleration is given by where Φ is the gravitational potential of the gas and g sink,i is the gravitational acceleration due to the i-th sink particle.Using the multigrid solver, we obtain Φ from the Poisson equation, with Newton's gravitational constant G.As for g sink,i , we directly evaluate Newton's inverse-square law, where x and x i are the positions of the gas and the i-th particle, respectively.Using the sink particle technique, we mask the neighborhoods of protostars with extremely short timescales to perform long-term simulations until the end of the accretion phase.Below, we briefly describe the implementation of sink particles in SFUMATO and refer the reader to Matsumoto et al. (2015) for details.
We create a new sink particle in a cell that satisfies the following conditions (Federrath et al. 2010): (i) the density is higher than the sink threshold density n sink .
(ii) the cell is a local minimum of gravitational potential Φ.
(iii) all the eigenvalues of the symmetric parts of the velocity gradient tensor ∇ i v j are negative.
(iv) the total energy of the gas within the sink radius r sink is negative (E thermal + E kin + E grav < 0).
Around each sink particle, we define a virtual sphere with the radius r sink (also called the sink sphere), which is used for the following purposes.First, the sink particle accretes the gas from the cells inside the sink sphere if the density exceeds n sink until the excess disappears.Second, the sink gravity is weakened inside the sink sphere.Finally, two sink particles merge when their sink spheres overlap.In our fiducial runs, we set r sink = 64 au (16 times the minimum cell size) and n sink = 2 × 10 11 cm −3 , as described in Section 2.3.
The AMR technique allows us to follow fine structures near multiple protostars at a relatively low computational cost.We refine the cells so that the local Jeans length is resolved with at least N J cells.In this work, we set N J = 16, which is much higher than the so-called Truelove condition N J ≥ 4 (Truelove et al. 1997) for avoiding artificial fragmentation.Although it has recently been claimed that higher resolution (N J ≳ 30) is needed to follow the turbulence amplification associated with gravitational collapse (Federrath et al. 2011;Higashi et al. 2021Higashi et al. , 2022)), we adopt the above value to perform our simulations until the end of the accretion phase.

Non-equilibrium chemistry and thermodynamics
One of major new additions of SFUMATO-RT to SFUMATO is a module for the non-equilibrium chemistry and thermodynamics of the primordial gas, which is essential to simulate Pop III star formation with realistic thermal evolution.We consider the chemical reactions among six species, H + , H, H 2 , H − , H + 2 , and e − , assuming that all the helium is neutral.To update y(i) and T consistently, we solve the coupled equations for the chemical and temperature evolution, adopting an implicit method for numerical stability.Our chemical and thermal models are summarized in Appendix A.
In short, we consider all the chemical reactions and cooling and heating processes that are relevant in the density range of n H < 10 13 cm −3 , as in H16.Major chemical reactions include the H − -channel H 2 formation, the 3-body H 2 formation, the collisional dissociation and the photodissociation of H 2 , the collisional ionization and the photoionization of H, the radiative recombination of H + .Major cooling and heating processes include the H 2 line cooling with line-trapping effect, the H Lyα cooling, the H − free-bound cooling, the H free-free cooling, the H free-bound cooling, the chemical cooling/heating, and the adiabatic compression heating and expansion cooling.
Inside sink spheres, we solve the non-equilibrium chemistry and thermodynamics as in the normal cells, but with neglected coupling to the radiation to avoid artificial enhancement of radiative feedback.We adopt this approach to prevent potential numerical artifacts arising from discontinuities between the interior and exterior of sink spheres.The radiation transfer model is described in the next section.

Radiation
Another feature of SFUMATO-RT is addition of a radiation module to handle the radiation from multiple sources.We here consider three types of radiation: extreme-ultraviolet radiation (EUV; hν > 13.6 eV), which photoionizes H; far-ultraviolet radiation (FUV; 11.2 eV < hν < 13.6 eV), which photodissociates H 2 ; and near-infrared radiation (NIR; hν < 11.2 eV), which photodetaches H − .We solve the radiation transfer for the EUV and FUV, while the NIR is assumed to be optically thin, which is a good approximation in the accreting envelope (see Hosokawa et al. 2011).Below, we describe the protostellar model, the radiation transfer method, the sink-scale treatments, and the coupling with chemistry and thermodynamics in this order.
Protostellar model -We assign the radiative properties of Pop III protostars by tabulated results of onedimensional (proto)stellar evolution calculations under constant accretion rates (Hosokawa & Omukai 2009;Hosokawa et al. 2010).The table gives the luminosity L * and stellar radius R * for a given set of the stellar mass M and accretion rate Ṁ .From L * and R * , we derive the effective temperature T eff , EUV emissivity S EUV , and FUV emissivity S FUV , assuming a blackbody spectrum (see Appendix B for more details).We also obtain the optically-thin H − photodetachment rate as a function of the distance from the protostar.
Following H16, we average the accretion rate over 300 years to account for the timescale of gas transport through unresolved parts of the accretion disk.We only consider radiation from protostars more massive than M rad,min = 5 M ⊙ to save computational cost since the UV emissivities before the onset of the Kelvin-Helmholtz contraction are extremely small (see Appendix B).
Our protostellar radiation model neglects the dependence of stellar properties on the accretion history.H16 found a case in which the intermittent plunge of fragments causes accretion bursts, and thus the accretion history matters indeed.In our simulations, however, we expect the dependence on the accretion history to be less important because the gas is mostly supplied to protostars by continuous gas accretion from the circum-stellar disks, with modest modulation due to binary interactions, except in early phases when radiation feedback is negligible.
Radiative transfer method -To solve the transfer of EUV and FUV radiation, we trace rays from each source by way of the ART method (Abel & Wandelt 2002;Krumholz et al. 2007;Wise & Abel 2011;Rosen et al. 2017;Kim et al. 2017).Below, we briefly explain our radiative transfer method and refer the readers to the literature above for details about the ART method.We present our specific implementation to SFUMATO-RT in Appendix C.
For the EUV radiation, we trace the direct photons from each protostar considering the attenuation due to photoionization, with the diffuse recombination photons treated by the on-the-spot approximation, following H16.At a distance r from a radiation source, the optical depth is obtained by the integration along a ray as with σ pi the effective photoionization cross-section depending on the source's spectrum.We omit the attenuation within the sink sphere in Eq. ( 8) by setting the lower bound of the integration to r sink , and separately consider the sink-scale attenuation, as described below.
With τ EUV along each ray, we evaluate the EUV photoionization rate in each cell as and the EUV photo-heating rate as where ϵ EUV = ⟨hν − 13.6eV⟩ is the mean energy deposited in the gas per photoionization.The bracket ⟨ ⟩ cell denotes the average over the rays penetrating the cell (e.g., Wise & Abel 2011).If rays from multiple sources reach the cell, we take the sum of each source's contribution.
For the FUV radiation, instead of the optical depth, we evaluate the H 2 column density along a ray, to model the self-shielding effect.We then estimate the FUV photodissociation rate as with the effective photodissociation cross-section σ pi and the self-shielding factor f shield .
In the ART framework, to achieve desired directional resolution, we split rays hierarchically using the HEALPix library (Górski et al. 2005).At each ray level l ray , the spherical 4π solid angle is sampled by N ray = 12 × 4 lray rays that cover equal-area pixels with 4π/N ray .We start ray tracing from each sink particle by injecting N ray,ini = 12 × 4 lray,ini rays at the initial ray level l ray,ini = 3.While tracing the rays, we split a parent ray into four daughter rays at one higher HEALPix level to ensure that at least N ray,min rays pass through each cell surface.In this study, we use N ray,min = 3, following the resolution study of Krumholz et al. (2007).We have examined the effect of this choice by performing a test run with N ray,min = 5 and have confirmed that its effect on the total mass evolution is insignificant, at least until middle phases of star formation.We randomly rotate the HEALPix orientation at each ART step to mitigate artifacts due to insufficient discretization (Krumholz et al. 2007).
Sink-scale treatments -We need separate treatments for radiation rays inside sink spheres, where the gas distribution is not resolved and can be artificial.For example, the inner geometrically thin part of an accretion disk in reality may be artificially thick in simulations due to the softening of the protostar's gravity inside the sink sphere.
To account for the shielding of central radiation by unresolved parts of a disk inside a sink sphere, we allow only rays with elevation angles greater than the disk thickness to cross the sink sphere.Here, we measure the disk thickness at the sink radius, assuming that the aspect ratio is an increasing function of the radius.1In practice, when a ray crosses a sink sphere, we check whether the density is higher than the disk threshold density n disk = 10 −2 n sink .If this is the case, assuming that the ray is traveling in the direction shielded by the disk, we terminate the ray before it has any effect on the surrounding gas.We have confirmed that the effect of radiative feedback is not sensitive to the particular choice of n disk , by comparing axisymmetric test simulations of accreting protostars with fiducial n disk and 10 times larger n disk .
In addition to the disk shielding, we model the EUV absorption in the polar directions following H16.For each ray emitted from a protostar, we measure the the density of the cell where the ray crosses the sink sphere and check whether the Strömgren radius r strm for a homogeneous medium with this density (the temperature is assumed to be T = 4 × 10 4 K, which is a typical value for HII regions around Pop III stars) is smaller than the sink radius.If this is the case, we assume that the EUV is consumed inside the sink sphere and set a very large τ EUV (effectively no EUV effect remains).Of course, the assumption of a homogeneous medium inside the sink sphere is not very accurate: the density along the ray increases inward in the case of a spherical free-falling flow as ∝ r −3/2 ; conversely, the density may decrease in the case that the polar regions are cleared due to the centrifugal force.
Alternatively, Jaura et al. (2022) distributed photons at the center of sink spheres and solve the radiation transfer inside those spheres.As a result, the effect of radiative feedback is strongly suppressed in their simulations, as all photons are absorbed by the thick disks inside the sink spheres.At present, the precise gas distributions surrounding protostars are unknown, introducing uncertainties into the sink-scale radiation model used in star formation simulations.We expect that highresolution RHD simulations, capable of resolving the star-disk interfaces and the inner part of accretion disks (see Kimura et al. 2023), are valuable in addressing this matter.
Coupling with chemistry -In each ART step, we trace rays on the fixed background of the density and chemical composition.We then update the chemistry and thermodynamics under the fixed radiation field reconstructed from the result of the previous ray tracing.We perform the ART plus chemistry/thermodynamics subcycles n sub times per hydrodynamics step.In this study, we choose n sub = 2 so that the maximum propagation speed of the dissociation and ionization fronts in the subcycles (two cells per one hydrodynamics step) is faster than that associated with the hydrodynamic flows (one cell per one hydrodynamics step).As a result, we expect that the fronts can (in principle) reach the positions determined by the ionization/recombination and dissociation/formation equilibria, respectively.We extract primordial star-forming clouds from cosmological 3D SPH/N-body simulations (Hirano et al. 2014(Hirano et al. , 2015) ) as the initial conditions for our 3D AMR RHD star-formation simulations.Specifically, we remap particle-based simulation data of primordial clouds early in the collapse phase (when the central density reaches 10 6 cm −3 ) to Cartesian grids to generate the initial conditions.

Initial conditions
Table 1 summarizes the properties of the three clouds studied in this paper.These clouds have a range of initial cloud-scale accretion rate Ṁcloud spanning from 10 −3 to 10 −2 M ⊙ /yr.We give names to these clouds based on their respective cloud-scale accretion rates: High-, Intermediate-, and Low-Ṁ clouds.The Intermediate-Ṁ cloud is the same one studied in Paper I. Our selection of clouds with different Ṁcloud enables us to study a wide variety of Pop III star forming environments, as Ṁcloud is known to correlate with the final stellar mass (Hirano et al. 2014, H16).
After the onset of our re-simulation, the pre-stellar collapse continues until the first protostar, i.e., sink particle, appears around the cloud center.Here, we show the state of the clouds just before the first protostar formation by the 3D rendering of the central cores in Fig. 1 and the 1D radial profiles of the entire clouds in Fig. 2.
Fig. 1 depicts the diverse morphology of central cores in the three clouds.The cores in the Low-and Intermediate-Ṁ cases are rotating and have disk-like shapes, with the former having a noticeable bar-spiral structure.In contrast, the core in the High-Ṁ case is turbulent and filamentary in shape.Those cores can be considered as the initial states of Pop III star formation via accretion, as described in Section 3.
As seen in the density and temperature profiles shown in the top two panels of Fig. 2, the clouds undergo socalled the runaway collapse, with a slightly increasing temperature characterized by an effective polytropic index γ eff ≈ 1.1 (Omukai & Nishi 1998).The third panel 10 6 10 8 10 10 The radial profiles of physical quantities just before the first protostar formation for the three clouds (Low-Ṁ , green dotted; Intermediate-Ṁ , orange dashed; and High-Ṁ , blue solid).From top to bottom, the density nH, temperature T , inflow rate Ṁ , degree of rotational support f Kepler , and turbulent Mach number M turb are plotted.presents the inflow rates, Ṁ , at a given radius, which are overall consistent with the values of Ṁcloud (Table 1).Note, however, that those two quantities do not exactly match at either radius because the measurements are taken at different epochs while Ṁ gradually increases over time (Hirano et al. 2015).The bottom two panels show the degree of rotational sup-port f kep = Ω/Ω Kepler and the turbulent Mach number M turb = v turb /c s , where Ω Kepler is the Keplerian angular velocity and v turb is the turbulent velocity.The High-Ṁ case exhibits stronger turbulence than the two other cases, with a lower degree of rotational support.The initial states of the clouds before the accretion phase are closely linked to the properties of final Pop III systems, as discussed in Section 4.

Simulation setup
We perform three runs, namely, High, Int, and Low runs, corresponding to the High-Ṁ , Intermediate-Ṁ , and Low-Ṁ clouds in Section 2.2, respectively.Table 2 summarizes our simulation setup and results.
We set the side length of our simulation box to 2.6 or 5.2 pc to encompass the collapsing region of the cloud with sufficient margin.The minimum cell size at the highest AMR level is ∆x min = 4 au, with the sink radius of r sink = 64 au, equivalent to 16 times ∆x min .Such a large number of cells per sink radius is necessary to resolve geometrically-thin accretion disks with the thickness of ∼ 10 au around the sink radius: inadequate resolution would artificially align the disk to one of the Cartesian axes since the gas can be advected only through the cell surfaces in grid simulations.We adopt a sink threshold density of n sink = 2×10 11 cm −3 , ensuring that the corresponding Jeans length matches the sink diameter of 2 r sink (with T = 10 3 K, representing a typical temperature of neutral gas at density n sink ).We terminate our simulations at ∆t = 80, 000 yr or 100, 000 yr, when accretion is quenched by radiation feedback and the total mass is almost fixed.

NUMERICAL RESULTS
In this section, we present our simulation results.We first describe the overall evolution in Section 3.1, and then examine noteworthy interesting phenomena in detail in Section 3.2.

Overall evolution
First, we provide an overview of the star formation process in the three clouds simulated in this paper.To visually depict the evolution, we show the face-on slice density in Fig. 3 and the edge-on slice temperature in Fig. 4 for the Low-(top), Intermediate-(middle), and High-Ṁ (bottom) runs at different evolutionary stages.Specifically, we plot the snapshots at ∆t = 3, 000 yr (left), 10, 000 yr (middle), and 30, 000 yr (right), with ∆t being the time since the first protostar formation.In all three cases, multiple protostars are formed as a result of the initial disk fragmentation (Fig. 3, left), in agreement with previous simulations of the Pop III star   1, respectively.b -Accretion continues when H16 stops their simulation.
formation that followed the early phase (e.g., Machida et al. 2008b;Clark et al. 2011;Greif et al. 2011Greif et al. , 2012;;Smith et al. 2011;Hirano & Bromm 2017;Susa 2019;Sharda et al. 2019).Subsequently, the protostars grow in mass via gas accretion (Fig. 3, middle) until the radiation feedback quenches it (Fig. 4, right).The resulting stellar system is either a binary of two stars (in the Low-and High-Ṁ cases) or a binary of a single star and a mini-triplet system (in the Intermediate-Ṁ case), as shown in Fig. 3 (right).
To see the evolution, we plot the mass (top), accretion rate (middle), and separation of selected pairs (bottom) of protostars, i.e., sink particles, in Fig. 5.As described above, multiple protostars are formed by the initial disk fragmentation (∆t ≲ 5, 000 yr), grow in mass via gas accretion (∆t ∼ a few × 10, 000 yr), and finally cease to grow due to radiative feedback, with the accretion rates dropping to such a low level that the total masses are almost fixed at the end of the simulations (∆t = 80, 000 or 100, 000 yr).Many protostars disappear in mergers induced by unstable gravitational interaction in a-fewbody systems, and only two to four stars survive at the end of the simulations.
We summarize the properties of the final stellar systems in Table 2.The total masses are M tot = 100 − 500 M ⊙ , with increasing masses toward higher cloudscale accretion rates (e.g., Hirano et al. 2014).The number of surviving protostars is N sink = 2 or 4. The masses of the most and the second most massive stars in the system, M 1 and M 2 , respectively, are all quite large (> 30 M ⊙ ) with the most massive one reaching even 370 M ⊙ .The mass ratio of q 12 = M 2 /M 1 is around 0.3-0.8, which means moderate mass hierarchy.The semi-major axis of their binary orbits is a 12 = 2 × 10 3 − 2 × 10 4 au at the end of the simulations.In all cases, the final outcomes are widely orbiting (> 2, 000 au) massive (> 30 M ⊙ ) multiple stellar systems.
Although only with very small sample size, it is worth comparing the distribution of our binary properties with the binary statistics in the literature, which have been used to estimate the rate of binary black hole mergers.Our sample of massive binaries, characterized by large primary masses (M 1 ≳ 60 M ⊙ ), wide orbits (a 12 ∼ 1, 000 − 10, 000 au), and moderate mass ratios (q 12 ∼ 0.3 − 0.8), roughly agrees with the binary statistics of Liu et al. (2021, see their Fig.9), who considered the effect of orbital expansion due to the conservation of angular momentum.The statistics of Liu et al. (2021) predict significantly fewer close binaries than other studies that assume simple models without the above expansion effect (e.g., Kinugawa et al. 2014;Belczynski et al. 2017).

Some characteristic phenomena
Next, we describe characteristic phenomena observed in our simulation runs.Some of them are common to all the runs, while others occur only in a limited run(s), partly explaining the reason for the observed similarity or diversity of the formed systems.In what follows, we describe those processes one by one in each run.

Low-Ṁ run
In the Low-Ṁ run, we identified the following two characteristic processes, both of which are also observed in the two other cases.The first is initial disk fragmentation, which seeds multiple protostars, and subsequent a-few-body interaction, which induces mergers between them.The second is photo-evaporation of a circumstellar disk by another protostar.
Initial disk fragmentation followed by a-few-body interaction among fragments -At the end of the gravitational collapse, the first protostar appears at the center of the rotating cloud (see Figs. 1 and 2).Around the protostar, an accretion disk soon forms thanks to the accumulation of the angular momentum from the rotating envelope.From top to bottom, the masses, accretion rates (averaged over 300 yr), and distances between selected pairs are plotted.In the top two panels, the same color is used for the same protostar, whose ID is indicated in the top panel, while a combination of the two colors of the member stars is used to indicate the pair in the bottom panel.The black dashed line in the bottom panel indicates the distance below which two sink particles are assumed to merge.In the Intermediate-Ṁ case, we do not follow the individual dynamics of the mini-triplet after t = 5.5 × 10 4 yr (dashed, see Paper I).In the High-Ṁ case, the second sink particle does not appear in the plot as it merges with the first one at ∆t = 100 yr, before the beginning of the plot.

Low-𝑴 15
Since the disk is relatively massive compared with the central protostar in the early phase (see Kimura et al. 2021), the disk promptly fragments and produces two new protostars at ∆t ≈ 2, 000 and 4, 000 yr (Fig. 5).A clump in a spiral arm (at lower-left in the upper left panel of Fig. 3 at ∆t = 3, 000 yr) will later develop into the third protostar.
After formation by the disk fragmentation, those three protostars interact with each other and two of them eventually merge, as seen in the top panel of Fig. 6, which depicts their trajectories of the before and after the merger.The trajectories of the three protostars change violently since a three-body system is generally unstable.The interactions tend to be strong in the earliest phase because the protostars have similar masses and short mutual distances, comparable to the disk radius.After the merger of two protostars at ∆t ≈ 20, 000 yr, a wide-orbit eccentric binary system is left, as shown in Fig. 6 (bottom).
Changing their orbits, the protostars grow in mass by gas accretion, as seen in Fig. 5 (upper left).The accretion rate onto each protostar remains roughly constant until the merger event, where the gas in the disks is stripped away by the gravitational interactions and the accretion rate drops suddenly (Fig. 5, middle left).Photo-evaporation by the central protostars causes mass loss from the accretion disk surface (see Fig. 4, top right) and reduces the accretion rate even more.Around this time, the protostars left are scattered into a low-density peripheral region and gas supply to the disks is also quenched (see Fig. 6, bottom).Although accretion is enhanced temporarily due to the perturbation by the other star on the circum-stellar disk around the pericentric encounter at ∆t ≈ 50, 000 yr (Park et al. 2023), associated mass growth is modest.
Disk photo-evaporation by a nearby protostar -After the pericentric encounter, the accretion rates soon drop again.One of the stars completes the photo-evaporation of its own accretion disk and then begins to photoevaporate the disk of another nearby star.The time sequence of this "external" photo-evaporation is shown in Fig. 7. Soon after the photo-evaporation of its own disk by the left star, the ionized region readily expands into a surrounding lower-density region and finally engulfs the disk on the right side.
We must remark that the external photo-evaporation in this case does not play a significant role in setting the final stellar mass.This is because, by this time, the accretion rate has already fallen to a value too low to change the mass.
Note that such external photo-evaporation is observed not only in the Low-Ṁ case, but also in all other cases.In our runs, its role in controlling the evolution of + 1,000 yr + 3,000 yr + 10,000 yr Δt = 64,000 yr + 0 yr Figure 7.The external photo-evaporation event shown as a time sequence in the Low-Ṁ case.From top to bottom, we show the gas (blue to red) and the ionization front (yellow) at the epochs of 0, 1, 000, 3, 000, and 10, 000 yr after ∆t = 64, 000 yr, when the left protostar begins to completely photo-evaporate its own disk from inside.protostellar systems is limited: although this may explain peculiar late-time orbital evolution seen in the Intermediate-Ṁ case (see Section 3.2.2), it does not affect the final stellar mass significantly in any runs.External photo-evaporation, however, may play an important role in setting the mass of low-mass stars that are unable to clear the disk material away by their own radiation alone.The reason that such stars are not observed in our current runs is partly due to resolution limitations.In any case, given the prevalence of external photo-evaporation, its role in Pop III star formation is worth further investigation.For example, in Galactic star-forming regions, photo-evaporating disks by external irradiation, also called "proplyds," are of interest because they not only leave interesting observational tracers but also have effects on the disk evolution, including planet formation (e.g., Yorke & Welz 1996;Haworth et al. 2017).In the Intermediate-Ṁ case, we also observe the initial disk fragmentation followed by a-few-body interaction, as well as the external disk photo-evaporation.The outcome of the former event is, however, different here: a circular, rather than eccentric, binary system is left, unlike in the Low-Ṁ run, possibly due to the chaotic nature of a-few-body systems.The binary then accretes gas from the circum-binary disk (Figs. 3, center) around the member stars' circum-stellar disks.Eventually, one of them fragments, leading to the formation of a minitriplet system (Figs.3, middle right).Here, we describe those two processes specific to the Intermediate-Ṁ case.

Intermediate-Ṁ case
Binary evolution by accretion from the circum-binary disk -After the initial interaction phase (∆t ≳ 6, 000 yr), the binary's mass and separation significantly change, as shown in Fig. 5.This is due to accretion from the circum-binary disk to the protostars through the circumstellar disks (Fig. 3, center).To see this clearly, we investigate the structure of the accretion flows on the circum-binary disk scale at ∆t = 10, 000 yr. Fig. 3 (center) illustrates the binary system with a separation of ≈ 1000 au embedded in a circum-binary disk spanning from ≈ 2000 au to the outer radius.This configuration is also evident in the 1D surface density profile presented in Fig. 18 (top) of Appendix D. Gases are transferred from the circum-binary disk to the circum-stellar disk of each star through the interfaces, i.e., on the left (right, respectively) side of the left (right) star, while gaps are created on the upper and lower sides.
The circum-binary disk rotates at a rate slightly lower than the Keplerian for the enclosed mass, i.e., the sum of protostars and gas (see Fig. 18 third row in Appendix D for the 1D angular velocity profile).Since its specific angular momentum is higher than the binary's orbital value, the accretion from the disk provides the binary with angular momentum.This explains the observed increase in binary separation with mass (e.g., Muñoz et al. 2019;Moody et al. 2019;Chon & Hosokawa 2019;Tiede et al. 2020;Heath & Nixon 2020).
Driving mechanism of accretion, i.e., of the angular momentum transfer, varies depending on the flow segment.In Fig. 8 (top), we present the face-on view of the Toomre's Q parameter, defined as 1/2 the epicyclic frequency, and Ω the angular frequency around the barycenter of the binary.This parameter is an indicator of gravitational stability: a region of the disk is stable if Q > 1 and unstable otherwise.Fig. 8 (top) shows that the circum-binary disk is marginally unstable with Q ∼ 1 at outer radii ≳ 2000 au (also see bottom panel of Fig. 18 in Appendix D).Within the circumbinary disk, the accretion is driven by the gravitational torque exerted by spiral arms, which are continuously generated and disrupted, as is commonly assumed in 1D models of a circum-stellar accretion disk (e.g., Matsukoba et al. 2021;Kimura et al. 2021).On the other hand, in the region between the circum-binary and the circum-stellar disks (≈ 1500−2000 au; see Fig. 3, center) Q > 1, i.e., the gas is stable against self-gravity, suggesting that accretion is primarily driven by the torque exerted by the binary's gravity.On the smaller scale within each circum-stellar disk, the accretion is again governed by disk's self-gravity (Q ∼ 1), as presented in Appendix E. The circum-stellar disks have a similar structure to the conventional (quasi-)axisymmetric one in the literature (e.g., Hosokawa et al. 2011;H16).
The vertical structure of the disks is shown in Fig. 8 (bottom).The circum-binary disk puffs up to > 1, 000 au in the vertical direction (also see fourth panel of Fig. 18 in Appendix D) due to its weak gravity, in contrast to the thin circum-stellar disks of a few hun-dred au (see Appendix E).At this moment, bipolar HII regions are still confined to the vicinity of the protostars and radiative feedback is insignificant (Fig. 4 center; see also Appendix E).
From the analysis above, a unified picture can be drawn regarding accretion flows, which cause the binary's mass growth, as well as orbit expansion.Accretion proceeds first from the circum-binary disk to the circum-stellar disks and then to the protostars, and the mechanism of angular momentum transfer depends on the region of the flow.
Finally, we comment on a possible effect of radiative feedback on the binary orbital evolution.The increase of the binary separation can be ascribed to the accretion of high angular momentum gas during the period ∆t ≲ 30, 000 yr, where the accretion onto stars is in fact vigorous.The separation, however, continues to increase thereafter despite no significant mass growth.This requires different explanation.One possibility is a backward/forward asymmetry in density around the star due to radiative feedback.Indeed, black holes moving in the ambient gas are known to be accelerated by such asymmetries (Park & Bogdanović 2017;Toyouchi et al. 2020;Sugimura & Ricotti 2020).Another possibility is reduction in the gravitational pull on the binary stars due to the photo-evaporative loss of gas in the region between the binary.The effect of radiative feedback on orbital evolution merits further investigation.
Late-time disk fragmentation and acceleration of photoevaporation -During binary accretion, one of the circumstellar disks fragments to form a mini-triplet system at ∆t = 20, 000 yr (see Fig. 5).The state of the disk before (∆t = 18, 000 yr) and after (∆t = 30, 000 yr) fragmentation is shown in Fig. 9.The fragmentation of a spiral arm in the circum-stellar disk (Fig. 9, top) provides new protostars (S6 and S7 in Fig. 5, upper middle), making a mini-triplet system (Fig. 9, bottom).The central protostar is far more massive than the companions, as in planetary systems.The mini-triplet inherits the small angular momentum of the circum-stellar disk and thus has a compact size of ≲ 1, 000 au, in contrast to the wide outer binary at a distance of 6, 000 au, which carries the large angular momentum of the circum-binary disk.
Here, one circum-stellar disk fragments and the other does not.Although they are similar before fragmentation (∆t = 18, 000 yr), the former is slightly more unstable with a few ×10 % more mass than the other and thus causes fragmentation.
Following the formation of the mini-triplet by disk fragmentation, the gas is quickly lost in the forming stars and their subsequent accretion, which accelerates the photo-evaporation of the disk (see Fig. 9, bottom).The photo-evaporation of the gas around the mini-triplet is completed much earlier than the counterpart in the wide outer binary, at ∆t = 50, 000 yr.After that, the hierarchical triplet system remains stable and shows no obvious evolution in our simulation, while the ionized region around the mini-triplet expands, reaches the other circum-stellar disk, and photo-evaporates it from the outside, as described in Section 3.2.1.

High-Ṁ case
In terms of event sequence, the High-Ṁ case resembles more with the Low-Ṁ case than the Intermediate-Ṁ case.We observe initial disk fragmentation and external photoevaporation (Section 3.2.1)but do not either accretion from a circum-binary disk or late-time disk fragmentation (Section 3.2.2).As a process unique to the High-Ṁ case, below we explain turbulent fragmentation preceding initial disk formation.
Turbulent fragmentation preceding disk formation -We present the gas distribution around protostars in Fig. 10 at the following three stages: initial turbulent fragmentation (∆t = 37 yr), the subsequent disk formation (∆t = 2, 000 yr), and just before the late-time merger of the inner binary (∆t = 18, 000 yr).
In the top panel, we can identify a fragment labeled with "(S3)" (not yet converted to a sink particle) on the upper side along a vertical filamentary structure, in addition to two protostars (S1 and S2) on the other side.Note that S2 soon merges with S1 at ∆t = 100 yr and does not appear in Fig. 5, which starts at ∆t = 1, 000 yr.The turbulent fragmentation occurs only in this run, presumably because the cloud has strong turbulence and weak rotation in the initial collapse phase (Figs. 1 and  2).
Eventually, the accumulation of angular momentum leads to the formation of a disk around the central protostar (S1), as shown in Fig. 10 (middle).The disk subsequently becomes gravitationally unstable and undergoes fragmentation, seeding several protostars (see Fig. 5, right).During the disk formation process, the gas originally belonging to the filament also accretes onto the disk.Simultaneously, the fragment, which transforms into the third sink particle (S3), is gravitationally captured by the central protostar.It then begins orbiting within the disk plane through the interactions with the disk gas, forming a close binary with ∼ 300 au separation.The short separation can be attributed to the relatively low initial angular momentum of the protostar formed through filament fragmentation, compared with those formed through disk fragmentation.
After the disk fragmentation, many protostars merge through gravitational interaction, leaving a three-body system wherein an outer protostar (S7) orbits an inner binary (S1-S3), as shown in Fig. 10 (bottom).Soon after the epoch shown in the figure, however, the longlasting inner binary, originated from the capture of a filament fragment around ∆t ≈ 2, 000 yr, eventually merges due to the perturbation by S7 at ∆t = 20, 000 yr.This demonstrates that in a hierarchical three-body system, the outer-orbit star can extract the angular momentum from the inner binary via gravitational torque and decrease its separation, rather than increasing the binary separation, as observed in the case of circum-binary disk accretion (Section 3.2.2).We here assumed that a binary has merged when the separation reaches 128 au in our criterion (see Section 2.1.1).In reality, however, our merger events could represent close-binary formation below the resolution limit.This suggests potential importance of turbulent fragmentation in close-binary formation.
After the merger of the inner binary, an eccentric binary of the merger product and the outer star is left (Fig. 3, lower right), similar to the Low-Ṁ case.We observe periodic modulation of the accretion rates onto the two stars with pericenter encounters of the binary in Fig. 5 (Park et al. 2023).While the accretion rates modulate at the same frequency, their amplitudes are different due to the difference in the stellar gravity.The lower-mass star has the accretion rate peaking as high as ∼ 10 −2 M ⊙ /yr, which results in temporal stellar swelling and significant suppression of UV radiation (H16; Park et al. 2023; see Appendix B).Later, with gradual decline of the peak, as well as average, accretion rate, the UV suppression ceases.

RELATION BETWEEN THE PROPERTIES OF THE FINAL STELLAR SYSTEM AND NATAL CLOUD
In all the three runs, we have observed the formation of massive and wide multiple-star systems, with quantitative variations among them.In Section 4.1, we see the relation between properties of the final stellar systems, i.e., the total mass and binary separation, and those of initial clouds.We then discuss why Pop III stars predominantly form as massive and wide multiple-stellar systems based on those relations in Section 4.2.

Dependence of the total mass and binary separation on the cloud properties
The relation between the final stellar mass and the initial cloud-scale accretion rate has been proposed based on 2D simulations of single Pop III star formation (Hirano et al. 2014).Our 3D simulations, however, have shown that multiple, rather than single, stellar systems Figure 11.The relation between the final total stellar mass, Mtot and the initial cloud-scale accretion rate, Ṁcloud .We plot the data from our 3D AMR simulations (blue), along with previous 3D spherical-grid simulations (H16, orange), and previous 2D simulations (Hirano et al. 2014, purple).The 2D data points are denoted by small symbols, except for those re-examined by 3D spherical-grid simulations.The dashed line shows the relation proposed by Hirano et al. (2015) based on the 2D results (Eq.13).
10 57 10 58 J enc [g cm 2 /s] are formed in reality.Below, we see whether similar correlation still exists in our case.Fig. 11 shows the relation between the final total stellar mass, M tot , and the initial accretion rate at the cloud scale, Ṁcloud (see Table 1).We plot our results together with those from previous 2D axisymmetric simulations (Hirano et al. 2014) and 3D spherical-grid simulations (H16).In both cases, only a single star forms in each cloud due to the numerical limitation, while accretion modulation induced by disk fragmentation is observed in the latter.We find that the total mass tends to be higher for a higher accretion rate also in our simulations, in agreement well with the fitting function proposed by Hirano et al. 2015 (their Eq.7): Formation of multiple protostars have two effects on the total mass.On the one hand, mass sharing among multiple stars leads to smaller mass of each star and thus weaker radiative feedback (see Appendix B).On the other hand, the displacement of protostars towards less dense off-centered regions leads to the suppression of accretion.In addition, the chaotic behavior of a-fewbody systems may also introduce extra scatters in each case (e.g., Susa 2019;Wollenberg et al. 2020).
The stellar orbits also correlate with the properties of the natal clouds.Fig. 12 shows the orbital angular momentum of the binaries, L orb , against the corresponding initial enclosed gas angular momentum, J enc .To determine these values, we first evaluate the binary's orbital angular momentum and total mass at ∆t = 30, 000 yr.We then measure the initial enclosed angular momentum within the radius containing the total mass (3,000, 4,000, and 9,000 au for the Low-Ṁ , Intermediate-Ṁ , and High-Ṁ cases, respectively) for the cloud profile just before the formation of the first protostar (see Fig. 2).Note that these radii are somewhat smaller than the cloud size, in which the cloud spin parameter λ cloud is measured in Table 1.The mini-triplet in the Intermediate-Ṁ case is treated as a single object put at their barycenter.This does cause little error in evaluating the angular momentum of the system as the angular momenta of the mini-triplet's internal orbits are negligible compared with that of the wider binary consisting of the virtual body and the other protostar.
We chose the epoch of ∆t = 30, 000 yr in order to eliminate the influence of radiative feedback on the binary's orbit and focus on quantities primarily determined by the conservation laws for mass and angular momentum.By this epoch, the accretion process is nearly completed, as shown in Fig. 5.As we see in the Intermediate-Ṁ case, radiative feedback possibly makes the binary's separation change even after the end of accretion.Although some more study on radiative feedback effect on the stellar orbits is needed, we speculate that its impact on the (already large) binary separation is rather limited.Fig. 12 shows that the late-time orbital angular momentum of the binary/multiple system, roughly agrees with the initial angular momentum of the cloud.This means that large angular momentum/separation of the binaries can be attributed to large angular momentum of the clouds.Slightly lower orbital angular momentum in the High-Ṁ case than is expected is probably due to strong turbulence in the initial cloud (Fig. 2), which transfers the angular momentum outward.

Reason for massive and wide multiple-star systems
When the first protostar is still in the infancy, a disk forms around it.Being relatively massive compared with the central star, the disk easily fragments and forms multiple protostars (see Kimura et al. 2021).Such smallbody systems are generally unstable by the gravitational interaction and most of them merge or are scattered away.Eventually, only a few protostars can stay within dense regions near the center and continue to grow in mass.In the following, we limit the case of binaries for simplicity although we expect similar conclusion in a-few body cases.
While the binary system gains mass by accretion from the circum-binary disk, its mass ratio approaches unity in accordance with previous dedicated studies (Bate & Bonnell 1997;Satsuka et al. 2017;Matsumoto et al. 2019;Chon & Hosokawa 2019).Note that accretion does not cause the merger of pre-existing binaries: instead it increases the binary separation (see Fig. 8 bottom), rather than decreasing it (Tiede et al. 2020;Heath & Nixon 2020).Consequently, Pop III stars tend to form as a system in which a few stars dominate the total mass, and presumably the total orbital angular momentum.
This system then becomes massive and wide.The total mass reaches as high as 80 − 500 M ⊙ , thanks to the high cloud-scale accretion rate (Fig. 11), which stems from the high temperature in the primordial gas due to the inefficient cooling (e.g., Stahler et al. 1986;Omukai & Nishi 1998).Simultaneously, the (outer) binary has a wide orbit with its separation reaching at least 2, 000 au due to the large initial angular momentum of the natal cloud (Fig. 12) without effective angular momentum extraction by such processes as magnetic braking or magnetically-driven outflows (e.g., Machida et al. 2008a;Sadanari et al. 2021Sadanari et al. , 2023)).In this way, Pop III stars predominantly form as multiple systems consisting of massive stars with wide orbits.Note, however, that this does not exclude the formation of close binaries, some of which could be progenitors of binary BH mergers observed by gravitational waves.Rather, our results have some implications for close binary formation.As in the Intermediate-Ṁ run, Pop III stellar systems can be hierarchical, in which outer wide orbit stars and mini-multiplet systems coexist.Whereas the separation of the outer binary tends to increase by accretion from the circum-binary disk, stars in the outer orbit can extract the angular momentum from the inner binary and shrink its separation of the inner one, as observed in the High-Ṁ case.Although the inner binary was regarded as merged in our simulation, the actual end product could be a close binary system below our resolution limit.Higher resolution simulation is needed in the future to answer whether such close binaries as gravitational event progenitors are formed or not among Pop III stars (e.g.Kirihara et al. 2023).
Note also that low-mass stars may still form although we have not found them in our simulations, possibly due to our limited spatial resolution and simplistic sink merger criteria (see also Sec. 5.2).Previous numerical studies on cosmological Pop III star formation have suggested that numerous stars remains low-mass (< 1 M ⊙ ) after ejection from the central region of the star-forming cloud (Greif et al. 2012;Latif et al. 2022).Whether or not such stars indeed coexist in the system, however, is unlikely to significantly affect the evolution of massive protostars.Nevertheless, their formation and survival are of substantial interest from observational perspectives as current observations have already ruled out the cases where an excessive number of low-mass stars survive up to the present (Hartwig et al. 2015;Ishiyama et al. 2016).
In summary, we argue that Pop III star systems are likely comprised of widely orbiting multiple massive stars.At the same time, these systems may include embedded close binaries, as well as numerous low-mass stars ejected from the center.

Effect of radiative feedback
The radiative feedback from protostars suppresses the accretion growth of protostars, as observed in our simulations.To see how this works more specifically, we here perform a series of additional runs for the Intermediate-Ṁ case with different numerical setups.Previous authors suggest that the accretion is suppressed either by FUV (e.g., Susa 2013;Susa et al. 2014) or EUV radiation (e.g., McKee & Tan 2008;Hosokawa et al. 2011).To discriminate the effect of each UV component, we perform runs with only FUV (no EUV) and without radiation (neither EUV nor FUV), in addition to the fiducial run both with EUV and FUV radiation.Table 3 summarizes the setups for the additional runs, along with another series of runs presented in Section 5.2.Fig. 13 gives the comparison of the protostar evolution in the fiducial (EUV and FUV), FUV-only, and noradiation runs.The total mass evolution (top panel) shows variation among the three runs.In the noradiation run, the mass increases linearly at a nearly constant accretion rate without any feedback.On the other hand, in the runs with FUV radiation, the mass growth slows down due to photodissociation feedback around ∆t ∼ 10, 000 yr.In the FUV-only run, the mass continues to grow at a reduced rate, while the EUV photoionization feedback nearly halts the accretion at ∆t ∼ 40, 000 yr in the fiducial run with EUV.The bottom two panels in Fig. 13 show the number of surviving sink particles, N sink , and the total number of sinks formed, N form . Due to a-few-body interaction, which reduces the number of protostars through mergers, N sink is similar in the range of 3 -5 in the three runs.On  the other hand, a larger difference is observed in N form .
While N form continues to increase with time in the no radiation run, it reaches a plateau around ∆t ∼ 20, 000 yr in both the fiducial and FUV-only runs.This plateau is likely caused by the reduction in the accretion rate onto the disks, which falls below the rate required to sustain the formation of new sink particles through disk fragmentation.
To see how radiative feedback affects the thermal state of gases, we present density-temperature phase diagrams in Fig. 14.The color in each bin represents either the mass (top) or the H 2 fraction (bottom).All the plots are for the epoch when the total mass is 130 M ⊙ .Around this time, FUV radiation, if included, has already suppressed the accretion and EUV radiation, again if in-cluded, begins to terminate the accretion (see Fig. 13, top).The corresponding epochs are ∆t = 30, 000 yr in the fiducial and FUV-only runs and ∆t = 15, 000 yr in the no-radiation run.
The presence of FUV feedback affects the envelope gas in the density range of 10 5 − 10 9 cm −3 .A substantial portion of the gas in the envelope is heated to T ≳ 1, 000 K in the runs with FUV, while most of the envelope gas remains at T ≲ 1, 000 K in the no radiation run (Fig. 14, top).The higher temperatures under FUV feedback can be attributed to the reduced H 2 fraction and thus cooling due to FUV photodissociation (Fig. 14, bottom).The suppression of mass growth at ∆t ≳ 10, 000 yr in the runs with FUV radiation (Fig. 13, top) is presumably caused by high temperature and thus pressure in the envelope.
EUV feedback leads to the formation of bipolar ionized bubbles around protostars, as depicted in Fig. 4.These ionized bubbles are observed as a hot (20, 000 K < T < 40, 000 K) and low-density (10 5 cm −3 < n H < 10 7 cm −3 ) component in Fig. 14.Of course, this component is observed only in the fiducial run with EUV radiation and not in the runs without EUV radiation.Consistent with previous studies (e.g., McKee & Tan 2008;Hosokawa et al. 2011), the EUV photo-evaporation of accretion disks appears to be responsible for the observed decline in accretion at ∆t ≳ 40, 000 yr in the fiducial run (Fig. 13).
The analysis above highlights the distinct roles played by FUV and EUV radiation in halting the accretion.The FUV feedback starts suppressing the accretion since an early phase by quenching the main coolant H 2 , but it alone is not able to completely halt the accretion.Later on, the EUV feedback becomes active and eventually terminates the accretion by the disk photo-evaporation, finally fixing the stellar mass.Given different roles of EUV/FUV radiation in Pop III star formation, we caution that ignoring either component may lead to unrealistic results.

Resolution effects
To accelerate the computations and follow the accretion evolution to the end, we introduced sink particles.Still sink particles at finite resolution give rise to some uncertainties.Here, we assess the resolution effect.
For this purpose, we perform additional runs for the Intermediate-Ṁ case with varying resolutions, as summarized in Table .3: two runs without radiation have two and four times higher resolution and one run with fiducial feedback model, including both EUV and FUV, has two times higher resolution.The ratio of the sink radius to the minimum cell size is fixed at r sink /∆x min = 16 for all runs, resulting in a sink radius of r sink = 32 (16) au in the two (four) times higher resolution.Furthermore, we ensure that the Jeans length (proportional to (n sink ) −2 for a constant temperature) matches 2 r sink by increasing n sink by a factor of four (16) when adopting two (four) times higher resolution.
Fig. 15 shows the resolution dependence of sink evolution in the runs without radiation, as in Fig. 13.The total mass, M tot , (top panel) is consistent among runs with varying resolutions.The number of surviving sinks (middle panel) is N sink ∼ 4 due to merger in all runs while in the bottom panel we observe somewhat earlier rise of the total number of sinks ever formed, N form , for higher resolution runs.This can be understood from the fact that, in lower resolution runs, it takes longer for the disks to acquire enough mass for fragmentation into sinks.From the consideration above, we conclude that our results are insensitive to the resolution within the examined ranges.
Similarly, we compare the different resolution runs with radiation in Fig. 16.Due to limited computational resources, it is not feasible to continue the high resolution run beyond ∆t = 3 × 10 4 yr, despite ongoing accretion (with significant impact of radiative feedback) at this stage.Apart from the earlier rise in N form in the higher resolution run, the evolution of M tot , N sink , and N form is in good agreement, as in the cases without radiation.This also confirms that the resolution dependence in simulations with radiative feedback is again insignificant within the examined range.
It is nonetheless plausible that some small-scale phenomena near protostars are not fully captured in our simulations.For instance, Prole et al. (2022a) reported the formation of more protostars in the initial disk fragmentation phase at higher resolution than ours.We, however, speculate that properties of massive stars would remain largely unaffected given only a few proto- stars can grow massive within the dense central region, as seen in Section 4.2.
Note that sink merger criterion is linked to the resolution issue.We employ a simplified approach where two sinks are assumed to merge if their sink spheres overlap.This criterion grossly simplifies reality, and the number of surviving sinks should be viewed as a conservative lower bound.More stars, and even close binaries below our resolution, could survive all the way through the end of the accretion phase.
To eliminate all those ambiguities associated with introduction of sink particles, it is desirable to directly resolve the protostars themselves.Although still in its infancy, such an attempt is currently underway (see Kimura et al. 2023).Nevertheless, we expect our conclusion of the formation of massive and wide binary/smallmultiple stellar systems remain unchanged, as discussed in Section 4.2.

SUMMARY AND CONCLUSIONS
We have studied the formation of Pop III stars by performing radiation hydrodynamics simulations for three different primordial clouds extracted from cosmological hydrodynamics simulations.To accurately follow the formation of multiple protostars through gas fragmen-tation and their radiative feedback on the surrounding gas at a reasonable computational cost, we have developed a new code SFUMATO-RT, which employs the adaptive-mesh-refinement and the adaptive-ray-tracing techniques.Starting from the cloud collapse stage, we follow the growth of protostars for ∼ 10 5 yr until their radiative feedback quenches the accretion and nearly fixes the stellar properties.Below, we summarize our findings.
(i) Pop III stars predominantly form as massive and wide binary/small-multiple-star systems in all three cases, with masses of 30 − 370 M ⊙ and separations of 2 × 10 3 − 2 × 10 4 au (Section 3).
(ii) The properties of the formed stellar systems correlate with those of their natal clouds (Section 4.1): the total mass of the system increases with the cloud-scale accretion rate, and the angular momentum of the binary orbit matches that of the cloud.
(iii) The total mass of the formed stellar system is consistent with that in previous simulations of singlestar formation (Hirano et al. 2014(Hirano et al. , 2015;;H16).The individual masses decrease due to mass sharing among the multiple stars.
(iv) The reason for the formation of massive and wide binary/small-multiple star systems is discussed in Section 4.2.The large total mass is due to the high accretion rate of the primordial gas, while the large angular momentum (i.e., large separation) is due to the absence of effective angular momentum transfer mechanisms via magnetic fields.
(v) The following processes are identified as characteristic of Pop III star formation: • A disk forms around the first protostar in the simulation.It then fragments and seeds multiple protostars.Subsequent unstable a-few-body interactions lead to the merger or scattering of most protostars from the central region, limiting the number of protostars that grow into massive stars.Consequently, Pop III stars form as a multiplestellar system, with a few stars dominating the system's total mass.
• Protostars that complete the photo-evaporation of their own disks earlier than the others subsequently externally photo-evaporate the disk around another protostar in their vicinity.While the impact of this process on the final stellar system requires further investigation, it may explain the observed late-time orbital evolution in one of the simulations.Moreover, external photoevaporation may influence the evolution of Pop III disk-star system, as explored in the context of present-day star formation.
• When binary protostars in a circular orbit form, they accrete gas from the circum-binary disk through their respective circum-stellar disks, resulting in an increase in their mass and separation due to the accretion of high angular momentum gas.The fragmentation of a circum-stellar disk leads to the formation of a mini-multiple-star system and accelerates the disk photo-evaporation.The mini-multiple stars have shorter separations than the original wide-orbit binary stars, suggesting that the late-time fragmentation of circumstellar disks may be one of the mechanisms for the formation of short-orbit systems.
• With strong initial turbulence, turbulent fragmentation occurs before disk formation.A protostar seeded by turbulent fragmentation is captured by and eventually merges with the central protostar due to its small angular momentum.The capture of stars originating from turbulent fragmentation may be another mechanism for the formation of short-orbit systems.
Our findings have two major implications.First, the reduction in the individual masses of Pop III stars should change their role in driving the subsequent evolution of the Universe through stellar radiation, supernovae, and seeding BHs that might rapidly grow into supermassive BHs (see, e.g., Sugimura et al. 2017Sugimura et al. , 2018;;Sugimura & Ricotti 2020, and references therein).Incorporating our findings into cosmological simulations of galaxy formation is crucial for unveiling the early evolution of the Universe (see, e.g., Garcia et al. 2023).
Second, the binary Pop III stars in our simulations are massive enough (> 30 M ⊙ ) for the merger of the remnant BHs to be detectable with current gravitational wave detectors if they merge (e.g., Kinugawa et al. 2014).Their separations are too large (∼ 1, 000 − 10, 000 au) for the remnant BHs to merge solely through the binary interaction within the age of the present Universe.However, if these binaries migrate into nuclear stellar clusters during the cosmic structure formation process, dynamical hardening can substantially accelerate the merger, as proposed by Liu & Bromm (2021, 2020).Furthermore, our simulations suggest the possible presence of mini-binaries embedded in a wider system.Such minibinaries may have separations short enough (≲ 1 au) for the remnant BHs to merge via binary interaction (e.g., Kinugawa et al. 2014;Tanikawa et al. 2021).In any case, further high-resolution simulations explicitly resolving such close binaries are necessary to clarify the origin of gravitational wave events.origin of gravitational wave events.
While we have consistently treated the radiation feedback from multiple protostars, we have not considered potentially important effects, such as a magnetic field and a background radiation field.While the strength of the primordial magnetic field is still unknown, it is proposed that the magnetic field is generated and amplified by small-scale dynamo during minihalo formation (Xu et al. 2008;Schleicher et al. 2010;Sur et al. 2010;Turk et al. 2012;Schober et al. 2012;McKee et al. 2020;Stacy et al. 2022).If this is the case, we need to consider magnetic field effects in Pop III formation simulations (see Machida et al. 2006;Sharda et al. 2020Sharda et al. , 2021;;Sadanari et al. 2021Sadanari et al. , 2023;;Prole et al. 2022b;Stacy et al. 2022;Saad et al. 2022;Hirano & Machida 2022).Some Pop III stars are thought to form under the influence of radiation fields from earlier stars.Various types of radiation, including X-rays, EUV, and FUV, can affect the Pop III star formation (e.g., Hummel et al. 2015;Park et al. 2021aPark et al. ,b, 2023)).
The recent launch of the James Webb Space Telescope (JWST) has accelerated the observational quest into the early Universe.To prepare for the wealth of upcoming observational discoveries, it has become more urgent than ever to unveil the early cosmic history through simulations, including Pop III star formation, which represents the very first step toward the formation of various objects in the Universe.
While our simulations have advanced our understanding of Pop III star formation, there is a need to increase the sample size to study the diversity in Pop III star formation.Furthermore, moving to higher resolutions is crucial for a more realistic depiction of Pop III star formation.In our future work, we will extend our research in the directions of higher resolutions and larger sample sizes.

A.2. Thermal processes
The heating and cooling processes considered in this paper are summarized in Table 5.We add to the thermal model of H16 the H free-free emission and the H free-bound emission cooling.Following H16, we do not consider the H 2 CIE cooling, which is ineffective in the density range of our simulations (n H < 10 12 cm −3 ), and becomes the dominant coolant only at higher densities (n H ≳ 10 13 cm −3 ; see Omukai 2001).

B. PROTOSTAR MODEL
We evaluate the radiation from Pop III protostars using the pre-calculated table based on one-dimensional stellar evolution calculations (Hosokawa & Omukai 2009;Hosokawa et al. 2010).The table gives the protostellar radius R * and luminosity L * for a given mass M * and accretion rate Ṁ .Using R * and L * , and assuming rays when either its job list is empty, more than 200 rays are stored in one of the lists for sending, or 2000 times ray advancements have been made.The recipient CPU receives the rays and stores them in its job list.
To check the completion of ray tracing, we define a total job size as N ray,tot = 12 × 2 2 lray,max , which is the number of rays if all rays are split to the maximum HEALPix level l ray,max .In this work, we set l ray,max = 15, which is large enough to ensure that rays are always split as long as the ray-splitting condition is satisfied.We terminate a ray either when it reaches a boundary of the computational box or when it is sufficiently attenuated (specifically, when τ EUV > 100 and N H2 > 10 22 cm −2 (f shield < 10 −8 ) in this work).When a ray at a level l ray is terminated, the corresponding job size, N ray = 12 × 2 2 lray , is added to each CPU's job completion counter, N ray,cpu .At some intervals, the CPUs check via MPI communication whether the sum of N ray,cpu over all CPUs is equal to N ray,tot .If this is the case, the ART procedure is completed.
We have performed several tests to check the validity of our ART implementation, including tests of static and expanding HII bubbles around single radiation sources and a test of static overlapping HII bubbles around two radiation sources.All test results are physically reasonable or agree well with analytical estimates, confirming the validity of our implementation.

D. 1D PROFILE OF A CIRCUM-BINARY DISK
To study more quantitatively the accretion flows on the scale of the circum-binary disk illustrated in Figs. 3 (center),4 (center), and 8, we present the 1D profile of angle-averaged disk quantities in Fig. 18.To obtain the disk quantities, we first set the barycenter as the center and the orbital axis of the binary as the vertical axis.Then, we define the disk surface based on the heights at which the density decreases by a factor of 0.01 from the equatorial plane.Finally, we obtain disk quantities through vertical integration or averaging between the disk surfaces, depending on the variable.
The column density shown in the top panel agrees with the 2D density distribution in Fig. 3 (center).It exhibits a peak around 700 − 1, 000 au, attributed to the circum-stellar disks, a valley between 1, 000 − 2, 000 au resulting from the binary-induced gaps, and a gradually decreasing slope outside, corresponding to the circumbinary disk.
The second panel displays the temperature, indicating that the gas within the circum-binary disk is roughly isothermal.Enhancement near the protostars is likely a result of the radiative feedback from them.In the third panel, we plot the rotational velocity, together with the two Keplerian velocities, considering only the mass of the protostars and considering both the mass of the protostars and the enclosed gas mass.Except for r ∼ 1, 500 au, where the binary's perturbation is significant, the rotational velocity lies between the two Keplerian velocities.This suggests that the gas self-gravity and pressure gradient play substantial roles in the radial force balance within the circum-binary disk.
The fourth panel illustrates the disk height, H disk , and the scale height, H scale = c s /Ω.The high pressure contributing to the radial force balance also influences the vertical force balance, causing the disk to inflate to an aspect ratio greater than unity.Note that the two heights satisfy the relation H disk ≈ 3 H scale for the hydrostatic equilibrium: the heights at which the density decreases by a factor of 0.01 is about 3 H scale in the hydrostatic profile ∝ e −z 2 /2 H 2 scale .In the last panel, Toomre's Q parameter has a marginally stable value of Q ≈ 1 across a wide range of the circum-binary disk.This quantitatively confirms our view in Section 3.2.2 that the accretion within the circum-binary disk is driven by the self-regulated gravitational instability.

E. STRUCTURE ON THE SCALE OF A CIRCUM-STELLAR DISK
Here, we examine the binary accretion flow on the scale of circum-stellar disks.We present the 2D gas distributions in Fig. 19 and the corresponding 1D profiles in Fig. 20.
In Fig. 19, we present the 2D snapshots around one of the circum-stellar disks inside the circum-binary disk shown in Figs. 3 (center), 4 (center), 8, and 18.The density (top), temperature (second), H 2 fraction (third), and H + fraction (last) all indicate that the gas is accreted onto the central protostar through the circum-stellar disk while radiative feedback producing bipolar photoionized and photodissociated regions above and below the disk, which agrees with the (quasi-)axisymmetric case studied in the literature (e.g., Hosokawa et al. 2011;H16).The gas supply mechanism to the disk differs from the axisymmetric case, as the gas enters from one side through a bridge structure, as observed in Fig. 3 (center) and explained in Section 3.2.2.However, this difference does not significantly affect the flow structure on the scale of the circum-stellar disk.Therefore, the knowledge gained from single Pop III star formation remains applicable in the context of binary accretion from the circum-binary disk.
Fig. 20 presents the 1D profiles of the circum-stellar disk illustrated in Fig. 19 in the same way as in Fig. 18.
Notes.Column 1: cloud name, Column 2: collapse redshift, Column 3: dark matter mass of the minihalo, Column 4: cloud radius, Column 5: gas mass of the cloud, Column 6: cloud-scale accretion rate, Column 7: cloud spin parameter.a The High-Ṁ , Intermediate-Ṁ , and Low-Ṁ clouds in this paper are identical to the clouds in Cases A, C, and D in H16, respectively.See Hirano et al. (2015) and H16 for further information about the cloud properties.The Intermediate-Ṁ cloud also corresponds to the cloud studied in Paper I. All cloud parameters are evaluated when the central density reaches nH = 10 7 cm −3 .b The radius at which the ratio of the enclosed mass to the local Jeans mass takes its maximum value (see Hirano et al. 2015).c The gas mass within R cloud .d The accretion rate through the spherical surface at R cloud e The spin parameter defined as λ cloud = J cloud / 2 G M 3 cloud R cloud , with J cloud the angular momentum within R cloud (see Hirano et al. 2015).

Figure 1 .
Figure 1.The state of the cloud cores just before the formation of the first protostars.We show the density rendered image for the Low-Ṁ (left), Intermediate-Ṁ (middle), and High-Ṁ (right) clouds.
Notes.Column 1: run name, Column 2: side length of the simulation box, Column 3: minimum cell size, Column 4: sink radius, Column 5: simulation duration since the first protostar formation, Column 6: number of sink particles surviving until the end of the simulation, Column 7: total mass of the stars, Column 8: mass of the most massive star, Column 9: mass of the second most massive star, Column 10: semi-major axis of the orbit of the most and the second most massive stars, Column 11: stellar mass in H16. a -Low, Int, and High runs are for the High-Ṁ , Intermediate-Ṁ , Low-Ṁ clouds in Table

Figure 3 .
Figure 3. Face-on slice density in the Low-Ṁ (top), Intermediate-Ṁ (middle), and High-Ṁ (bottom) cases at ∆t = 3, 000 (left), 10, 000 (middle), and 30, 000 yr (right) since the first protostar formation.The projected positions of the sink particles are shown as open circles, with their masses indicated in M⊙.Note that the size of the sink particles shown is larger than the actual size for better visibility.In the right column, we show the trajectories of the most and second most massive protostars.

Figure 4 .Figure 5 .
Figure 4. Same as Fig. 3 but for the edge-on slice temperature.The bipolar photoionization and photodissociation fronts are demarcated with solid and dashed lines.The gas velocity is indicated by arrows whose length is proportional to its amplitude.

Figure 6 .
Figure 6.Trajectories of protostars before (∆t = 10, 000 yr, top) and after (∆t = 17, 000 yr, bottom) merger in the Low-Ṁ case.We plot trajectories on the background of a whiteshaded density map.Three-body interaction induces the merger, leading to the formation of an eccentric binary.The white cross in the upper panel indicates the point of merger.

Face
Figure 8.The Toomre's Q parameter (top, face-on view) and number density (bottom, edge-on slice) on the scale of the circum-binary disk at ∆t = 10, 000 yr for the Intermediate-Ṁ case.The dashed contours in the upper panel indicate Q = 1, below which the disk is gravitationally unstable.The positions of the sink particles are indicated with open circles.

Figure 9 .
Figure9.The density distribution of the circum-stellar disk before (∆t = 18, 000 yr, top) and after (∆t = 30, 000 yr, bottom) fragmentation, which leads to mini-triplet formation in the Intermediate-Ṁ case.The gas distribution is shown as a face-on slice together with the sink particles marked by filled circles.In the lower panel, the trajectories of the satellites are overplotted, in the frame where the barycenter of the inner two protostars is fixed, with the masses of the protostars indicated in M⊙.

Figure 10 .
Figure10.The time sequence of initial turbulent fragmentation (∆t = 37 yr, top), subsequent disk formation (∆t = 2, 000 yr, bottom), and just before the merger of the inner binary induced by a star in the outer orbit (∆t = 18, 000 yr) in the High-Ṁ case.The column density distribution is shown in the face-on view.We also show the projected positions of the sink particles with their ID's.At the epoch of the top panel, the fragment labeled with "(S3)" is yet to be converted into a sink particle.The secondary protostar, S2, in the first epoch merges with the primary, S1, and disappears before the second epoch.In the bottom panel, we show the trajectories of protostars in the frame where their barycenter is fixed.The inner binary later merges due to angular momentum extraction by the outer protostar.The axis perpendicular to each image plane is fixed to the same direction.

Figure 12 .
Figure12.The relation between the angular momentum of the binary orbits, L orb , and the corresponding initial enclosed angular momentum of gas clouds, Jenc.The binary's orbital angular momentum and total mass are evaluated at ∆t = 30, 000 yr, before stellar radiative feedback drives orbital evolution.The initial enclosed angular momentum Jenc is calculated for the radius Renc, whose enclosed mass is equal to the total mass Mtot, i.e., Jenc ≡ J<R enc with Renc satisfying M<R enc = Mtot.The plotted points correspond to the Low-Ṁ , Intermediate-Ṁ , and High-Ṁ cases from left to right.

Figure 13 .
Figure13.Dependence of protostellar evolution on the prescription of radiative feedback.From top to bottom, we plot the total mass Mtot, the number of surviving sinks N sink , and the number of sinks formed N form .The line types indicate the runs: the fiducial run with both EUV and FUV (blue solid), the FUV-only run (orange dotted-dash), and the run without radiation (green dashed).

Figure 14 .
Figure14.The density-temperature phase diagram in the fiducial run with both EUV and FUV (left), the FUV-only run (middle), and the run without radiation (right).We take the data when Mtot ≈ 130 M⊙ in all the runs (∆t = 30, 000 yr in the first two runs and ∆t = 15, 000 yr in the last run).The color indicates the mass in the top row and the H2 fraction in the bottom row.In the bottom row, hatched is the region where the gas is ionized (y(H + ) > 0.5)

Figure 15 .
Figure15.Same as Fig.13but for the resolution dependence in no-radiation runs.We compare the runs with different sink radii of r sink = 64 (blue), 32 (orange), and 16 au (green).

Figure 16 .
Figure16.Same as Fig.13but for the resolution dependence in full-feedback runs.We compare the runs with different sink radii of r sink = 64 (blue) and 32 au (orange).

Figure 18 .
Figure18.One-dimensional profile of the circum-binary disk presented in Fig.8.From top to bottom, we plot the angle-averaged column density Σ, temperature T , angular velocity v ϕ , height H disk and scale height H scale (see text for definitions), and Toomre's Q parameter Q.In the third panel, we show the two Kepler velocities considering only the mass of the protostars (green) and considering both the mass of the protostars and the enclosed gas mass (orange).Black dashed lines indicate the aspect ratio of unity, H/R = 1, in the fourth panel and the marginally stable state, Q = 1, in the last panel.The two sink particles are located around 800 au from the center.

Figure 19 .
Figure19.The edge-on slice snapshots for one of the circum-stellar disks at ∆t = 10, 000 yr in the Intermediate-Ṁ case.From the top to bottom, we plot the density, temperature, H2 fraction, and H + fraction.In the lower three panels, the bipolar photoionization and photodissociation fronts are demarcated with solid and dashed lines.

Table 1 .
Summary of initial clouds.

Table 3 .
Summary of simulation setup for additional runs.

Table 4 .
Chemical ReactionsOur chemical network is summarized in Table4.It is the same as in H16, except different treatment of H 2 photoionization.The model of H 2 photoionization has an ambiguity since H 2 direct photoionization requires an energy of at least 15.2 eV but we have only one frequency bin for EUV photons (hν > 13.6 eV).In contrast to H16, we forbid H 2 direct photoionization and allow H 2 to be photoionized by a two-step process: H 2 photodissociation H 2 + FUV → 2 H followed by H photoionization H + EUV → H + + e.