The Birth of a Massive First-star Binary

The first generation of stars in the universe, also known as the Population III stars, are born in minihalos of 10⁵–10⁶ M_⊙ around the redshift z ∼ 20–30 in the framework of the standard Lambda cold dark matter (ΛCDM) cosmology (e.g., Abel et al. 2002; Bromm et al. 2002; see also Glover 2013; Greif 2015 for review). Small embryonic protostars formed by the gravitational collapse of the natal clouds (e.g., Omukai & Nishi 1998; Yoshida et al. 2008) grow in mass by accretion of the surrounding gas (e.g., Omukai & Palla 2003; Tan & McKee 2004), finally reaching the mass of 10–1000 M_⊙ (e.g., Susa et al. 2014; Hirano et al. 2015; Hosokawa et al. 2016, hereafter H16) when the radiation feedback quenches the gas supply (e.g., Omukai & Inutsuka 2002; McKee & Tan 2008; Hosokawa et al. 2011). If born as binaries, such massive stars can be the progenitors of the binary black holes (BHs) with masses ≳30 M_⊙ recently observed in gravitational waves (e.g., Kinugawa et al. 2014; Abbott et al. 2016; Hartwig et al. 2016, but also see Belczynski et al. 2017), as well as the X-ray binaries that affect the thermal history of the intergalactic medium (Mirabel et al. 2011; Jeon et al. 2014), which is a target of future 21 cm line observations (e.g., Dewdney et al. 2009).

Recent numerical simulations have shown that the disks around Population III protostars are gravitationally unstable and fragment into clumps, in which other protostars can form (e.g., Stacy et al. 2010, 2012, 2016; Clark et al. 2011; Greif et al. 2011, 2012; Smith et al. 2011; Hirano & Bromm 2017; Sharda et al. 2019). Although most of these protostars merge with another protostar soon after their formation, some of them may survive for a longer period and may end up with binary or multiple systems (e.g., Chon & Hosokawa 2019; Susa 2019). Unfortunately, however, those simulations covered only early phases of the star formation process, and thus the fate of the protostars remains unknown. To reveal the nature of resulting stellar systems, we need to follow not only the formation of multiple protostars by disk fragmentation but also their long-term evolution under the influence of the protostellar radiation.

In this work, we perform simulation of Population III star formation, self-consistently taking into account the radiation from multiple protostars, by using a newly developed radiation hydrodynamics code, SFUMATO-RT. We follow the entire star formation process until the protostellar radiation feedback terminates the gas accretion to the newborn stellar system. In this Letter, we describe results of a run with a typical initial condition for the primordial star formation. Results for various cases with further analysis, together with a detailed description of our code, will be presented in a forthcoming publication (K. Sugimura et al. 2020, in preparation).

2.1. Radiation Hydrodynamics Simulation

2.2. Initial Condition

Figure 1 shows the time evolution of the sink particles, i.e., protostars, appearing in our simulation. In the figure, masses, accretion rates, and separations are plotted from top to bottom. We give IDs S1–S7 to the protostars in order of the formation time. According to qualitative transitions of the system, we define five evolutionary phases: (a) gravitational collapse, (b) initial fragmentation, (c) binary accretion, (d) late-time fragmentation, and (e) photoevaporation. For each evolutionary phase, we show a snapshot of the face-on surface density and edge-on temperature in Figure 2. In the late-time fragmentation phase, we also present a 3D rendering view in Figure 3. Below, we will describe how star formation proceeds in our simulation, tracing the five evolutionary stages.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

(a) Gravitational collapse. We start our simulation from the initial condition of the cloud 1.5 × 10⁵ yr before the first protostar formation (Figure 2(a)). The gravitational collapse proceeds in a self-similar fashion (Larson 1969; Penston 1969), with the central core having an oblate shape due to the rotation of gas.

(b) Initial fragmentation. As the maximum density increases as a result of the self-similar collapse, the first protostar, S1, forms at the center of the rotating core. Subsequently, the gas with finite angular momentum falls to the vicinity of S1 and a circumstellar disk is formed. The disk is highly gravitationally unstable because of its relatively high mass compared with the newly forming central protostar and vigorously fragments into several protostars (Figure 2(b)). Most of them, however, do not survive as a result of accretion on the central star after inward migration through the disk or merger with others by three-body scattering.

