Chemically Homogeneous Evolution: A rapid population synthesis approach
Jeff Riley, Ilya Mandel, Pablo Marchant, Ellen Butler, Kaila Nathaniel, Coenraad Neijssel, Spencer Shortt, Alejandro Vigna-Gomez
MMNRAS , 1–15 (2020) Preprint 2 October 2020 Compiled using MNRAS L A TEX style file v3.0
Chemically Homogeneous Evolution: A rapid populationsynthesis approach
Jeff Riley , , Ilya Mandel , , , , Pablo Marchant , Ellen Butler ,Kaila Nathaniel , Coenraad Neijssel , , Spencer Shortt , Alejandro Vigna-G´omez School of Physics and Astronomy, Monash University, Clayton, Victoria 3800, Australia ARC Centre of Excellence for Gravitational Wave Discovery – OzGrav, Australia Birmingham Institute for Gravitational Wave Astronomy and School of Physics and Astronomy,University of Birmingham, B15 2TT, Birmingham, UK ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions – ASTRO 3D, Australia Institute of Astronomy, KU Leuven, Celestijnenlaan 200D, 3001 Leuven Argelander-Institut f¨ur Astronomie, Universit¨at Bonn, Auf dem H¨ugel 71, 53121 Bonn, Germany, Department of Mathematics, University of Colorado Boulder, Boulder, CO, USA DARK, Niels Bohr Institute, University of Copenhagen, Jagtvej 128, 2200, Copenhagen, Denmark
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
We explore chemically homogeneous evolution (CHE) as a formation channel for mas-sive merging binary black holes (BBHs). We develop methods to include CHE in arapid binary population synthesis code, Compact Object Mergers: Population Astro-physics and Statistics (COMPAS), which combines realistic models of binary evolutionwith cosmological models of the star-formation history of the Universe. For the firsttime, we simultaneously explore conventional isolated binary star evolution under thesame set of assumptions. This approach allows us to constrain population proper-ties and make simultaneous predictions about the gravitational-wave detection ratesof BBH mergers for the CHE and conventional formation channels. The overall massdistribution of detectable BBHs is consistent with existing gravitational-wave observa-tions. We find that the CHE channel may yield up to ∼
70% of all gravitational-wavedetections of BBH mergers coming from isolated binary evolution.
On September 14th, 2015 the first direct observation of grav-itational waves was made by the Advanced Laser Interfer-ometer Gravitational-wave Observatory (aLIGO) (Abbottet al. 2016). The detected signal, now known as GW150914,was also the first observation of two black holes merging,thus confirming the existence of binary stellar-mass blackhole systems and providing evidence that they can mergewithin the current age of the Universe. Based on 10 binaryblack hole (BBH) detections during the first two observingruns of aLIGO and advanced Virgo, Abbott et al. (2019)estimate a local BBH merger rate of 25 – 109 Gpc − yr − at90% confidence.How the BBH sources of these gravitational wave sig-nals form remains an open question. To be the source ofgravitational waves detected at aLIGO, which is sensitiveto signals with frequencies of tens to hundreds of Hz, com-pact objects orbiting each other must spiral in as they loseenergy through the emission of gravitational waves. Orbitalenergy loss through gravitational wave emission is not ef-ficient at wide separations, and the timescale for gravita-tional wave emission to drive a binary to merger scales as the fourth power of the orbital separation (Peters 1964).In order for two 30 M (cid:12) black holes to merge within ≈ (cid:46) R (cid:12) . And therein lies a problem, as thisis smaller than the radial extent reached by typical slowlyrotating massive stars during their evolution. The differentastrophysical channels proposed for forming merging BBHsgenerally fall into two categories (see, e.g., Mandel & Farmer2018; Mapelli 2018, for reviews):(i) isolated binary evolution, in which two stars may in-teract through tides and mass transfer, but are dy-namically decoupled from other stars (e.g., Tutukov& Yungelson 1973, 1993; van den Heuvel 1976).(ii) dynamical formation, where dynamical interactions ina dense environment and/or a hierarchical triple sys-tem play a key role in forming and hardening a com-pact BBH (Sigurdsson & Hernquist 1993; Miller &Lauburg 2009; Rodriguez et al. 2015; Antonini et al.2016; Bartos et al. 2017; Stone et al. 2017).A variant of the isolated binary evolution channel re-lies on rotationally-induced chemical mixing in massive starsto prevent the establishment of a strong chemical gradient c (cid:13) a r X i v : . [ a s t r o - ph . S R ] S e p J. Riley et al. (Maeder 1987; Maeder & Meynet 2000; Heger et al. 2000).As long as the star continues to rotate at a sufficiently highrate, it will remain quasi-chemically homogeneous (Maeder1987; Langer 1992). Contrary to the core-envelope struc-ture exhibited by conventional, more slowly rotating stars,and the characteristic expansion of the envelope as the corecontracts, the radius of quasi-chemically homogeneous starswill shrink or remain constant as they become hotter andbrighter (Yoon et al. 2006; Mandel & de Mink 2016).As we discuss below, previous work on the chemicallyhomogeneous evolution (CHE) channel for BBH formation(Mandel & de Mink 2016; Marchant et al. 2016; de Mink &Mandel 2016; du Buisson et al. 2020) explored this channelindependently of the usual isolated binary evolution chan-nel. In this paper, we present our rapid population synthesismodel for the CHE of binary systems, allowing for a directcomparison of the rates and properties of CHE and non-CHE BBHs under the same set of assumptions. Our CHEmodel is implemented in the rapid binary population syn-thesis code COMPAS (Stevenson et al. 2017; Vigna-G´omezet al. 2018), with thresholds on CHE evolution computed us-ing rotating stellar models in the Modules for Experimentsin Astrophysics (MESA) code (Paxton et al. 2015; Paxtonet al. 2013, 2011).The remainder of this paper is organised as follows. Sec-tion 2 is a brief outline of CHE and previous work on theformation of BBHs through this channel. Section 3 presentsa description of our CHE model and the implementation ofthe model in COMPAS. We present our results in Section 4.We provide some concluding remarks in Section 5.
Stars evolving on the main sequence (MS) typically developincreasingly helium-rich cores and hydrogen-rich envelopesas radial mixing is inefficient. However, von Zeipel (1924)showed that rotating stars cannot simultaneously be in hy-drostatic and thermal equilibrium if the rotational velocityis a function of radius only, which has been argued to resultin meridional currents in the radiative layers of a rotatingstar (Sweet 1950; Eddington 1929). In massive rapidly rotat-ing stars in low-metallicity environments, these currents canmix material from the convective core throughout the radia-tive envelope, leading to chemically homogeneous evolutionfor rapidly rotating stars (Maeder 1987).Due to strong chemical mixing, chemically homoge-neous stars do not maintain a hydrogen-rich envelope – thusavoiding the dramatic expansion exhibited during the post-main sequence phase by non-chemically homogeneous stars.The radius of a chemically homogeneous star remains stable,or shrinks slowly, as the star becomes increasingly heliumrich over the course of the main sequence, with the star con-tracting to a massive naked helium star post-main sequence.Chemically homogeneous components of a very close binarysystem can thus avoid overfilling their Roche lobes, masstransfer, and probable merger.de Mink et al. (2009) modelled the evolution of rotatingmassive stars using the hydrodynamic stellar evolution codedescribed by Yoon et al. (2006) and Petrovic et al. (2005),which includes the effects of rotation on the stellar structure and the transport of angular momentum via rotationally-induced hydrodynamic instabilities (Heger et al. 2000). Thebinary models developed by de Mink et al. (2009) and Songet al. (2016) show that constituent stars in very tight binarysystems can achieve rotational frequencies sufficient to in-duce CHE. de Mink et al. (2009) proposed CHE as a viableformation channel for high-mass black-hole X-ray binaries.VFTS 352 (Almeida et al. 2015) and HD 5980 (Koenigs-berger et al. 2014) are examples of observed binary systemsthought to have undergone CHE (de Mink & Mandel 2016).Mandel & de Mink (2016) and Marchant et al. (2016)introduced and investigated CHE as a channel for formingmerging BBHs. They concluded that for sufficiently highmasses and sufficiently low metallicities, a narrow range ofinitial orbital periods (short enough to allow rapid rota-tion necessary for CHE, but not so short that the binarywould immediately merge) could allow this channel to pro-duce merging BBHs.Mandel & de Mink (2016); de Mink & Mandel (2016)used approximate thresholds for CHE based on the models ofYoon et al. (2006) to investigate the rates and properties ofBBHs formed through the CHE channel. They estimated amerger rate of ∼
10 Gpc − yr − in the local Universe for thischannel, subject to a number of evolutionary uncertainties,which they explored in a population-synthesis-style study.Marchant et al. (2016) used the MESA code to conductdetailed simulations of the CHE channel, which were fol-lowed until the BBH stage. The simulations were conductedfor close binaries with component masses above ∼
20 M (cid:12) ,and included the over-contact phase in a majority of CHEBBH progenitors. Marchant et al. (2016) suggested that aslong as material does not overflow the L2 point in over-contact binaries, co-rotation can be maintained, and a spiral-in due to viscous drag can be avoided. In this scenario, closebinary systems typically enter the over-contact phase in theearly stages of core hydrogen burning, and then equilibratetheir masses through mass transfer between the constituentstars. du Buisson et al. (2020) extended the results of theMESA simulations performed by Marchant et al. (2016) andcombined them with the cosmological simulations of thechemical and star-formation history on the universe by Tay-lor & Kobayashi (2015). Their population synthesis study in-vestigated the population properties, cosmological rates andaLIGO detection rates of BBHs, including the dependenceon the early-Universe star formation rate (SFR), which theyfind to be mild for moderate variations in the high-redshiftSFR.
In this section, we describe the implementation of CHEwithin the COMPAS rapid binary population synthesis code.Using COMPAS allows us to rapidly evolve a large syntheticpopulation of binaries, which includes binaries whose compo-nent stars evolve conventionally (i.e. along a redwards trackon the Hertzsprung-Russell diagram (HRD)), and otherswhose components evolve via CHE (i.e. along a bluewardstrack on the HRD), thus providing data for both pathwaysthat can be compared directly. Below, we summarise thekey physics implemented in COMPAS, starting with our ap-proximate model of quasi-chemically homogeneous evolution
MNRAS , 1–15 (2020) hemically Homogeneous Evolution: A rapid population synthesis approach based on MESA experiments, as well as the choices made forthe metallicity-specific star formation history. The basics of stellar and binary evolution and BBH popula-tion modelling in COMPAS are described by Stevenson et al.(2017); Vigna-G´omez et al. (2018); Neijssel et al. (2019).Here, we provide a brief summary and describe differencesfrom previous COMPAS studies.
We used a set of MESA models of single stars with a fixedrotational frequency and no mass loss to determine the min-imal angular frequency ω necessary for chemically homoge-neous evolution as a function of mass and metallicity. Ourfits to these angular frequency thresholds are provided inAppendix A. COMPAS uses these fits to determine whethera star is evolving chemically homogeneously.Stellar evolution in COMPAS follows the analytical fitsof Hurley et al. (2000) to the stellar models from Polset al. (1998). In order to address CHE, we introduce a newchemically homogeneous (CH) stellar type to the Hurleyet al. (2000) collection of stellar types. In our simplifiedmodel, we neglect the very limited radial evolution of aCH star and set its radius equal to the zero-age main se-quence (ZAMS) radius of a non-rotating star of the samemass and metallicity. This slightly over-estimates the radiusrelative to MESA models of CH stars, and therefore con-servatively narrows the parameter space for CHE. We com-pute the mass loss rate for CH stars in the same way as forregular MS stars, but with this fixed rather than evolvingradius. As a consequence, the total mass lost over the MS byCH stars is generally within (cid:46)
10% of that lost by non-CHstars of the same ZAMS mass and metallicity, except for themost massive stars in our simulations, with initial massesabove 100 M (cid:12) , where the absence of radial expansion leadsto significantly reduced MS mass loss estimates for CH stars.Finally, we assume that if a star evolves chemically homo-geneously through the main sequence, it contracts directlyinto a naked helium star at the end of the main sequence,retaining its full mass at that point. Future evolution followsthe Hurley et al. (2000) models of helium stars.Tides are very efficient at ensuring circularisation andsynchronisation in very close binaries through tidal locking(e.g., Hut 1981). We therefore assume that all potential can-didates for CHE are tidally synchronised at birth, so thattheir rotational angular frequency equals the orbital angu-lar frequency. We check this angular frequency at birth todetermine whether a star belongs to the CH type and con-tinue to check it at every time step on the main sequence. Ifthe angular frequency ever drops below the threshold valuefor CHE, e.g., because of binary widening as a consequenceof mass loss through winds, the star is henceforth evolvedas a regular main sequence star (in our simplified treat-ment, it immediately jumps to the track of a regular mainsequence star of the same mass). We assume that once achemical gradient is formed, it is very challenging to over-come and ensure efficient mixing, so in our model, a star thatis not evolving chemically homogeneously cannot become a CH star (cf. BPASS models, which allow quasi-chemicallyhomogeneous evolution through accretion-induced spin-upEldridge et al. 2017). Although we assume perfect tidal syn-chronisation for CH stars, we disregard the angular momen-tum stored in the stellar rotation when considering binaryevolution with mass loss. Each binary system in a COMPAS simulation is describedat birth (i.e. at zero-age main sequence (ZAMS)) by its ini-tial conditions: constituent star masses, separation, eccen-tricity and metallicity. Initial conditions for our experimentswere chosen using statistical distribution functions from theliterature that were themselves based on observations. Wedescribe the most important of these, and some importantparameters that affect the evolution of the constituent starsas well as the binary system, in the following paragraphs.The mass of the primary star in the binary system (themore massive star at ZAMS) m , i is described by the Kroupa(2001) initial mass function (IMF), the distribution functionof which is given by p ( m ,i ) ∝ m − α ,i , (1)where α = 2 . m , i ∈ [5 , (cid:12) . We assume that the IMF is the same forall metallicities.The mass of the secondary star (less massive at ZAMS) m , i is determined by drawing a mass ratio between theconstituent stars q i ≡ m , i / m , i that follows a flat distri-bution, p ( q i ) = 1 (Sana et al. 2012; Kobulnicky et al. 2014).Since we are interested in BBH formation, we explore only m , i ≥ (cid:12) here. However, for both the primary and sec-ondary mass, we consider the full mass range to normalisethe simulation results to a given star-forming mass or starformation rate (e.g., Neijssel et al. 2019).The initial separation is drawn from a flat-in-log distri-bution independently of the masses (see Moe & Di Stefano2017 for coupled initial conditions): p ( a i ) ∝ a i , (2)where a i ∈ [0 . , e i = 0),see Vigna-G´omez et al. (2020) for further discussion. Closebinaries are tidally circularised at birth, so this has no im-pact on potential CHE systems.We simulate thirty different metallicities spaced uni-formly in the logarithm across the range − ≤ log Z ≤− . We use the mass loss rates as prescribed by Hurley et al.(2000, 2002) and references therein for cooler stars withtemperatures of 12 ,
500 K and below. For stars hotter than12 ,
500 K we use the wind mass loss rates from Vink et al.(2001), as implemented in Belczynski et al. (2010).The luminous blue variable (LBV) stars (Maeder 1989;Pasquali et al. 1997), located close to the Humphreys-Davidson limit in the HRD (Humphreys & Davidson 1979),
MNRAS000
MNRAS000 , 1–15 (2020)
J. Riley et al. are treated differently. For these stars we use the LBV windmass loss rate prescribed by Belczynski et al. (2010): dMdt = f LBV × − M (cid:12) yr − , (3)where f LBV = 1 . dMdt = f WR × − L . ( Z Z (cid:12) ) m M (cid:12) yr − , (4)where L is the luminosity, m = 0 .
86 (Vink & de Koter 2005),we take Z (cid:12) = 0 .
014 (Asplund et al. 2009) and f WR = 1 . f WR in order to study the impact ofWR mass loss on CHE: f WR ∈ { . , . , . , . } .All mass lost in winds is assumed to promptly departthe binary without further interaction with the companionin so-called ‘Jeans mode’ mass loss, carrying away the spe-cific angular momentum of the donor. We use the prescriptions described in Vigna-G´omez et al.(2018); Neijssel et al. (2019) to determine the dynami-cal stability of mass transfer through Roche-lobe overflow(RLOF), the fraction of mass accreted onto the compan-ion and the specific angular momentum carried away bynon-conservative dynamically stable mass transfer, and theoutcome of common-envelope evolution. For non-CHE bi-naries that go through a common-envelope phase, we as-sume that Hertzsprung-gap donors do not survive (the ’pes-simistic’ prescription of Belczynski et al. 2007; Neijssel et al.2019) and we assume that immediate post-common-envelopeRLOF indicates a merger.We deviate from previous COMPAS models in the treat-ment of binaries that experience RLOF at ZAMS. Unlikeprevious work, we now allow such binaries to equilibratetheir masses. The new separation of the equal-mass binarywith a conserved total mass is determined by angular mo-mentum conservation. Binary components are allowed toover-fill their Roche lobes, creating over-contact systems.However, if the components extend past the L2 Lagrangepoints after equilibration, we assume that the binary losesco-rotation and promptly merges (Marchant et al. 2016). Forequal-mass circular binaries, the volume-equivalent radiusfor half of the volume within the L2 equipotential surfaceequals half the orbital separation. Therefore, our criterionfor avoiding a prompt merger is equivalent to demandingthat the sum of the unperturbed stellar radii is smaller thanthe orbital separation a . Stellar evolution models predict that single stars with he-lium cores in the range ∼
60 – 130 M (cid:12) can become unstabledue to electron-positron pair production, leading to pair-instability supernovae (PISNe) which disrupt the star, leav-ing no remnant behind (e.g., Fowler & Hoyle 1964; Barkatet al. 1967; Fraley 1968; Woosley et al. 2002; Woosley 2019;Farmer et al. 2019). Stars with helium cores more mas-sive than 130 M (cid:12) also experience a rapid collapse drivenby pair production, but in these stars photodisintegrationprevents a subsequent explosion; such stars may again pro-duce merging BBHs (e.g., Marchant et al. 2016; du Buissonet al. 2020), but are not explored in our models, which havemaximum initial stellar masses of 150 M (cid:12) . Meanwhile, starswith somewhat lower helium core masses, between ∼
35 and ∼
60 M (cid:12) , are predicted to eject significant fractions of theirtotal mass over several episodes (e.g., Yoshida et al. 2016;Woosley 2017; Marchant et al. 2019; Renzo et al. 2020).Such pulsational pair-instability supernovae (PPISNe) leavebehind a black hole remnant, albeit with a reduced mass.PISNe and PPISNe are expected to produce a PISN massgap in the distribution of remnant masses from single stel-lar evolution – a dearth of black holes with masses between ∼ M (cid:12) and ∼ M (cid:12) .Some superluminous supernovae have been identified asPISN candidates (Gal-Yam 2012, and references therein),while iPTF2014hls has been identified as a PPISN candi-date (Arcavi et al. 2017). Furthermore, the distribution ofmasses of gravitational-wave observations appeared consis-tent with a cutoff due to (P)PISNe (Abbott et al. 2019),though GW190521 is a BBH merger with at least one com-ponent in the predicted PISN mass gap (Abbott et al. 2020).Here, we follow the Stevenson et al. (2019) fit to theMarchant et al. (2019) models for predicting the range ofPISN masses and the PPISN remnant masses from themasses of the progenitor helium cores. We apply the entirePPISN mass loss in one time step. Moreover, in our treat-ment both supernovae happen in one timestep for equal-mass stars. This over-estimates the post-supernova periodand eccentricity of binaries whose components lose signifi-cant mass in a PPISN. The local merger rate of BBHs depends on their formationrate at higher redshifts due to the possibly significant timedelays between formation and merger, and is therefore sensi-tive to the star formation rate as a function of redshift. Fur-thermore, the yield of BBHs per unit star forming mass, theBBH mass distribution, and the distribution of delay timesbetween formation and merger are all sensitive functionsof the metallicity of progenitor stars, both for CHE (e.g.Marchant et al. 2016) and non-CHE (Neijssel et al. 2019;Chruslinska et al. 2019) systems. We must therefore specifya metallicity-specific star formation rate (MSSFR) in orderto estimate the merger rate and properties of BBHs. We usethe preferred model of Neijssel et al. (2019) for the MSSFR.Figure 1 shows the contribution of different ranges of starformation metallicities to the total star formation rate. Thismodel has higher star formation metallicities in the local
MNRAS , 1–15 (2020) hemically Homogeneous Evolution: A rapid population synthesis approach Universe than the Taylor & Kobayashi (2015) model usedby du Buisson et al. (2020) (cf. their Figure 2).
Figure 1.
The total star formation rate as a function ofredshift (red) and subdivided into different ranges of metal-licity, following the preferred model of Neijssel et al. (2019).The dark blue and green curves are most relevant for BBHformation.
We evolved a total of 12 million binaries as described insection 3. These were equally divided into 30 metallicity binsand 4 choices of the WR mass loss rate multipliers f WR , fora total of 100 ,
000 binaries for each of 120 combinations of Z and f WR . Binaries are evolved until a double compactobject is formed, or until an event happens which makesthis outcome impossible (e.g., the stars merge or the binarybecomes unbound), or the system reaches 14 Gyr in age.Our simulations are based on a Monte Carlo samplingof binaries. We estimate the sampling uncertainty on all de-rived quantities via bootstrapping: we uniformly resample,with replacement, a new population of 12 million binariesfrom the original, evolved, population of 12 million. Errorbars on plots, where shown, correspond to the 5 th and 95 th percentiles from bootstrapping. The population statistics are shown in Table 1. From a pop-ulation of 11,025,296 binaries that survived beyond ZAMS(i.e. did not merge at ZAMS), 16,891 were composed of twoCH stars at ZAMS, with a further 10,419 composed of oneCH star and one MS star at ZAMS. Furthermore, in all of thebinaries with only one CH star at ZAMS it was, as we wouldexpect, the primary, more massive, star that was chemicallyhomogeneous. A total of 261,741 BBHs were formed in thesimulation, but only 43,625 of these were close enough tomerge within 14 Gyr, the current age of the Universe. Amongthe 13,644 simulated binaries that evolved chemically ho-mogeneously throughout the main sequence, 9,670 went onto form BBHs, the vast majority of which, 9,062, merged within 14 Gyr (the few non-merging ones are those whichlost significant mass in PPISNe).
Figure 2 presents a visual summary of the evolutionary out-comes for each of the 12 million binary systems synthesised,with each point on the plot representing a single binary sys-tem, and the colour indicating the initial parameters andthe outcome of the evolution (per the legend). We are par-ticularly interested in systems for which both stars evolvechemically homogeneously and eventually collapse to forma BBH, so we have agglomerated some of the less interest-ing progenitor types and outcomes into groups so that theplot is not overly busy. Because Figure 2 is a summary overthe entire grid of metallicities and WR mass loss rate mul-tipliers synthesised, it allows us to see on a broad scale theevolutionary outcomes for both CHE systems and non-CHEsystems. The COMPAS models for the formation of non-CHE BBHs have been discussed by (Neijssel et al. 2019), sowe will focus our discussions hereafter one the CHE channel.Figure 3 shows the parameter space in which CHE is ex-pected to occur in synchronously rotating binaries accordingto our CHE threshold. The darkest grey area in the lowerpart of the diagram indicates the region in which L2 overflowoccurs and the stellar components merge; the lighter greyarea in the upper part indicates the region in which the stel-lar components do not rotate rapidly enough to induce CHE.The central, lighter, area of the diagram indicates the regionin which we expect CHE to occur, with the darker, lower,part of the central area indicating the important region ofbinaries whose components overflow their Roche lobes butavoid L2 overflow, occupied by the over-contact systems de-scribed by Marchant et al. (2016). This over-contact regionis responsible for much of the BBH formation through CHE(cf. Figure 2).As expected (given our PPISN and PISN mass limits,see section 3.1.5), we see BBHs from PPISNe begin to ap-pear at a total ZAMS mass of (cid:38)
70 M (cid:12) while PISN eventsappear at a total ZAMS mass of (cid:38)
120 M (cid:12) . A few unboundCHE systems correspond to simultaneous PISNe that in-stantaneously removed more than half the mass of the bi-nary in our treatment (see Section 3.1.5). In practice, suchsystems will undergo a series of pulsations leading to non-simultaneous mass loss and may survive, but at separationstoo large to merge within the current age of the Universe.The horizontal band of PISNe just above the CH binaries inFigure 2 are hybrid systems comprised of a CH star and aMS star, whereas the vertical band of PISNe at the upperright of the plot are systems comprised of two MS stars.
The initial system total masses and orbital periods of CHEsystems that go on to form BBHs merging within 14 Gyr areshown in Figure 4. We show binaries evolved with the WRmass loss multiplier f WR = 1. Each point on the plot rep-resents a simulated binary shaded according to its metallic-ity. Higher metallicity binaries are shifted toward the top ofthe plot. This is consistent with Figure 3, which shows thathigher-metallicity stars have greater stellar radii and hence MNRAS000
The initial system total masses and orbital periods of CHEsystems that go on to form BBHs merging within 14 Gyr areshown in Figure 4. We show binaries evolved with the WRmass loss multiplier f WR = 1. Each point on the plot rep-resents a simulated binary shaded according to its metallic-ity. Higher metallicity binaries are shifted toward the top ofthe plot. This is consistent with Figure 3, which shows thathigher-metallicity stars have greater stellar radii and hence MNRAS000 , 1–15 (2020)
J. Riley et al.
Table 1.
Population Statistics Population f WR = 0 f WR = 0 . f WR = 0 . f WR = 1 . Total
Number of binaries evolved 3,000,000 3,000,000 3,000,000 3,000,000 12,000,000L2 overflow at ZAMS 243,717 243,638 243,419 243,930 974,704Surviving binaries 2,756,283 2,756,362 2,756,581 2,756,070 11,025,296Surviving PopulationAt least one star experiencing RLOF at ZAMS 75,530 75,410 75,328 75,707 301,975Both stars in binary CH at ZAMS 4,193 4,281 4,201 4,216 16,891Primary only CH at ZAMS 2,607 2,615 2,593 2,604 10,419Secondary only CH at ZAMS 0 0 0 0 0Post-ZAMS Merger 618,001 616,749 618,425 620,250 2,473,425BBHs formed 68,231 67,200 66,016 60,294 261,741BBHs merging in 14 Gyr 11,004 11,048 10,926 10,647 43,625Both stars CH at ZAMSAt least one star experiencing RLOF at ZAMS 3,661 3,761 3,715 3,715 14,852Both stars remained CH on MS 3,444 3,461 3,379 3,360 13,644Primary only remained CH on MS 43 89 116 160 408Secondary only remained CH on MS 0 0 0 0 0Neither star remained CH on MS 706 731 706 696 2,839BBHs formed 2,152 2,370 2,527 2,621 9,670BBHs merging in 14 Gyr 2,057 2,322 2,377 2,306 9,062Primary only CH at ZAMSAt least one star experiencing RLOF at ZAMS 0 0 0 0 0Primary remained CH on MS 1,405 1,353 1,341 1,337 5,436BBHs formed 0 2 3 7 12BBHs merging in 14 Gyr 0 0 0 0 0 greater minimal separation, as well as lower CHE thresholdrotational frequency.Binaries with reduced WR winds have similar initialdistributions, but show a clear-cut maximum total mass of ≈ M (cid:12) , which matches the mass threshold of 60 M (cid:12) forindividual He star masses beyond which PISNe occur andleave no remnants. At higher f WR , high-metallicity systemscan lose a significant fraction of their mass, so binaries withinitial total masses above 120 M (cid:12) can avoid PISNe.To illustrate this, we plot the mass lost by a CH starwith a ZAMS mass of 40 . M (cid:12) over the naked helium phasein Figure 5, for a range of WR mass loss multipliers andmetallicities. At f WR = 1 and Z = Z (cid:12) , this star loses nearlyhalf of its mass in WR winds. Meanwhile, at low metallici-ties, which are typical for high formation redshifts, the totalmass lost in WR winds is very low, except at artificiallyenhanced f WR values of 5 and 10, which disagree with ob-servational constraints and are not considered in this study.Consequently, we do not expect to see a significant impactof f WR on low-metallicity BBH formation, which matchesour findings as discussed below.Table 1 shows, across all simulated metallicities and WRmass loss multipliers, ∼
80% of binaries composed of twoCH stars at ZAMS retain two CH stars at the end of themain sequence. For binaries composed of one CH star andone MS star at ZAMS, the CH star will remain chemicallyhomogeneous by the end of the MS in only ∼
50% of simu-lations. Since we assume tidal locking in the CHE model im-plemented in COMPAS, as a binary widens due to mass lossand the orbital frequency of the binary slows, the rotationalfrequency of the constituent CH stars slows commensurably.Binaries in which only the primary is CH at ZAMS avoidedRLOF and are typically wider, so further widening through winds is more likely to spin down the primary sufficiently toevolve off the CHE track.Figure 6 shows the distribution of the BBH total massesand orbital periods just after BBH formation for systemsevolving through the CHE channel. As in Figure 4, we selectonly BBHs that will merge in 14 Gyr and shade binaries bymetallicity. On this plot, we select f WR = 0 .
2. This allowsus to show not only the sharp disappearance of BBHs withtotal masses above ≈ M (cid:12) due to PPISN mass loss andcomplete disruption in PISNe, but also their reappearance atmasses above ≈ M (cid:12) , on the other side of the ‘PISN massgap’. There are only very few such high-mass binaries in oursimulations because, with our ZAMS mass upper limit of150 M (cid:12) , they require very low mass loss. Consequently, thereare no such binaries in our f WR = 1 . Figure 7 shows the merging BBH yield: the formation rateper unit star forming mass as a function of metallicity forsystems that will merge within 14 Gyr. The solid lines arethe rates for the entire population — both CHE and non-
MNRAS , 1–15 (2020) hemically Homogeneous Evolution: A rapid population synthesis approach Figure 2.
Initial parameters and final outcomes for each of the binary systems synthesised, showing the initial orbital period T ZAMS (in days) vs the initial total mass (in M (cid:12) ). The population represents a grid of 30 metallicities evenly spaced over the range-4 ≤ log ( Z ) ≤ -1.825, and 4 WR mass loss multipliers, f WR ∈ { . , . , . , . } . Regions shaded in black represent all systemsthat experienced L2 overflow at ZAMS; pale green represents systems that did not form BBHs. Systems for which both stars were chem-ically homogeneous at ZAMS and remained so throughout their main sequence lifetime are represented by regions shaded cyan if theyformed BBHs via regular core-collapse supernovae, blue if they formed BBHs after undergoing PPISNe, and magenta if they explodedas PISNe. Systems in which at least one of the stars did not evolve chemically homogeneously for its entire main sequence lifetime arerepresented by areas shaded light green if they formed BBHs via core-collapse supernovae, dark green if they formed BBHs followingPPISNe, and maroon if either star exploded as a PISN. CHE binaries — while the dashed lines are the rates for CHEbinaries only. WR mass loss multipliers are differentiated bythe colour of the lines.The overall yield of merging BBHs is quantitatively sim-ilar to the simulations of Neijssel et al. (2019), who predicteda yield of ∼
6, 4, and 1 merging BBHs per 10 M (cid:12) of starformation at Z = 0 . Z (cid:12) , 0 . Z (cid:12) , and 0 . Z (cid:12) , respectively.The small differences are due partly to the inclusion of theCHE channel as well as PISNe and PPISNe in this work,which were not included in Neijssel et al. (2019).Meanwhile, the low-metallicity CHE channel yield ofslightly less than 1 merging BBHs per 10 M (cid:12) of star for-mation is similar to both the Mandel & de Mink (2016) back-of-the-envelope estimate and the Marchant et al. (2016) de-tailed models which indicate ∼ . M (cid:12) of star formation at Z = 0 . Z (cid:12) .The paucity of CHE BBHs at high metallicity, Z (cid:38) . Z (cid:12) , is due primarily to a combination of the up-ward shifting of the allowed initial periods at higher metal-licities (see Figures 3 and 4) and greater orbital wideningby stronger high-metallicity winds. The increase in orbitalperiod at BBH formation increases the delay times, prevent-ing the BBHs from merging within 14 Gyr. The wideningby mass loss is ameliorated by reduced WR mass loss rates. However, the WR mass loss multipliers have negligible effectat low metallicities because the total mass loss rate is toolow even for f WR = 1 (see Figure 5 and associated discus-sion). Neijssel et al. (2019) discuss the impact of metallicityon the non-CHE BBHs yield, highlighting the contributionsof wind-driven widening and stellar evolutionary stage atmass transfer.Figure 8 shows the BBH formation rate per unit co-moving volume per unit source time as a function of red-shift. The formation rate for CHE BBHs peaks at z ≈ . f WR = 1 .
0, and at z ≈ . Figure 9 indicates the distribution of delay times betweenstar formation and BBH mergers. This figure combines allmetallicities with equal weights, without considering their
MNRAS000
MNRAS000 , 1–15 (2020)
J. Riley et al.
Figure 3.
Parameter space for equal-mass binary systems with the indicated companion mass at which chemically homogeneous evolutionis expected to occur at ZAMS. Solid lines show the thresholds for CHE implemented in COMPAS (see Appendix A), dotted lines areRLOF thresholds, and dashed lines are L2 overflow thresholds. Colours differentiate metallicities. Shading corresponds to Z = 0 . Figure 4.
Initial total masses and orbital periods for CHE sys-tems that go on to form BBHs that will merge within 14 Gyr.Each point represents one simulated binary, evolved with WRmass loss multiplier f WR = 1, shaded according to its metallicity. contribution to the observable systems, so should be viewedas an indicative sketch.Non-CHE binaries in Figure 9 have a very broad distri-bution of delay times. Some are very short, less than 10 Myr,due to significant hardening during mass transfer episodes,including through dynamically unstable mass transfer andcommon-envelope ejection, as well as fortuitously directed Figure 5.
Total mass lost by a WR star with a ZAMS mass of40 . (cid:12) as a function of metallicity. Line colour indicates the WRmass loss rate multiplier (solid lines). Also shown are the mass(on the same scale as the mass loss curves) and luminosity at thestart of the WR phase as a function of metallicity (dashed lines). supernova natal kicks. Meanwhile, there is an almost flattail of long delay times on this logarithmic plot, correspond-ing to a p ( τ delay ) ∼ /τ delay distribution.On the other hand, binaries formed through CHE areseen to have a more strongly clustered delay time distribu-tion, with typical delay times of between 100 Myr and 1 Gyr. MNRAS , 1–15 (2020) hemically Homogeneous Evolution: A rapid population synthesis approach Figure 6.
Total masses and orbital periods immediately afterBBH formation for CHE systems that will merge within 14 Gyr.Each point represents a simulated binary, evolved with WR massloss multiplier f WR = 0 .
2, shaded according to its metallicity. Theempty area between ≈
80 M (cid:12) and ≈
250 M (cid:12) is a consequenceof systems that lost mass as PPISNe or left no remnants afterexploding as PISNe.
Figure 7.
Yield of BBHs that will merge within 14 Gyr per unitstar forming mass as a function of metallicity. The solid lines arethe rates for the entire population — both CHE and non-CHEbinaries — while the dashed lines are the rates for the CHE bina-ries only. Colours indicate WR mass loss rate multipliers. Errorbars indicate 90% confidence intervals from sampling uncertainty.
There are no ultra-short delay times because, with the ex-ception of RLOF at ZAMS, such binaries do not undergomass transfer that could harden the binary. Moreover, thehigh masses of CHE stars imply that they do not experi-ence asymmetric supernovae and associated natal kicks inthe
COMPAS model.
Figure 8.
BBH formation rate per Gpc of comoving volume peryear as a function of redshift for BBHs that will merge within 14Gyr. Error bars indicate sampling uncertainty. Figure 9.
Distribution of delay times between formation andmerger for BBHs. All metallicities from the simulation are com-bined with equal weights and arbitrary counts per uniform binsin log delay time are shown. Error bars indicate sampling uncer-tainty.
The smallest time delay between formation and mergerfor CHE systems in our simulations ranges from ∼ .
025 Gyrfor f WR = 0 . ∼ .
033 Gyr for f WR = 1 .
0. The combina-tion of lower metallicities and reduced mass loss rates yieldsthe shortest delay times, allowing binaries to start evolutionfrom closer separations while avoiding L2 overflow and toavoid subsequent widening through mass loss. This is con-sistent with the minimal delay times found in other studies.Mandel & de Mink (2016), who consider only Z = 0 . Z (cid:12) ,estimate minimum delay times of ∼ . ∼ . MNRAS000
0. The combina-tion of lower metallicities and reduced mass loss rates yieldsthe shortest delay times, allowing binaries to start evolutionfrom closer separations while avoiding L2 overflow and toavoid subsequent widening through mass loss. This is con-sistent with the minimal delay times found in other studies.Mandel & de Mink (2016), who consider only Z = 0 . Z (cid:12) ,estimate minimum delay times of ∼ . ∼ . MNRAS000 , 1–15 (2020) J. Riley et al.
Figure 10.
BBH merger rate per Gpc of comoving volume peryear of source time as a function of redshift. Error bars indicatesampling uncertainty. the metallicity dependence. du Buisson et al. (2020) considerthe lowest metallicities among these studies, Z = 10 − , andfind the shortest delay times, ∼ .
02 Gyr.Some CHE binaries will be significantly widened bymass loss, potentially losing up to a factor of ∼ a M − (Peters 1964), soa factor of 2 each in mass decrease and semi-major axis in-crease would yield a factor of 2 ∼
100 increase in the delaytime. This explains the long delay time tail of the CHE BBHdistribution, as well as the decrease in the prominence of thistail as the WR wind mass loss multiplier is reduced. Evenwhen f WR = 0, some CHE BBHs will have long delay timesdue to the mass lost in PPISNe. Figure 10 shows the BBH merger rate per Gpc of comov-ing volume per year of source time as a function of red-shift. The merger rate for CHE BBHs peaks at z ≈ f WR = 1 .
0, and at z ≈ z ≈ f WR = 0 CHE BBHs, which explainstheir suppressed merger rate in the local Universe.The merger rates of BBHs that could be observed byaLIGO operating at final design sensitivity merger ratesare shown in Figure 11. Binaries formed through CHEhave higher average masses than non-CHE binaries, whichincreases the range within which they are detectable by Figure 11.
The merger rate of BBHs detectable by aLIGO atfinal design sensitivity, as a function of merger redshift. Errorbars indicate sampling uncertainty.
Figure 12.
Cumulative BBH detections as a function of mergerredshift, per year of observing at aLIGO O1 sensitivity. Error barsindicate sampling uncertainty. aLIGO. Therefore, CHE BBHs make up a higher fractionof all detections at greater redshifts.
Figures 12 and 13 show the predicted cumulative detectionrates per year of observing time as a function of redshift foraLIGO O1 and final design sensitivities, respectively.Figure 12 shows that the total expected detection rateat O1 sensitivity is 38–55 detections per year, depending onthe assumed value of the WR mass loss rate multiplier. Thiswould correspond to 17–25 detections over the 166 days ofcoincident data over the first two advanced detector observ-
MNRAS , 1–15 (2020) hemically Homogeneous Evolution: A rapid population synthesis approach Figure 13.
Cumulative BBH detections as a function of mergerredshift, per year of observing at aLIGO final design sensitivity.Error bars indicate sampling uncertainty. ing runs, O1 and O2. In fact, only 10 BBHs were observedduring this time (Abbott et al. 2019).The increased detection rate relative to the preferredmetallicity-specific star formation rate model of Neijsselet al. (2019), who predicted 22 detections per year, in agree-ment with the O1 and O2 observations, is due to the con-tribution of CHE BBHs. CHE BBHs may constitute up to ∼
70% of all BBH detections at both the O1 sensitivity andthe final design sensitivity of aLIGO.The star formation history model of Neijssel et al.(2019) was tuned to the gravitational-wave observations, andexplaining the relatively high masses of observed BBHs re-quired significant high-redshift, low-metallicity star forma-tion. The inclusion of CHE BBHs naturally yields a popu-lation of high-mass sources, allowing for the high-mass starformation rate to be reduced in line with the Madau & Dick-inson (2014); Madau & Fragos (2017) models. This wouldnaturally bring rate predictions in line with the O1 and O2observations and correspondingly reduce predicted detectionrates for future detectors.Using the preferred cosmic metallicity star formationmodel of Neijssel et al. (2019), as we do here, and assuming f WR = 1, we predict a total BBH detection rate of ≈ ≈
37 at O1 sen-sitivity), with ≈
470 ( ≈
27) of these coming from the CHEchannel. The CHE BBH detection rates are a factor ∼ ≈
250 ( ≈
13) CHE BBHs per year may bedetected at aLIGO design (O1) sensitivity. The differencesin the assumed metallicity-specific star formation rates inthese studies are responsible for much of this difference.We note that in both Figures 12 and 13 the order ofthe lines with respect to the number of detections does notmatch the order of the WR mass loss multipliers. This is dueto the interplay between the formation rate of BBHs andtheir delay times as a function of f WR , which are describedin Figures 8 and 9 and associated discussion. For example,in the absence of WR winds ( f WR = 0), reduced delay times Figure 14.
Chirp mass posteriors for the 10 BBH mergers de-tected during the first and second aLIGO observing runs (Abbottet al. 2019) are shown in colour at the bottom, with labels at top.These are randomly sampled to construct the cumulative densityfunctions shown in lavender (each curve corresponds to a cumu-lative distribution through 10 samples, one from each posterior).Cumulative density functions for
COMPAS chirp mass predic-tions based on f WR = 1 . due to a lack of binary widening relative to simulations withhigher WR mass loss rates mean that very few CHE BBHs,which predominantly form at lower metallicities and thushigher redshifts, merge in the local Universe, where theywould be detectable. The cumulative distribution functions for the modelled chirpmass distribution of detectable BBH mergers are shown inFigure 14. The dark blue lines indicate the chirp mass dis-tribution of all BBHs while the light blue lines indicate thechirp mass distribution of CHE BBHs. In both cases, resultsfor the WR mass loss multiplier f WR = 1 . COMPAS models – the number of BBHs detectedduring O1 and O2 – in order to indicate the variation due tosampling fluctuations. Some of the granularity in
COMPAS models is due to the discreteness of the metallicity grid (seeDominik et al. 2015; Neijssel et al. 2019, for a discussion)that can be addressed with improved sampling.As mentioned previously, CHE BBHs are more massivethan typical non-CHE BBHs. The initial masses of CHEBBHs must be high to allow for CHE (see Figure 3). More-over, CHE in our model allows stars to convert all of theirmass to helium, whereas non-CH massive stars typicallyhave (cid:38)
50% of their mass in hydrogen-rich envelopes, whichthey lose prior to collapse into black holes in the course of
MNRAS000
MNRAS000 , 1–15 (2020) J. Riley et al.
Figure 15.
The fraction of BBHs formed through the CHE chan-nel among all BBHs detectable at aLIGO O1 sensitivity, plottedas a function of chirp mass. binary evolution. This is highlighted in Figure 15, which in-dicates the fraction of all BBHs detectable at aLIGO O1sensitivity that formed through the CHE channel, as a func-tion of chirp mass. The CHE channel dominates the produc-tion of BBHs at high chirp masses, particularly for reducedWR mass loss models, when it yields increasingly large chirpmasses ( (cid:38) M (cid:12) in the absence of WR winds).Figure 14 allows for a direct comparison between themodelled chirp mass distribution and the aLIGO observa-tions from the first two observing runs. The individual pos-terior samples from the 10 aLIGO BBH detections duringthose observing runs are plotted at the bottom of the plot.Randomly sampled cumulative distribution functions of thechirp mass of observed events are constructed by taking 10random samples, one from each of the 10 aLIGO observationposteriors and displayed as light lavender curves. The over-lap of the lavender and dark blue lines in Figure 14 showsthat the COMPAS model of BBH formation, which includesthe contribution of CHE, yields a chirp mass distribution ofdetectable BBH mergers that is consistent with detectionsduring the first two aLIGO observing runs.
We described the model of chemically homogeneous evolu-tion (CHE) that we implemented in the rapid populationsynthesis code
COMPAS . We used MESA models to deter-mine the critical rotation thresholds for CHE, and providedfits that can be used in other rapid binary population synthe-sis applications. We synthesised 12 million binary systemsover a range of metallicities (30 metallicities evenly spacedacross the range − ≤ log Z ≤ − . f WR ∈ { . , . , . , . } ). We confirmedthat our simplified models match detailed binary evolutionsimulations (Marchant et al. 2016; du Buisson et al. 2020)well.We investigated the contribution of CHE and non-CHE channels to BBH formation under a single set of assump-tions. We found that the CHE channel may contribute morethan half, and perhaps as much as three quarters, of allaLIGO BBH detections arising from isolated binary evolu-tion. CHE BBHs may represent (cid:38)
80% of detectable sourceswith the highest chirp masses of (cid:38) M (cid:12) . A comparison be-tween our model population and the population of detectedbinaries from the first two advanced detector observing runsindicates that the current model over-predicts the total num-ber of sources by a factor of ∼
2, but matches the observedchirp mass distribution.We made a number of simplifying assumptions in thisstudy that can be investigated and improved on in the fu-ture. We generally erred on the side of being conservativeabout CHE predictions: • We used Hurley et al. (2000) MS models to set the radiiand mass loss rates of CH MS stars. At high metallicity,the Hurley et al. (2000) ZAMS radii are larger than theradii from MESA models which we used to determine criticalrotation rates required to keep the stars on CHE, so usingconsistent radii may expand the parameter space for CHE. • We used simplified tidal interaction assumptions underwhich CH stars are immediately tidally synchronised, yetdo not store angular momentum. Accounting for the angu-lar momentum stored in stars – and the additional angularmomentum carried away by winds from a rotating star –impacts the response of the binary’s orbit to mass loss, andreduces the amount of orbital widening by wind mass lossin close binaries. • Contrary to our simplified assumptions, winds may interactwith binary companions. This is particularly true in closebinaries, when the wind speeds are comparable to the orbitalspeeds, and wind interactions may produce additional dragand reduce the amount of orbital widening (e.g., Brookshaw& Tavani 1993; MacLeod & Loeb 2020). • We ignored the possibility of initially non-CH stars switch-ing to CHE in response to mass accretion. • We assumed that all mass loss in PPISNe happens instan-taneously, rather than over several pulsations (although thefirst pulsation is likely to be dominant, so this approximationmay not be especially problematic).Our predicted BBH merger rate at redshift zero of50 Gpc − yr − (including 20 Gpc − yr − from the CHEchannel) for the default WR mass loss rate f WR = 1 . MNRAS , 1–15 (2020) hemically Homogeneous Evolution: A rapid population synthesis approach CHE sources should be ∼ .
4, much lower than the super-critical spins expected at WR star formation. The fractionof stellar angular momentum lost in winds during the WRphase can be estimated as∆ LL ∼ (cid:18) R WR R WR , g (cid:19) ∆ MM , (5)where the ratio of the radius of the WR star to its gyrationradius is R WR /R WR , g ∼
10. Thus, Wolf-Rayet winds couldlose the overwhelming bulk of the angular momentum thatCHE stars have, as long as ∆
M/M > .
01, which is trueeven at Z = 0 . Z (cid:12) if WR mass loss is not suppressed (seeFigure 5).CHE BBH progenitors could yield interesting observa-tional candidates. Systems such as WR20a (Rauw et al.2004) and BAT99-32 (Shenar et al. 2019) may belong inthis category. The metallicity of the Galaxy is too high toallow for merging CHE BBHs according to our models, butwe expect them to be formed at a rate of ∼ × − peryear in the Magellanic clouds. Given the typical MS and WRphase lifetimes of 3 × and 3 × years, respectively, wemay hope to detect ∼
10 MS CH binaries and ∼ ACKNOWLEDGEMENTS
Simulations in this paper made use of the
COMPAS rapidpopulation synthesis code, which is freely available at http://github.com/TeamCOMPAS/COMPAS . The version of
COM-PAS used for these simulations was v02.11.01a, built specif-ically for these simulations; functionality in this release wasintegrated into the public
COMPAS code base in v02.11.04.The authors thank Selma de Mink and other colleaguesin Team COMPAS, as well as Morgan MacLeod, for help-ful discussions. We also thank Tim Riley for assistance inrunning COMPAS simulations. IM is a recipient of the Aus-tralian Research Council Future Fellowship FT190100574.AVG acknowledges funding support by the Danish NationalResearch Foundation (DNRF132).
DATA AVAILABILITY
The data underlying this article will be available via https://zenodo.org/communities/compas/ . REFERENCES
Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, Classicaland Quantum Gravity, 34, 044001, doi: —. 2019, ApJ, 882, L24, doi:
Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019, PhysicalReview X, 9, doi:
Abbott, B. P., et al. 2016, Phys. Rev. Lett., 116, 061102, doi:
Abbott, R., Abbott, T. D., Abraham, S., Acernese, F., et al. 2020,arXiv e-prints, arXiv:2009.01075. https://arxiv.org/abs/2009.01075
Abt, H. A. 1983, ARA&A, 21, 343, doi:
Almeida, L. A., Sana, H., Mink, S. E. d., et al. 2015, The As-trophysical Journal, 812, 102, doi:
Antonini, F., Chatterjee, S., Rodriguez, C. L., et al. 2016, TheAstrophysical Journal, 816, 65, doi:
Arcavi, I., Howell, D. A., Kasen, D., et al. 2017, Nature, 551, 210,doi:
Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009,ARA&A, 47, 481, doi:
Barkat, Z., Rakavy, G., & Sack, N. 1967, Physical Review Letters,18, 379, doi:
Bartos, I., Kocsis, B., Haiman, Z., & M´arka, S. 2017, ApJ, 835,165, doi:
Belczynski, K., Bulik, T., Fryer, C. L., et al. 2010, The Astrophys-ical Journal, 714, 1217–1226, doi:
Belczynski, K., Taam, R. E., Kalogera, V., Rasio, F. A., & Bulik,T. 2007, ApJ, 662, 504, doi:
Brookshaw, L., & Tavani, M. 1993, ApJ, 410, 719, doi:
Chruslinska, M., Nelemans, G., & Belczynski, K. 2019, MNRAS,482, 5012, doi: de Mink, S. E., Cantiello, M., Langer, N., et al. 2009, A&A, 497,243, doi: de Mink, S. E., & Mandel, I. 2016, Monthly Notices of the RoyalAstronomical Society, 460, 3545–3553, doi:
Dominik, M., Berti, E., O’Shaughnessy, R., et al. 2015, ApJ, 806,263, doi: du Buisson, L., Marchant, P., Podsiadlowski, P., et al. 2020,arXiv e-prints, arXiv:2002.11630. https://arxiv.org/abs/2002.11630
Eddington, A. S. 1929, Monthly Notices of the Royal Astronom-ical Society, 90, 54, doi:
Eldridge, J. J., Stanway, E. R., Xiao, L., et al. 2017, Publ. Astron.Soc. Australia, 34, e058, doi:
Farmer, R., Renzo, M., de Mink, S. E., Marchant, P., & Justham,S. 2019, ApJ, 887, 53, doi:
Fowler, W. A., & Hoyle, F. 1964, ApJS, 9, 201, doi:
Fraley, G. S. 1968, Ap&SS, 2, 96, doi:
Gal-Yam, A. 2012, Science, 337, 927, doi:
Heger, A., Langer, N., & Woosley, S. E. 2000, The AstrophysicalJournal, 528, 368–396, doi:
Humphreys, R. M., & Davidson, K. 1979, ApJ, 232, 409, doi:
Hurley, J. R., Pols, O. R., & Tout, C. A. 2000, Mon. Not. Roy. As-tron. Soc., 315, 543, doi:
Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, Mon. Not. Roy. As-tron. Soc., 329, 897, doi:
MNRAS000
MNRAS000 , 1–15 (2020) J. Riley et al.
Hut, P. 1981, A&A, 99, 126Kobulnicky, H. A., Kiminki, D. C., Lundquist, M. J., et al. 2014,ApJS, 213, 34, doi:
Koenigsberger, G., Morrell, N., Hillier, D. J., et al. 2014, TheAstronomical Journal, 148, 62, doi:
Kroupa, P. 2001, Mon. Not. Roy. Astron. Soc., 322, 231, doi:
Langer, N. 1992, A&A, 265, L17MacLeod, M., & Loeb, A. 2020, arXiv e-prints, arXiv:2007.07252. https://arxiv.org/abs/2007.07252
Madau, P., & Dickinson, M. 2014, Annual Review of As-tronomy and Astrophysics, 52, 415–486, doi:
Madau, P., & Fragos, T. 2017, ApJ, 840, 39, doi:
Maeder, A. 1987, A&A, 178, 159—. 1989, Astrophysics and Space Science Library, Vol. 157, Onthe Evolutionary Status and Instability Mechanism of the Lu-minous Blue Variables, ed. K. Davidson, A. F. J. Moffat, &H. J. G. L. M. Lamers, 15Maeder, A., & Meynet, G. 2000, Annual Review of Astronomyand Astrophysics, 38, 143–190, doi:
Mandel, I., & de Mink, S. E. 2016, Monthly Notices of the RoyalAstronomical Society, 458, 2634–2647, doi:
Mandel, I., & Farmer, A. 2018, ArXiv e-prints. https://arxiv.org/abs/1806.05820
Mapelli, M. 2018, arXiv e-prints. https://arxiv.org/abs/1809.09130
Marchant, P., Langer, N., Podsiadlowski, P., Tauris, T. M., &Moriya, T. J. 2016, Astron. Astrophys., 588, A50, doi:
Marchant, P., Renzo, M., Farmer, R., et al. 2019, The Astrophys-ical Journal, 882, 36, doi:
Miller, M. C., & Lauburg, V. M. 2009, The Astrophysical Journal,692, 7, doi:
Moe, M., & Di Stefano, R. 2017, ApJS, 230, 15, doi:
Neijssel, C. J., Vigna-G´omez, A., Stevenson, S., et al. 2019,arXiv e-prints, arXiv:1906.08136. https://arxiv.org/abs/1906.08136 ¨Opik, E. 1924, Publications of the Tartu Astrofizica Observatory,25, 1Pasquali, A., Langer, N., Schmutz, W., et al. 1997, ApJ, 478, 340,doi:
Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3,doi:
Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4,doi:
Paxton, B., Marchant, P., Schwab, J., et al. 2015, The Astro-physical Journal Supplement Series, 220, 15, doi:
Peters, P. C. 1964, Physical Review, 136, B1224Petrovic, J., Langer, N., Yoon, S. C., & Heger, A. 2005, A&A,435, 247, doi:
Pols, O. R., Schr¨oder, K.-P., Hurley, J. R., Tout, C. A., &Eggleton, P. P. 1998, MNRAS, 298, 525, doi:
Punturo, M., Abernathy, M., Acernese, F., et al. 2010, Classicaland Quantum Gravity, 27, 084007, doi:
Rauw, G., De Becker, M., Naz´e, Y., et al. 2004, A&A, 420, L9,doi:
Renzo, M., Farmer, R., Justham, S., et al. 2020, arXiv e-prints,arXiv:2002.05077. https://arxiv.org/abs/2002.05077
Rodriguez, C. L., Morscher, M., Pattabiraman, B., et al. 2015, Phys. Rev. Lett., 115, 051101, doi:
Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Science, 337,444, doi:
Shenar, T., Sablowski, D. P., Hainich, R., et al. 2019, A&A, 627,A151, doi:
Sigurdsson, S., & Hernquist, L. 1993, Nature, 364, 423, doi:
Song, H. F., Meynet, G., Maeder, A., Ekstr¨om, S., & Eggenberger,P. 2016, A&A, 585, A120, doi:
Stevenson, S., Sampson, M., Powell, J., et al. 2019, ApJ, 882, 121,doi:
Stevenson, S., Vigna-G´omez, A., Mandel, I., et al. 2017, NatureCommunications, 8, 14906, doi:
Stone, N. C., Metzger, B. D., & Haiman, Z. 2017, MNRAS, 464,946, doi:
Sweet, P. A. 1950, Monthly Notices of the Royal AstronomicalSociety, 110, 548, doi:
Taylor, P., & Kobayashi, C. 2015, MNRAS, 448, 1835, doi:
Tutukov, A., & Yungelson, L. 1973, Nauchnye Informatsii, 27, 70Tutukov, A. V., & Yungelson, L. R. 1993, MNRAS, 260, 675,doi: van den Heuvel, E. P. J. 1976, in IAU Symposium, Vol. 73, Struc-ture and Evolution of Close Binary Systems, ed. P. Eggleton,S. Mitton, & J. Whelan, 35Vigna-G´omez, A., MacLeod, M., Neijssel, C. J., et al. 2020,arXiv e-prints, arXiv:2001.09829. https://arxiv.org/abs/2001.09829
Vigna-G´omez, A., Neijssel, C. J., Stevenson, S., et al. 2018, MN-RAS, 481, 4009, doi:
Vink, J. S., & de Koter, A. 2005, Astronomy & Astrophysics, 442,587–596, doi:
Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, Astron-omy & Astrophysics, 369, 574–588, doi: von Zeipel, H. 1924, MNRAS, 84, 665, doi:
Woosley, S. E. 2017, ApJ, 836, 244, doi: —. 2019, arXiv e-prints, arXiv:1901.00215. https://arxiv.org/abs/1901.00215
Woosley, S. E., Heger, A., & Weaver, T. A. 2002, Reviews ofModern Physics, 74, 1015, doi:
Yoon, S.-C., Langer, N., & Norman, C. 2006, Astronomy & As-trophysics, 460, 199–208, doi:
Yoshida, T., Umeda, H., Maeda, K., & Ishii, T. 2016, MNRAS,457, 351, doi:
APPENDIX A: CHE THRESHOLDS
We evolved single stars over a range of masses, metallici-ties, and rotational frequencies with MESA in order to findthe boundary between CHE and regular non-CH stellar evo-lution . The stars were evolved at a constant angular fre-quency until the end of the main sequence without massloss. The star was considered to evolve chemically homoge-neously if the difference between the helium fraction acrossthe star did not exceed 0 . The complete set of MESA input files necessary to reproducethese simulations will be made available after acceptance of themanuscript. MNRAS , 1–15 (2020) hemically Homogeneous Evolution: A rapid population synthesis approach Figure A1.
Rotational frequency threshold for chemically homogeneous evolution as a function of mass and metallicity. Downwardand upward triangles represent the slowest-rotating CHE model and fastest-rotating non-CHE model at the given mass and metallicity,respectively. The curves indicate the fits of Eqs. (A1) and (A2). which the star remains non chemically homogeneous (up-ward triangles), and the minimum rotational frequency atwhich the star becomes chemically homogeneous (downwardtriangles), for a grid of masses ranging from 10 to 150 M (cid:12) and three metallicities, Z = 0 . , . , . COMPAS and shown in Figure A1: ω M,Z = ω M , Z . .
09 ln ( Z0 . ) +1 (A1) where ω M , Z . = (cid:40)(cid:80) i =0 a i M i M . rad s − , M ≤
100 M (cid:12) (cid:80) i =0 a i i M . rad s − , M > (cid:12) (A2) and a = 5.7914 × − a = − × − a = − × − a = 1.0150e × − a = − × − a = 2.9051 × − We expect these fits to be valid over the range wherethey are constructed (10 M (cid:12) ≤ M ≤ M (cid:12) , 10 − ≤ Z ≤ .
01) but caution should be exercised if the fits are extrap-olated significantly beyond these boundaries.
MNRAS000