On the Formation and Interaction of Multiple Supermassive Stars in Cosmological Flows

Supermassive primordial stars with masses exceeding ∼105 M ⊙ that form in atomically cooled halos are the leading candidates for the origin of high-redshift quasars at z > 6. Recent numerical simulations, however, find that multiple accretion disks can form within a halo, each of which can potentially host a supermassive star. We investigate the formation and evolution of secondary supermassive stars in atomically cooled halos, including strong variations in their accretion histories driven by gravitational interactions between their disks and those surrounding the primary supermassive stars in each halo. We find that all secondary disks produce long-lived supermassive stars under sustained rapid accretion. We also find, however, that the majority of secondary supermassive stars do undergo at least one protracted quiescent accretion phase, during which time they thermally relax and may become powerful sources of ionizing feedback. In many halos, the two satellite disks collide, suggesting that the two stars can come into close proximity. This may induce additional mass exchange between them, leading to a great diversity of possible outcomes. These range from coevolution as main-sequence stars to main sequence—black hole pairs and black hole—black hole mergers. We discuss the likely outcome for these binary interactions based on the evolutionary state of both supermassive stars at the end of our simulations, as well as prospects for their future detection by current and next-generation facilities.


INTRODUCTION
More than 200 quasars powered by supermassive black holes have now been discovered at z > 6 (Bosman 2021), including nine at z > 7 (Mortlock et al. 2011;Wu et al. 2015;Bañados et al. 2018;Wang et al. 2021).A natural explanation for the seeds of these quasars are direct collapse black holes (DCBHs), which arise due to the formation of supermassive stars (SMSs) in atomically cooled halos at z 20 (Bromm & Loeb 2003;Lodato & Natarajan 2006;Wise et al. 2008;Regan & Haehnelt 2009;Latif et al. 2013).In this picture, a primordial halo grows to ∼10 7 M without forming stars because, for example, it is immersed in strong Lyman-Werner (LW) cerra et al. 2015, 2018;Ardaneh et al. 2018;Luo et al. 2018), or for the longer times required to evolve the flows over many dynamical times but at the cost of resolving fragmentation deep in the accretion disk of the star (e.g., Chon et al. 2018;Regan & Downes 2018;Suazo et al. 2019;Latif et al. 2021).The first simulations to evolve atomically cooled flows for the entire life of a SMS were Latif et al. (2020), which showed that some could form in binaries or small clusters.Most recently, Patrick et al. (2020) followed the collapse of a variety of atomicallycooled halos for 2 − 4 Myr and found that multiple accretion disks could form in a halo, each of which could harbor a SMS.These disks experienced multiple close encounters that gravitationally torqued gas into their centers, triggering brief bouts of massive accretion onto their respective stars.These stars thus interacted indirectly with one another via encounters between their accretion disks.Although the simulations suggested a variety of final outcomes for these objects, they did not evolve the stars themselves.
Recently, Woods et al. (2021) modeled the growth of the stars of the primary disks in Patrick et al. (2020) and found that they collapse at final masses of 1.1 × 10 5 M − 1.9 × 10 5 M .Here, we investigate the co-evolution of the stars in the most massive and stable of the companion disks in Patrick et al. (2020) with the Kepler Lagrangian stellar evolution and hydrodynamics code.Our Kepler models use the time-dependent accretion rates that were tallied at the centers of the disks in the cosmological simulations.We follow the evolution of the SMSs in the companion disks until the disks were subsumed into the primary disks or torn apart by tidal interactions with them.In Section 2 we describe our cosmological and stellar evolution models.The evolution of the SMSs is examined in Section 3. From the interaction histories of the disks and the formation and collapse times of the stars within them we determine if they form SMS -SMS, SMS -DCBH or binary DCBH systems and discuss possible consequences for their detection in Section 4. We discuss prospects for binary interactions in these systems on much smaller scales at later times in Section 5.