(c) Binary accretion. After the initial fragmentation phase, only two stars, S1 and S3, survive and make up a binary system. They are surrounded by a circumbinary disk and each of them has its own circumstellar disk (Figure 2(c)). The gas accretes from the circumbinary disk to the circumstellar disks, and each star acquires the mass from its own circumstellar disk. The accretion drives the binary evolution in mass and orbit: the total mass increases keeping the mass ratio around unity, and the separation increases due to accretion of higher angular momentum gas from the circumbinary disk. At this moment, double bipolar ionized bubbles with T ≳ 3 × 10⁴ K grow around the binary stars as their EUV emissivities rise (Figure 2(c), right). Note that this phenomenon could not be captured by previous simulations, which could treat only the central radiation source (e.g., Stacy et al. 2016; H16), and become tractable for the first time with our multi-source RT with the ART method. The accretion rate begins to decrease as the flows in the polar directions are quenched by the bubbles and the gas supply continues only in the equatorial directions.

(d) Late-time fragmentation. Due to imbalance between accretion rate from the circumbinary disk to the circumstellar disks and that from circumstellar disks to the stars, the gas accumulates on the circumstellar disks. As a result, the circumstellar disk of S1 fragments into a few objects (Figures 2(d) and 3). After a series of scatter and merger events, a mini-triplet system emerges, with newly formed companion stars (S6 and S7) orbiting around the massive central star (S1). Rapid accretion onto the companions embedded in the disk, together with the enhanced accretion onto the central object by non-axisymmetries, quickly exhausts the disk material. Due to the hierarchical structure of the mini-triplet system both in mass and distance, the orbits of the component stars remain nearly stable. The ionized bubbles continue expanding with more ionizing radiation being emitted as the stars become more massive and with the density in the surrounding envelope decreasing as the collapse of natal cloud proceeds.

(e) Photoevaporation. Thereafter, the disk in the mini-triplet system, which has already depleted significantly by the accretion onto the stars, is lost by the photoevaporation due to the radiation from the central massive star. As a consequence, the ionized bubble around the mini-triplet expands and merges with that around the star S3 (Figure 2(e)). The accretion rate onto S3 is diminished gradually by the radiation not only by S3 itself but also by the mini-triplet. The distance between S3 and the mini-triplet becomes larger partly because the gravitational binding becomes weaker due to the loss of the gas in between them by the photoevaporation. At 5.5 × 10⁴ yr after the first protostar formation, we replace the mini-triplet system with a single sink particle that represents the gravity center, to save the computational cost. Since the properties of the mini-triplet system, namely, the masses and separations among member stars, hardly changed for the last 10⁴ yr, we assume that they remain unchanged in the rest of the simulation.

We stop our simulation at 1.2 × 10⁵ yr after the first protostar formation, when the mass of S3 almost reaches its final value with the accretion rate already reduced to $\dot{M}\sim 3\times {10}^{-5}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$. Extrapolating this decreasing trend of $\dot{M}$ in time, we find the mass of S3 increases by at most ≲M_⊙. We thus expect that further mass growth of S3 is insignificant. This means that the star formation process has been practically completed at the end of the simulation.

The end product at the final time step is a massive binary system consisting of stars with 56 and $67\,{M}_{\odot }$ orbiting each other at a wide separation of 2 × 10⁴ au. The 56 M_⊙ star is also a member of a mini-triplet system with smaller mass companions with 12 and 13 M_⊙ orbiting at 200 and 800 au, respectively. We have observed the formation of massive Population III star binaries starting from the cosmological initial condition.

In this Letter, we have investigated the formation of a Population III stellar system with an initial condition of a typical primordial star-forming cloud, by way of radiation hydrodynamics simulation. When the radiation feedback from the forming protostars quenches the gas supply, a binary system consisting of nearly equal-mass massive stars emerges, with one of the stars making up its own mini-triplet system with less-massive close companion stars. Although we have examined only one case in this Letter, our results suggest that Population III stars are commonly formed as massive binary/multiple systems.