enzo
Accretion rates for the companion stars in our study were taken from cosmological simulations with the enzo adaptive mesh refinement code by Patrick et al. (2020).These simulations resolved accretion at the centers of the disks on 0.01 pc scales with six-species primordial gas chemistry (H, He, e − , H + , He + , and He 2+ ) to approximate very high LW backgrounds that ensured isothermal atomic cooling in the halos.The simulations included collisional ionization and excitation cooling by H and He, recombination cooling by H and He, bremsstrahlung cooling, and inverse Compton cooling by the cosmic microwave background (Bryan et al. 2014).The accretion rates were tabulated by computing mass fluxes through a 0.134 pc sphere at the center of each disk at 10 kyr intervals.
The eight halos in Patrick et al. (2020) were chosen to span a range of assembly histories and spin parameters and had masses of ∼ 1 -9 × 10 7 M at collapse at redshifts z = 13.9 -20.4.Resolution at the smallest scales on which fragmentation is known to occur in the disks had to be sacrificed in these simulations to follow their evolution for the times required for the stars to form, evolve and collapse (e.g., Becerra et al. 2015Becerra et al. , 2018)).Our accretion rates therefore exclude these fragments, but they are expected to fall onto the protostar prior to the main sequence (Inayoshi & Haiman 2014).We show accretion rates only for the largest companion disks that form in the halos in Figure 1, not smaller fragments that could undergo partial collapse but are short lived because they are taken up by the main disk or are tidally disrupted by it.

kepler
We follow the evolution of the stars with kepler (Weaver et al. 1978;Woosley et al. 2002).kepler typically partitions each star into a few thousand zones, over which it solves the energy and angular momentum equations: where r is the radius, v is the velocity, P is the pressure, u is the internal energy per unit mass, is the local energy generation rate per unit mass, and L is the rate of energy flow through a shell of radius r.We include a first-order post-Newtonian correction to the gravitational constant, G rel : In Equation 1, the factor Q in the viscous term is where η ν is the dynamic viscosity as defined in Weaver et al. (1978): which includes both the real viscosity, η R , and the artificial viscosity, which can be modified arbitrarily to dampen acoustic oscillations during quiescent phases of the evolution of a star.In this study, we take the standard values of l 1 = 0.1∆r and l 2 = 2∆r, where ∆r is the width of a grid zone.At each time step, Kepler initially attempts to make an arbitrarily large jump ∆t before iterating to find the maximum time step permitted by preset restrictions on the change in radius, temperature, density, luminosity, or velocity between zones, allowing us to model long-lived, quiescent evolutionary phases and to follow the emergence of shocks or the onset of collapse on short timescales.The masses of the stars at collapse triggered by the the post-Newtonian instability (Chandrasekhar 1964) in Kepler are consistent with analytic predictions (Haemmerlé 2021).
To close these equations, nuclear-burning and energy generation are coupled to the hydrodynamics and solved with an adaptive network (Woosley et al. 2004), convection is treated in a time-dependent manner with heat transport following the Ledoux criterion (Weaver et al. 1978), and a Helmholtz-like equation of state is used that incorporates electron-positron pair production, relativistic and non-relativistic degenerate and non-degenerate electrons, and radiation (Timmes & Swesty 2000).We neglect mass loss because wind and pulsational mass losses are expected to be negligible relative to the accretion rates in these stars (Vink et al. 2001;Baraffe et al. 2001;Hosokawa et al. 2013).For simplicity, we neglect rotation although it could affect the stellar structure (Haemmerlé et al. 2018b).
As in previous studies (Hosokawa et al. 2013;Woods et al. 2017;Haemmerlé et al. 2018a;Woods et al. 2021), we initialize all models as 10 M , n = 3 polytropes with central densities ρ c = 10 −3 g cm −3 and temperatures T c = 1.2 × 10 6 K, i.e., at the onset of deuterium burning.We treat accretion onto the star and the associated advection of entropy as described in Woosley et al. (2004) and Woods et al. (2017).The star is evolved until the onset of collapse, the end of the enzo simulation, or if the evolution time corresponds to a moment in the enzo simulation when the companion disk interacted strongly with or merged with the primary disk.