In spite of the same initial condition, we have observed the formation of multiple stars while only a single star was formed in the previous simulation of H16. This difference may come from a difference in resolution. In H16, only a single star survives possibly because of artificial mergers due to the insufficient resolution around companion stars. Besides, the total mass in stars (∼150 M_⊙) in our simulation is smaller than that of the single star (∼300 M_⊙) in H16, which may be related to the fact that the single star always resides at the center of the cloud, where the accretion is not easily quenched by the protostellar feedback because of high density. The total mass being shared with several stars, the mass of an individual object is even smaller in our simulation. This result should be taken into account in discussing the subsequent evolution of the universe as the intergalactic-scale feedback from the Population III stars, e.g., types of SNe (e.g., Woosley et al. 2002), largely depends on their masses.

The number of stars obtained in our simulation should be regarded as the lower limit because of our usage of the sink particle technique and ignorance of what may happen inside the sink particles, such as fragmentation. Moreover, although we assume that overlapping sink particles merge, two or more stars may continue orbiting each other at a distance shorter than the sink radius. Each sink particle, however, is likely to contain at least one massive star because a system consisting only of a large number of low-mass stars is unstable due to mutual gravitational scatterings. Therefore, our conclusion that Population III stars typically form as massive binary/multiple systems is not altered by the possible existence of unresolved objects.

The binary system with 60 and 70 M_⊙ stars found in our simulation is massive enough to be a progenitor of the binary BH mergers observed in gravitational waves (e.g., Kinugawa et al. 2014; Abbott et al. 2016; Hartwig et al. 2016; Belczynski et al. 2017), while the separation is too wide for any interaction to shrink the orbit in the late stage of stellar evolution and the remnants are unlikely to merge within the age of the universe. To examine whether Population III stars can be the progenitors of the observed BH binaries, the distribution of binary separations needs to be known. This is also crucial for predicting the abundance of Population III X-ray binaries, which contribute to heating up the intergalactic medium in the epoch of reionization. For this purpose, we plan to perform similar simulations with a large number of samples in the future. Note that both the masses and separation of our binary system are located near the upper end of the statistical distributions of Population III binaries given by Stacy & Bromm (2013), who took their samples 5 × 10³ yr after the first protostar formation. The masses and separation of our system are already on the large side in the distributions by Stacy & Bromm (2013) at the same timing and become even larger during the later evolution (see Figure 1).

Aside from the massive wide binary, we have also seen the formation of the mini-triplet system. Its separation is much smaller than that of the massive binary because it is formed from the circumstellar disk, which has lower angular momentum compared with the original cloud. On the contrary, the massive binary has a large separation as a result of the accretion of high angular momentum gas from the cloud. In fact, we have seen that the separation of the binary increases as it accretes the gas from the circumbinary disk, as suggested in recent simulations of binary accretion (Duffell et al. 2019; Moody et al. 2019; Muñoz et al. 2019). The late-time disk fragmentation leading to such multiple systems may play some roles in formation of close binaries, which can evolve to the progenitors of BH merger events or X-ray binaries. Additionally, other mechanisms, such as a-few-body scatterings of protostars and angular momentum extraction by magnetic fields, if any, might help shape the close binaries.

We have succeeded in seeing a new evolutionary aspect of the Population III star formation. This is, however, just the beginning of the attempts toward its thorough understandings as more processes, such as sub-sink-scale physics, are still to be clarified. Those issues need to be addressed in future studies.

The authors thank Sunmyon Chon, Michiko Fujii, Masahiro Machida, Massimo Ricotti, Hajime Susa, Hidekazu Tanaka, Ataru Tanikawa, Masauyuki Umemura, Hidenobu Yajima, and Naoki Yoshida for fruitful discussions and comments. K.S. and S.H. appreciate the support by the JSPS Overseas Research Fellowship and JSPS Research Fellowship, respectively. This work is supported in part by MEXT/JSPS KAKENHI grant Nos. 17H02863, 17K05394, 18H05436, 18H05437 (T.M.), 16H05996, 19H09134 (T.H.), 18J01296 (S.H.), and 17H01102, 17H02869, 17H06360 (K.O.). The numerical simulations were performed on the Cray XC50 at CfCA of the National Astronomical Observatory of Japan, the computer cluster Draco at Frontier Research Institute for Interdisciplinary Sciences of Tohoku University, and the Cray XC40 at Yukawa Institute for Theoretical Physics in Kyoto University.