EVOLUTION OF MULTIPLE SUPERMASSIVE STARS IN PRIMORDIAL HALOS
We find that all the secondary disks in the halos in Patrick et al. (2020) form long-lived nuclear-burning stars.Like the stars in the primary disks (Woods et al. 2021), we find a striking variety of internal structures for these secondary SMSs as shown in Figure 2, in marked contrast to stars that evolve under constant accretion rates (e.g., Hosokawa et al. 2013;Woods et al. 2017;Haemmerlé et al. 2018a).The stars all initially exhibit a deep radiative envelope corresponding to the surge in accretion associated with the formation of their natal disks, which proceeds on much shorter times than the star's thermal timescale and leads to the buildup of a steep entropy gradient (e.g., Begelman 2010;Hosokawa et al. 2013;Woods et al. 2017;Haemmerlé et al. 2018a).In several cases, such as Halos 01 and 16, however, almost all of mass of the star eventually lies within its convective core because there is a long ( 100 kyr) quiescent phase ( Kelvin-Helmholtz time-scale) in accretion in which the star can thermally relax and its structure can approach that of an n = 3 polytrope (Chandrasekhar 1964; Woods et al. 2020).In other models, we see the formation of both transient and long-lived convective cells in the otherwise deep, high-entropy radiative envelopes of some rapidly-accreting stars, similar to those in some constant-accretion rate models (e.g., Woods et al. 2017;Haemmerlé et al. 2018a).
In Halo 01, a secondary disk forms at 400 kyr by breaking off from one of the spiral arms of the primary disk.Throughout its evolution, this companion disk remains in an elliptical orbit around the first with typical separations of 0.5 -1 pc.The sharp spike in accretion ∼300 kyr after formation is due to a merger with a clump in the halo.Accretion continues until the second disk merges with the first 650 kyr after formation.At this time, the SMS in the second disk is a somewhat evolved main sequence star with a central hydrogen fraction of ∼0.4.The SMS in the first disk is only slightly more evolved, with a central hydrogen fraction of 0.38.Although the subsequent evolution of these stars is beyond the scope of this work, Woods et al. (2021) assumed that the collision of the disks would trigger a surge in accretion onto the first star but the presence of the second SMS may lead to other outcomes, as discussed in § 4.1.
Halo 02 is an example of a particularly turbulent, chaotic system in which three clumps form and interact with the primary and companion disks (Patrick et al. 2020).Interactions between these clumps and the second disk lead to two large bursts of accretion at 130 kyr and 250 kyr that create a SMS that grows to ∼60 kM by 275 kyr and has a deep radiative envelope.At this point, there is a brief quiescent phase in accretion due to a close encounter between the two disks.The two disks then merge 378 kyr after the formation of the second disk.The first and second SMSs are 62 kM and 143 kM , respectively, and have central hydrogen fractions of 0.58 and 0.34.Notably, the collision of the disks produces a large spike in accretion onto the first SMS that   1. Final masses, ages, evolutionary stages (ES), and central hydrogen fractions for the primary (1) and secondary (2) SMSs in each halo, followed by whether or not they form a close pair by the end of the simulation.Note that the SMS in the primary disk in Halo 02 collapses to a black hole during the merger with the secondary disk.
quickly brings it up to the post-Newtonian instability and causes it to collapse (Woods et al. 2021).
Only one satellite disk forms in Halo 08, which lasts from 685 kyr to about 1 Myr after the formation of the first disk.Accretion onto the star is highly variable because of the eccentric orbit of the satellite disk, and the fluctuations correlate with the smallest and largest distances between the two disks, which are ∼ 0.3 pc and 2.0 pc, respectively.This accretion history produces the distinctive step-like structure in the Kippenhahn diagram of the star in Figure 2.About 1 Myr after the formation of the primary disk, the star in the second disk has an age of 503 kyr, a mass of 99 kM , and a central hydrogen fraction of 0.56.The SMS in the primary disk, however, has by this time collapsed to a black hole, having encountered the post-Newtonian instability at an age of 950 kyr when it reached a mass of 186 kM and a central hydrogen fraction of 0.38.
The first disk in Halo 10 is stable and does not begin to fragment for 1.14 Myr.The second disk forms ∼ 1.2 pc from the center of the first in an initially highly elliptical orbit.Two major episodes of accretion drive the growth of the second SMS that in each case build up a deep, high-entropy envelope.Between these episodes there is a long (nearly 100 kyr) quiescent phase during which the star thermally relaxes without becoming entirely convective.This star reaches a mass of 53 kM and a central hydrogen fraction of 0.55 by the end of the simulation.At this time the secondary disk is still in an elliptical orbit around the first with separations that vary from 0.5 -1 pc.The SMS in the first disk, however, has collapsed because of the post-Newtonian Chandrasekhar instability late on the main sequence at an age of 1.95 Myr, a final mass of 132 kM , and a central hydrogen fraction of 0.06.
The primary disk in Halo 12 is highly turbulent and frequently fragments, but most of the clumps soon migrate to the center of the disk.The longest-lived of the clumps forms at ∼ 1 Myr and produces its own disk that survives for 380 kyr.Rapid accretion in this disk creates a SMS with a deep radiative envelope that grows to a mass 79 kM with a central hydrogen fraction of 0.57 by the time the two disks merge.About 200 kyr before the destruction of the second disk, the SMS in the primary disk collapsed via the Chandrasekhar instability while still on the main sequence at a final mass of 178 kM and central hydrogen fraction of 0.32.
Like Halo 10, Halo 16 has a relatively stable disk that does not fragment for ∼ 1.5 Myr.A key difference, however, is that the SMS in the primary disk collapses to a BH at about the time the second disk forms because it encounters the post-Newtonian instability at the very end of the main sequence at a mass of 109 kM , when its core hydrogen is exhausted.It is unclear from the Enzo simulation how X-rays from this black hole would affect the evolution of the companion disk but we discuss some possible outcomes in §4.2.As the primary disk continues to fragment, 3-body interactions fling the companion disk into a highly elliptical orbit ∼ 300 kyr after its formation.The initial burst of accretion is followed by a long quiescent phase (∼ 0.5 Myr) during which the second SMS, which is now ∼ 40 kM , becomes almost completely thermally relaxed and almost fully convective (Woods et al. 2020).The rapid accretion beginning at ∼ 750 kyr builds up a massive radiative envelope on top of this convective core and the mass of the star grows to 186 kM by the end of the simulation.At this point the SMS is still on the main sequence in a disk on a long orbit around the first disk, but it appears to be on the verge of collapse.With a mass nearing the upper limit for SMSs with similar accretion rates (e.g., Woods et al. 2017;Haemmerlé et al. 2018a;Woods et al. 2021), the star is unlikely to survive for much longer so the system may soon produce a DCBH binary.
As the primary disk in Halo 19 begins to fragment, two clumps merge and form a stable companion disk after ∼ 800 kyr.The initial burst of accretion due to the formation of this disk is particularly strong and builds up a star with a massive convective envelope in ∼ 500 kyr.Accretion in the disk is so rapid that the evolution of the second SMS overtakes that of the first.By the time the disks merge, the first SMS has reached a mass of 46 kM and a central hydrogen fraction of 0.28 but the second is at a mass of 133 kM and a central hydrogen fraction of 0.54.Both stars are still on the main sequence.
The primary disk in Halo 20 begins to fragment soon after formation and quickly forms a single, massive companion disk.Accretion rates in the disk are fairly high, and are rejuvenated at one point when a clump collides with the disk.These high rates quickly build up a SMS with a deep, high-entropy envelope.The companion disk later exchanges mass with the primary disk over a number of orbits that mostly halts accretion onto the second star and allows it to thermally relax.In the meantime, the SMS in the primary disk collapses at 1.48 Myr to a BH via the post-Newtonian/Chandrasekhar instability late in the main sequence at a core hydrogen fraction of 0.14 and a mass of 189 kM .At 1.77 Myr, the end of the simulation, the companion disk is in a relatively long elliptical orbit (∼ 1.8 pc separation) that is growing in radius.At this point, the companion is still on the main sequence but, as in Halo 16, appears to be about to collapse, having reached a mass of 154 kM .

DISCUSSION
In most of the cases above, the accretion history and evolution of the second SMS ends with the merger of its host disk with the primary disk.Although the subsequent evolution of the stars or DCBHs in the disks cannot be determined here, a number of outcomes can be inferred for them for follow-up with future simulations, as we outline below.

MS-MS Pairs
The primary and secondary disks in Halo 01, Halo 19, and nominally Halo 02 merge while their central stars are still on the main sequence (MS).With a maximum physical resolution of 0.014 pc (∼ 3 kAU), the Enzo simulations cannot determine if they later interact.Drag forces in the inner disk could bring about a swift merger between the two SMSs or the stars could carve out a gap in the circumbinary disk.Just before the merger of the two disks in Halo 02 one of the SMSs appears to be on the verge of collapse via the post-Newtonian instability.Examples such as this make it clear that a more careful treatment of stellar structure during the resulting enormous surge in accretion is essential (Menon & Heger 2017) and could be informed by smoothed-particle hydrodynamics calculations (e.g., Glebbeek et al. 2013).Such mergers may also produce unique observational transients that signal the formation of SMSs in the early Universe.If instead these supermassive stellar pairs evolve into relatively long-lived binaries, they could later become BH -MS or BH -BH pairs.

MS-BH Pairs
In Halos 08 and 12 the primary and companion disks merge after the star in the primary disk has collapsed to a DCBH while the SMS in the other disk is still on the main sequence.Given that the masses of both objects are comparable in both halos, later interactions between the two are unlikely to produce tidal disruption events.The tidal disruption radius in such cases is of the order of the stellar radius and would therefore only happen during a collision.However, fragmentation on scales smaller than those resolved here may lead to the formation of much less massive stars that could then be torn apart by the DCBH and produce a powerful radio or NIR transient (Kashiyama & Inayoshi 2016;Regős et al. 2021).
If these objects form a long-lived binary, they will continue to evolve subject to additional torques from the circumbinary disk.Depending on the orbital evolution, the companion SMS may overflow its Roche lobe while still on the main sequence or as it evolves, in which case its deep radiative envelope could sustain long-lived mass exchange and, perhaps, the formation of a "supermassive" X-ray binary (e.g., Soberman et al. 1997).Eventually, the collapse of the second SMS due to exhaustion of core nuclear fuel or the Chandrasekhar instability will produce a DCBH binary as discussed below.

BH-BH Pairs
Some of our models may produce long-period DCBH binaries.At the end of the simulation in Halo 10 the second SMS was still on the main sequence but its host disk was in a highly elliptical orbit with the first disk, whose star had collapsed to a DCBH.This long-period orbit suggests that the SMS will collapse to a DCBH before interacting closely with the other BH.Although the secondary SMSs in Halos 16 and 20 are also still on the main sequence at the end of the simulation they are both on the verge of collapse via the post-Newtonian instability unless accretion stops.Both of these halos could thus produce long-period (∼ 1 pc separation) binary DCBHs.The timescale for a merger between the two BHs would depend on their dynamics in the halo and drag forces due to gas flows therein, and the strength of their gravitational wave (GW) signal depends on their masses when they finally do merge.Here, we estimate the GW emission from the merger of the DCBH binary in Halo 20 (Robson et al. 2019) 1 , considering a sky and polarization averaged transfer function and assuming that the orbit has circularized by prior GW emission.The expected signal is shown in Figure 3 for the redshift at which the halo collapses, z = 17.7.This merger will be detectable by LISA, as would mergers of the other BH -BH pairs in our halos at their respective redshifts.

CONCLUSIONS
A growing number of recent numerical simulations indicate that the formation of multiple very massive stars (e.g., Wise et al. 2019;Regan et al. 2020) or SMSs (e.g., Latif et al. 2020;Patrick et al. 2020) was common in primordial, atomic-cooling halos, and that these could be the origin of the first quasars in the Universe.Here, we have shown that secondary disks in these halos form SMSs whose accretion histories can be coupled to those of the stars in the primary disks via tidal interactions.Most satellite disks merge with the primary disk in just 1 -2 Myr, which could yield supermassive stellar merg- ers, supermassive X-ray binaries, massive black hole seed mergers, and SMS binaries.Such interactions could pro-foundly affect the evolution of the stars, the SMBH seeds they produce, and their observational signatures.Simulations that follow the dynamics of these binaries and mass exchange between them with stellar evolution calculations will be the aim of future studies.

Figure 1 .
Figure 1.Accretion rates from Patrick et al. (2020) for the secondary stars.Note that times here are measured from the formation of the secondary disk, not the beginning of the simulation.

Figure 2 .
Figure 2. Kippenhahn diagrams of the internal structures of the stars in the satellite disks over time.The blue colors indicate energy generation rates from nuclear-burning, green lines denote convective regions, red lines mark semi-convective regions, and light blue lines indicate radiative/convectively neutral regions.

Figure 3 .
Figure 3. Sky-averaged LISA sensitivity curve (orange) and the observable characteristic strain (blue) from the merger of the DCBHs in Halo 20.
T.E.W. acknowledges support from the National Research Council Canada's Plaskett Fellowship.S.P. was supported by STFC grant ST/N504245/1 and D.J.W. was supported by the Ida Pfeiffer Professorship at the Institute of Astrophysics at the University of Vienna.A.H. was supported by the Australian Research Council (ARC) Centre of Excellence (CoE) for Gravitational Wave Discovery (OzGrav) through project number CE170100004, by the ARC CoE for All Sky Astrophysics in 3 Dimensions (ASTRO 3D) through project number CE170100013, and by the National Science Foundation under Grant No. PHY-1430152 (JINA Center for the Evolution of the Elements, JINA-CEE).