COSMIC Variance in Binary Population Synthesis

The main output from COSMIC is the fixed population, a collection of binary systems that contains enough binaries to capture the underlying shape of the population's parameter distribution functions resulting from a user-specified SFH and binary evolution model. The fixed population can be convolved with more complex SFHs and scaled to a large number of astrophysical populations. These populations can then be synthetically "observed" and used to explore the variance in synthetic catalogs associated with the simulated binary population and binary evolution model.

The general process to simulate a fixed population is as follows:

• 1.  
  Select a binary evolution model and SFH.
• 2.  
  Generate an initial population based on user-selected models for the SFH and initial binary parameter distributions.
• 3.  
  Evolve the initial population according to the user-specified binary evolution model.
• 4.  
  If on the first iteration, compare a subset containing half of the simulated population to the total population to determine how closely the binary parameter distributions match one another (we introduce a quantitative "match" condition for making this assessment, described in Section 2.1). If on the second or later iteration, compare the population from the previous iteration to the population containing both the current and previous iterations.
• 5.  
  If the binary parameter distributions have converged, the population is called the "fixed population," which contains a large-enough sample to describe the essential statistical features of a binary evolution model.
• 6.  
  Scale the fixed population to astrophysical populations, weighted either by mass or by number, by sampling the fixed population with replacement.
• 7.  
  Apply a synthetic observation pipeline to generate a set of synthetic catalogs from the astrophysical populations.

Figure 1 illustrates the structure and process that COSMIC uses to generate the fixed population. The fixed population is simulated once for each binary evolution model; thus, a study with, for example, 10 binary evolution models will have 10 associated fixed populations. Astrophysical populations can then be sampled from the fixed population, convolved with an SFH, and used to generate a statistical set of synthetic catalogs for each model (i.e., a full study with 10 binary evolution models and one SFH may contain 10,000 synthetic catalogs).

Zoom In Zoom Out Reset image size

2.1.1. Initializing a Population

2.1.2. Convergence of Parameter Distributions

There is no formulaic way to a priori predict the required number of binaries to be evolved for a fixed population, because each population depends on a different binary evolution model. The ideal number of simulated systems in a fixed population is such that the population adequately describes the final parameter distribution functions while not simulating so many systems as to be inefficient. To quantify this number, we develop a discrete match criteria, based on the cross ambiguity function, or overlap function, used in matched filtering techniques (e.g., Equation (6) of Chatziioannou et al. 2017).

We use independently generated histograms for each binary parameter, with bin widths determined using Knuth's rule (Knuth 2006), implemented in astropy (Astropy Collaboration et al. 2013), to track the distribution of each parameter as successive populations are simulated and cumulatively added to the fixed population. Before each histogram is generated, we ensure that similar bin widths are used for each parameter by transforming each set of binary parameters to lie between 0 and 1. We then enforce the physical limits of the simulated systems (e.g., positive definite values for mass and orbital period and eccentricities between 0 and 1) by taking the inverse logistic transform, or logit, which redistributes the transformed data to have limits between −inf and inf.

We define the match as

$$\begin{eqnarray}&&\mathrm{match}=\displaystyle \frac{{\displaystyle \sum }_{k=1}^{N}{P}_{k,i}{P}_{k,i+1}}{\sqrt{{\displaystyle \sum }_{k=1}^{N}({P}_{k,i}{P}_{k,i}){\displaystyle \sum }_{k=1}^{N}({P}_{k,i+1}{P}_{k,i+1})}},\end{eqnarray} \tag{ 1 }$$

where P_k,i denotes the probability for the kth bin for the ith iteration. The match is limited to values between 0 and 1 and tends to unity as the parameter distributions converge to a distinct shape. The match is specified by the user and can be set to any value. For the study in Section 4, we set match ≥1–10⁻⁵ as a fiducial choice, but caution that the results of each simulated population should be carefully checked to confirm that the population does not contain artificial gaps due to low-number statistics in cases of very rare subpopulations.

Figure 2 uses the solar metallicity NS + NS population from Section 4 to illustrate how distributions of the semimajor axis converge to a distinct shape as more binaries are simulated. The evolution of the match (see Equation (1)) as a function of total simulated mass is also shown in terms of Log₁₀(1-match) to quantify the convergence of the distributions. As the total simulated mass increases, more NS + NS systems are added to the population, which fills in the normalized distribution and consequently drives the match toward one.

Zoom In Zoom Out Reset image size

Once the fixed population satisfies the user-specified convergence criteria, an astrophysical population can be sampled, and each binary in the astrophysical population can be assigned a position and orientation. The number of sources in each astrophysical population can be calculated by normalizing the size of the fixed population by the ratio of the mass of the astrophysical population to the mass of stars formed to produce the fixed population,

$$\begin{eqnarray}&&{N}_{\mathrm{astro}}={N}_{\mathrm{fixed}}\displaystyle \frac{{M}_{\mathrm{astro},\mathrm{tot}}}{{M}_{\mathrm{fixed},\mathrm{stars}}},\end{eqnarray} \tag{ 2 }$$

or by the ratio of the number of stars in the astrophysical population to the total number of stars formed to produce the fixed population,

$$\begin{eqnarray}&&{N}_{\mathrm{astro}}={N}_{\mathrm{fixed}}\displaystyle \frac{{N}_{\mathrm{astro},\mathrm{tot}}}{{N}_{\mathrm{fixed},\mathrm{stars}}}.\end{eqnarray} \tag{ 3 }$$

For each astrophysical realization, N_astro binaries are sampled with replacement from the fixed population. These sampled binaries can also be assigned a three-dimensional position and an inclination (i), an argument of periapsis (ω), and a longitude of the ascending node (Ω).

COSMIC allows for several general Galactic position distributions. For the axisymmetric thin and thick disks, the radial and vertical distributions are assumed to be independent. We adopt Galactic position distributions and Galactic component masses from McMillan (2011) for the thin disk, thick disk, and bulge as a fiducial model. However, COSMIC can distribute binaries according to exponential distributions in the radial and vertical directions with any scale height, as well as in spherically symmetric distributions.

Following McMillan (2011), we assume the mass in the thin and thick disks to be 4.32 × 10¹⁰ M_⊙ and 1.44 × 10¹⁰ M_⊙, respectively, while we assume the mass of the bulge to be 8.9 × 10⁹ M_⊙. Finally, we emphasize that while COSMIC has utility functions implemented to distribute binaries in the Milky Way, the spatial distributions are independent of the binary evolution prescriptions. Thus, multiple spatial distributions can be assumed for a single fixed population, and more sophisticated spatial distribution models can be integrated into future releases of COSMIC.

Mass loss through stellar winds plays an important role in compact object formation because this determines the mass of the star just before it undergoes a supernova (SN) explosion (e.g., Mapelli et al. 2009; Belczynski et al. 2010; Fryer et al. 2012; Giacobbo et al. 2018). In recent years, models of stellar winds have been revised to reflect updates in our understanding of various relevant physical processes. We have amended the original BSE stellar wind prescriptions with several of these prescriptions, summarized below.

Recent work has shown that line-driven winds exhibit a strong dependence upon metallicity, both on the main sequence (MS) and post-MS (e.g., Vink et al. 2001, 2011; Meynet & Maeder 2005; Gräfener & Hamann 2008). We have updated the original BSE prescription to include metallicity-dependent winds for O and B stars as well as for Wolf–Rayet stars. The winds for O and B stars are treated according to the prescription of Vink et al. (2001), which considers stars with temperatures 12,500 K < T_eff < 22,500 K and 27,500 K < T_eff < 50,000 K separately. As in Rodriguez et al. (2016), we adopt the methods of Dominik et al. (2013), where the high- and low-temperature prescriptions are extended to T_eff = 25,000 K. Wolf–Rayet star winds are treated according to Vink & de Koter (2005).

Additionally, recent models suggest that as stars approach the limit imposed by the electron-scattering Eddington factor, ${{\rm{\Gamma }}}_{e}={\chi }_{e}L/(4\pi {cGM})$ where χ_e is electron scattering opacity, winds may become insensitive to metallicity (e.g., Gräfener & Hamann 2008; Vink et al. 2011; Chen et al. 2015). Thus, we include Eddington-limited winds using the prescriptions described in Gräfener & Hamann (2008) and Giacobbo et al. (2018).

If the primary star loses mass through a stellar wind, the secondary may accrete some of the ejected material as it orbits through it. The accretion rate onto the secondary can be estimated according to a Bondi–Hoyle type mechanism (Bondi & Hoyle 1944), which is sensitive to the velocity of the wind lost from the primary, v_W² = 2β_WGM/R, where β_W is a constant depending on stellar type (see Equation (9) of Hurley et al. 2002). We have implemented β_W values from Belczynski et al. (2008) and choose this as a default. In this case, β_W is a function of stellar type. For H-rich MS stars, β_W = 0.125 for masses below 1.4 M_⊙, β_W = 7.0 for masses above 120 M_⊙, and it is linearly interpolated between. For H-rich giants, β_W = 0.125 regardless of mass. For He-rich stars, including MS and giants, β_W = 0.125 for masses below 10 M_⊙, β_W = 7.0 for masses above 120 M_⊙, and it is linearly interpolated between.

The stability of Roche-lobe overflow mass transfer is determined through the same process as in Hurley et al. (2002), using critical mass ratios determined from radius-mass exponents (Webbink 1985) where q_crit = m_donor/m_accretor. COSMIC allows both models for critical mass ratios for giant star donors detailed in Section 2.6.1 of Hurley et al. (2002). We have added new models that are consistent with binary_c, following Claeys et al. (2014), where we have modified the critical mass ratios of MS/helium-MS stars and binaries with degenerate accretors based on de Mink et al. (2007). For the models that are consistent with Claeys et al. (2014), we adopt the Roche-lobe overflow mass transfer rates in their Equation (10) in the case of stable mass transfer. We have also added a model that is consistent with StarTrack following Belczynski et al. (2008). The most notable difference in the critical mass ratios between StarTrack and BSE or binary_c is the treatment of helium-MS donors where q_crit = 1.7 and for helium stars on the Hertzsprung gap or giant branch where q_crit = 3.5 based on the findings of Ivanova et al. (2003). Finally, we note that COSMIC allows for user-specified critical mass ratios based on evolutionary stages to allow easy integration with future models as they arise. As a fiducial model for critical mass ratios, we choose the standard model reported in Hurley et al. (2002), which does not differentiate between models for degenerate and nondegenerate accretors.

We employ the standard αλ model for common envelope (CE) evolution as done in Hurley et al. (2002). In this case, systems that undergo unstable mass transfer enter into a CE, which can be expelled by the injection of orbital energy from the binary. In this formalism, λ is a factor that determines the binding energy of the envelope to its stellar core, while α is the efficiency factor for injecting orbital energy into the envelope. As with the currently available version of BSE, COSMIC defaults to a variable λ that depends on the evolutionary state of the star following the description in the appendix of Claeys et al. (2014). However, constant λ is also allowed as an option.

As a fiducial CE model, we use the variable binding energy parameter in conjunction with a constant CE efficiency parameter α = 1.0 following previous results (e.g., Nelemans et al. 2001a; Dominik et al. 2012). However, previous studies of post-CE binaries point to an efficiency as low as α = 0.2 (Zorotovic et al. 2010; Toonen & Nelemans 2013; Camacho et al. 2014). On the other hand, detailed modeling of the CE phase for double NS (DNS) progenitors suggests CE efficiencies may be as high as α ≈ 5 (Fragos et al. 2019), which may also reduce the tension between rate predictions from CE channels (e.g., Mapelli & Giacobbo 2018) and the empirical DNS merger rate derived from LIGO/Virgo (Abbott et al. 2017).

We also include an option for a "pessimistic CE" scenario, in which unstable mass transfer from a donor star without a well-developed core–envelope structure is always assumed to lead to a merger (Belczynski et al. 2008). This assumption applies to hydrogen-rich and helium stars on the main sequence and Hertzsprung gap, as well as all white dwarfs. When this option is set, we also assume mergers to occur when unstable mass transfer is triggered from donor stars without a clear entropy jump at the core–envelope boundary (Ivanova & Taam 2004), which includes stars on the Hertzsprung gap. Note, however, that the outcomes of CEs in this case are still uncertain (e.g., Deloye & Taam 2010).

Mass transfer involving an evolved He-star donor and a compact object is a critical phase of binary evolution for forming hardened compact binaries that can merge within a Hubble time as GW sources. For progenitors of DNSs, detailed modeling of this phase finds that mass transfer typically proceeds stably (Ivanova et al. 2003; Tauris et al. 2013, 2015) and does not lead to the onset of a CE. Stable mass transfer is also corroborated by studies that compare population modeling to the properties of Galactic DNSs (Vigna-Gómez et al. 2018). Post-mass-transfer separations are typically wider if this phase is modeled stably versus unstably. We allow for the stable mass-transfer evolution of He-star donors with compact object companions to be approximated using the fitting formulae in Tauris et al. (2015), which can fit for either the post-mass-transfer separation and remaining He-star envelope mass, or just the post-mass-transfer separation.

One important change implemented in COSMIC is the treatment of supernovae that directly follow a CE phase. By default, if a supernova immediately follows a CE (which typically occurs for evolved He-star donors), BSE determines the post-SN barycentric velocity, orbital properties, and survival using the post-CE orbital separation and pre-CE stellar mass. Therefore, the mass loss in the supernova includes the mass of the donor-star envelope that formed the CE—a self-inconsistent treatment because the ejection of this envelope is used to determine the hardening during the CE phase. This inconsistency leads to artificially amplified mass-loss kicks during the supernova, which are particularly apparent in extremely tight binaries. Therefore, by default, we use the post-CE separation and post-CE mass (which does not include the mass of the envelope) when determining the effect of the subsequent supernova. The effect of this change, particularly on DNS population properties, is explained in more detail in Zevin et al. (2019).

Finally, we note that contact systems can occur if both stellar components overflow their Roche radii or if a stellar component's radius is larger than the binary's periastron distance. In this case, a common envelope is always triggered, regardless of the mass transfer stability criteria. This simplification does not allow for long-lived contact binaries, such as W UMa, to be studied in detail, but their formation rates and progenitor populations can be investigated.

3.3.1. Standard Core-collapse Supernovae

3.3.2. Electron-capture Supernovae

3.3.3. Ultrastripped Supernovae

3.3.4. Pair-instability and Pulsational Pair-instability Supernovae

We have added new prescriptions for the spins of newly formed BHs from collapsing massive stars. Due to the uncertain values associated with BH natal spins, we allow the spins of all BHs to be set to a specific Kerr value (specified by the user) or drawn from a uniform distribution whose bounds are also user specified. In addition, we have included the prescriptions for BH spin based on the precollapse CO core mass of the progenitor from Belczynski et al. (2020), which result in high spins for low-mass BHs (≲30M_⊙, depending on metallicity) and low spins for high-mass BHs.

The latter prescription, based on stellar models computed by the Geneva stellar evolution code (Eggenberger et al. 2008), assumes that angular momentum is transported in massive stars via meridional currents and is not efficient enough to spin down the core prior to collapse. This is in contrast to newer work (e.g., Fuller & Ma 2019) suggesting that the Taylor–Spruit magnetic dynamo may allow for extremely efficient angular momentum transport through the envelopes of massive stars, producing BHs with dimensionless Kerr spin parameters a/M ∼ 0.01. We do not include the latter prescription in COSMIC, nor do we allow for the BH spin to be increased by accretion or mergers. These effects will be added in a future release.

We have updated BSE to implement the NS magnetic field and spin period evolution following Kiel et al. (2008) and Ye et al. (2019). All NSs are born with magnetic fields and spin periods that match the observed young pulsars (Manchester et al. 2005); their initial magnetic fields and spin periods are randomly drawn in the range of 10^11.5–10^13.8 G and 30–1000 ms, respectively.

For single NSs or NSs in detached binaries, we assume magnetic dipole radiation for the spin period evolution. The spin-down rate is calculated by

$$\begin{eqnarray}&&\dot{P}=K\displaystyle \frac{{B}^{2}}{P},\end{eqnarray} \tag{ 5 }$$

where P is the spin period, B is the surface magnetic field, and $K=9.87\times {10}^{-48}\,\mathrm{yr}\,{{\rm{G}}}^{-2}$. We also assume that the magnetic fields follow exponential decays in a timescale of τ = 3 Gyr (Kiel et al. 2008),

$$\begin{eqnarray}&&B={B}_{0}\exp \left(\displaystyle \frac{-T}{\tau }\right),\end{eqnarray} \tag{ 6 }$$

where B₀ is the initial magnetic field, and T is the age of the NSs.

On the other hand, binary evolution will affect the NS magnetic fields and spin periods. For NSs in nondetached binaries, their magnetic fields can change significantly during mass accretion on short timescales. We assume "magnetic field burying" (e.g., Bhattacharya & van den Heuvel 1991; Rappaport et al. 1995; Kiel et al. 2008; Tauris et al. 2012) for the magnetic field decay during mass accretion:

$$\begin{eqnarray}&&B=\displaystyle \frac{{B}_{0}}{1+({\rm{\Delta }}M/{10}^{-6}{M}_{\odot })}\exp \left(-\displaystyle \frac{T-{t}_{\mathrm{acc}}}{\tau }\right).\end{eqnarray} \tag{ 7 }$$

Here, ΔM is the mass accreted, and t_acc is the accretion duration. These NSs are spun up according to the amount of angular momentum transferred (Hurley et al. 2002, Equation (54)).

Furthermore, when a neutron star merges with another star (e.g., main-sequence star, giant, or WD) and the final product is an NS, the magnetic field and spin period of this NS are reset by drawing a new magnetic field and spin period from the same initial ranges. If a millisecond pulsar (MSP) is involved in the collision or merger, however, the newly formed NSs are assigned different ranges of initial magnetic fields and spin periods to match those of MSPs, and the newborn NSs will remain MSPs. In this case, the magnetic fields and spin periods are randomly drawn from ranges 10⁸–10^8.8 G and 3–20 ms, respectively. In addition, we assume a lower limit of 5 × 10⁷ G for the NS magnetic field (Kiel et al. 2008). No lower limit is assumed for the spin periods.

When two stellar cores spiral in toward one another, the outcome of the subsequent merger depends upon the internal structures (i.e., density profiles) of the two objects. If one of the cores is much denser than the other (for example, if a BH/NS merges with a giant star), a common-envelope-like event ensues. However, if the two stars have comparable compactness (for example, in the case of a roughly equal-mass MS–MS merger), then the merger may result in efficient mixing of the two objects. In either case, the stellar age of the merger product must be specified. Here we adopt the prescriptions of Hurley et al. (2002) to determine outcomes of stellar mergers, with one exception, outlined below.

In the case of an MS–MS merger, we assume the stellar age of the MS merger product is given by

$$\begin{eqnarray}&&{t}_{3}={f}_{\mathrm{rejuv}}\displaystyle \frac{{t}_{\mathrm{MS}3}}{{M}_{3}}\left(\displaystyle \frac{{M}_{1}{t}_{1}}{{t}_{\mathrm{MS}1}}+\displaystyle \frac{{M}_{2}{t}_{2}}{{t}_{\mathrm{MS}2}}\right)\end{eqnarray} \tag{ 8 }$$

where M₁ and M₂ are the masses of the two merger components; M₃ = M₁ + M₂ is the mass of the merger product (assuming the stars merge without mass loss; Hurley et al. 2002); t_MS1, t_MS2, and t_MS3 are the MS lifetimes of the two merger components and the merger product, respectively; and t₁ and t₂ are the stellar ages of the two merger components at the time of merger. Also, f_rejuv is a factor that determines the amount of rejuvenation the merger product experiences through mixing. This factor of course depends upon the internal structure of the two stars, as well as the nature of the merger (i.e., the relative velocity of the two objects at coalescence). In an original BSE, a fixed value of 0.1 is assumed for f_rejuv. However, in many instances, this likely leads to overrejuvenation of the merger product. Here we include f_rejuv as a free parameter and adopt f_rejuv = 1 as our default value.

The outcome of stellar mergers and collisions is expected to play a critical role in dense star clusters where dynamical interactions lead to a pronounced increase in stellar mergers or collisions relative to isolated binaries (e.g., Hills & Day 1976; Bacon et al. 1996; Lombardi et al. 2002; Fregeau & Rasio 2007; Leigh et al. 2011). The details of these merger products have important implications for blue straggler stars (e.g., Sandage 1953; Chatterjee et al. 2013), as well as the formation of massive BHs (e.g., Portegies Zwart & McMillan 2002; Gürkan et al. 2004; Kremer et al. 2020). Thus when modeling stellar and binary evolution in collisional environments like globular clusters, care must be taken when assigning the ages of stars upon collision or merger.

We adopt the merger products from the collision matrix of Hurley et al. (2002, Table 2), but we caution that the unmodified version of BSE contains a typo that causes the merger of two He-MS stars to produce an MS star. The collision matrix covers all merger types, from MS to compact remnants, and depends on the relative compactness of each component star's core. In the case of WD mergers with companion stars, the WD can mix completely with the stellar core; for example, a merger between a He WD and a star on the first giant branch will produce a giant star with a more massive helium core. If the cores have different compactness, for example, in a merger between a CO WD and a star on the first giant branch, the product will distribute the helium core of the giant around the CO core of the white dwarf, producing an early AGB star. For a detailed discussion of stellar merger products, see Section 2.7.2 of Hurley et al. (2002).

The characteristic strain of a GW source, as well as the signal-to-noise ratio (S/N) for a given GW detector, can be calculated in several different approximations, depending upon the properties of the source of interest. Specifically, the most general case of a sky- and polarization-averaged eccentric and chirping source can be approximated if the source is circular or stationary (nonchirping). In this section, we describe the computation of the characteristic strain and LISA S/N in four different regimes: eccentric and chirping, circular and chirping, eccentric and stationary, and circular and stationary.

In the most general case of an eccentric chirping source, the characteristic strain at the nth harmonic can be written as (e.g., Barack & Cutler 2004)

$$\begin{eqnarray}&&{h}_{c,n}^{2}=\displaystyle \frac{1}{{\left(\pi {D}_{L}\right)}^{2}}\left(\displaystyle \frac{2G}{{c}^{3}}\displaystyle \frac{{\dot{E}}_{n}}{{\dot{f}}_{n}}\right).\end{eqnarray} \tag{ 9 }$$

Here, D_L is the luminosity distance to the source, and f_n is the source-frame GW frequency of the nth harmonic given by f_n = n f_orb, where f_orb is the source-frame orbital frequency; f_n is related to the observed (detector frame) GW frequency, f_n,d, by ${f}_{n}={f}_{n,d}(1+z)$.

Here, ${\dot{E}}_{n}$ is the time derivative of the energy radiated in GWs at source-frame frequency f_n, which to lowest order is given by (e.g., Peters & Mathews 1963)

$$\begin{eqnarray}&&{\dot{E}}_{n}=\displaystyle \frac{32}{5}\displaystyle \frac{{G}^{7/3}}{{c}^{5}}{\left(2\pi {f}_{\mathrm{orb}}{{ \mathcal M }}_{c}\right)}^{10/3}\,g(n,e)\end{eqnarray} \tag{ 10 }$$

where ${{ \mathcal M }}_{c}$ is the source-frame chirp mass, which is related to detector-frame chirp mass, ${{ \mathcal M }}_{c,d}$, by

$$\begin{eqnarray}&&{{ \mathcal M }}_{c}=\displaystyle \frac{{{ \mathcal M }}_{c,d}}{1+z}=\displaystyle \frac{{\left({M}_{1}{M}_{2}\right)}^{3/5}}{{\left({M}_{1}+{M}_{2}\right)}^{1/5}}\displaystyle \frac{1}{1+z}.\end{eqnarray} \tag{ 11 }$$

Then ${\dot{f}}_{n}=n\,{\dot{f}}_{\mathrm{orb}}$ is given by

$$\begin{eqnarray}&&{\dot{f}}_{n}=n\displaystyle \frac{96}{10\pi }\displaystyle \frac{{\left(G{{ \mathcal M }}_{c}\right)}^{5/3}}{{c}^{5}}{\left(2\pi {f}_{\mathrm{orb}}\right)}^{11/3}\,F(e),\end{eqnarray} \tag{ 12 }$$

where $F{(e)=[1+(73/24){e}^{2}+(37/96){e}^{4}]/(1-{e}^{2})}^{7/2}$. Combining Equations (9), (10), and (12), we obtain the characteristic strain in the detector frame:

$$\begin{eqnarray}&&{h}_{c,n,d}^{2}=\displaystyle \frac{2}{3{\pi }^{4/3}}\displaystyle \frac{{\left(G{{ \mathcal M }}_{c}\right)}^{5/3}}{{c}^{3}{D}_{L}^{2}}\displaystyle \frac{1}{{f}_{n,d}^{1/3}{\left(1+z\right)}^{2}}{\left(\displaystyle \frac{2}{n}\right)}^{2/3}\displaystyle \frac{g\left(n,e\right)}{F\left(e\right)}.\end{eqnarray} \tag{ 13 }$$

For z ≈ 0, the distinction between the detector-frame and source-frame quantities becomes negligible. In this case, ${h}_{c,n,d}\approx {h}_{c,n}$, f_n,d ≈ f_n, and ${{ \mathcal M }}_{c,d}\approx {{ \mathcal M }}_{c}$, allowing us to write Equation (13) as

$$\begin{eqnarray}&&{h}_{c,n}^{2}=\displaystyle \frac{2}{3{\pi }^{4/3}}\displaystyle \frac{{\left(G{{ \mathcal M }}_{c}\right)}^{5/3}}{{c}^{3}{D}_{L}^{2}}\displaystyle \frac{1}{{f}_{n}^{1/3}}{\left(\displaystyle \frac{2}{n}\right)}^{2/3}\displaystyle \frac{g(n,e)}{F(e)}.\end{eqnarray} \tag{ 14 }$$

Henceforth, we ignore dependence on redshift for simplicity, an appropriate approximation because z ≪ 1 for the Galactic sources of interest in this analysis. For binaries with e = 0, all GW power is emitted in the n = 2 harmonic. Thus, the characteristic strain for circular and chirping binaries is

$$\begin{eqnarray}&&{h}_{c,2}^{2}=\displaystyle \frac{2}{3{\pi }^{4/3}}\displaystyle \frac{{\left(G{{ \mathcal M }}_{c}\right)}^{5/3}}{{c}^{3}{D}_{L}^{2}}\displaystyle \frac{1}{{f}_{2}^{1/3}},\end{eqnarray} \tag{ 15 }$$

where we have simply taken n = 2 in Equation (14).

For eccentric and stationary sources where $\dot{f}\lt f/{T}_{\mathrm{obs}}$, the dimensionless GW strain of the nth harmonic is given by

$$\begin{eqnarray}&&{h}_{n}^{2}=\displaystyle \frac{128}{5}\displaystyle \frac{{\left(G{{ \mathcal M }}_{c}\right)}^{10/3}}{{c}^{8}}\displaystyle \frac{{\left(\pi {f}_{\mathrm{orb}}\right)}^{4/3}}{{D}_{L}^{2}}\displaystyle \frac{g(n,e)}{{n}^{2}}.\end{eqnarray} \tag{ 16 }$$

Here we have divided Equation (15) by two times the number of cycles, N, observed for a source within a given frequency bin: $N={f}_{n}^{2}/\dot{{f}_{n}}$ (see, e.g., Moore et al. 2015). The strain for a stationary and circular source is obtained from Equation (16) and by adopting n = 2. The amplitude spectral density (ASD) for a stationary source is

$$\begin{eqnarray}&&{\mathrm{ASD}}_{n}={h}_{n}\sqrt{{T}_{\mathrm{obs}}},\end{eqnarray} \tag{ 17 }$$

where the power spectral density (PSD) is simply

$$\begin{eqnarray}&&{\mathrm{PSD}}_{n}={\mathrm{ASD}}_{n}^{2}.\end{eqnarray} \tag{ 18 }$$

Having written the expressions for characteristic strain in the four different circular/eccentric and stationary/chirping regimes, we now show how to calculate the LISA S/N. For chirping and eccentric sources, S/N is computed following Smith & Caldwell (2019) as

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{\rm{S}}}{{\rm{N}}}\right)}^{2}=\displaystyle \sum _{n=1}^{\infty }{\int }_{{f}_{\mathrm{start}}}^{{f}_{\mathrm{end}}}{\left[\displaystyle \frac{{h}_{c,n}({f}_{n})}{4{h}_{f}({f}_{n})}\right]}^{2}d\mathrm{ln}{f}_{n},\end{eqnarray} \tag{ 19 }$$

where ${h}_{c,n}$ is given by Equation (9), and h_f is the characteristic LISA noise curve, which we take from Robson et al. (2019). Here, f_start = nf_orb is the GW frequency emitted at the nth harmonic at the start of the LISA observation, and f_end is either the GW frequency at merger or the GW frequency of the nth harmonic of the orbital frequency of the binary at the end of the LISA observation time. The characteristic noise can be expressed as

$$\begin{eqnarray}&&{h}_{f}({f}_{n})=\sqrt{\displaystyle \frac{{f}_{n}\,{P}_{n}({f}_{n})}{{ \mathcal R }({f}_{n})}},\end{eqnarray} \tag{ 20 }$$

where ${P}_{n}({f}_{n})$ is the noise power spectral density of the detector, and ${ \mathcal R }({f}_{n})$ is the frequency-dependent signal response function.

For a chirping circular source, h_c,n in Equation (19) is replaced by h_c,2 (Equation (15)) so that the S/N is given by

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{\rm{S}}}{{\rm{N}}}\right)}^{2}=\displaystyle \sum _{n=1}^{\infty }{\int }_{{f}_{\mathrm{start}}}^{{f}_{\mathrm{end}}}{\left[\displaystyle \frac{{h}_{c,2}({f}_{n})}{4{h}_{f}({f}_{n})}\right]}^{2}d\mathrm{ln}{f}_{n}.\end{eqnarray} \tag{ 21 }$$

For stationary sources, it is no longer necessary to integrate over frequency space. In the stationary and eccentric regime, we obtain

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{\rm{S}}}{{\rm{N}}}\right)}^{2}\approx \displaystyle \sum _{n=1}^{\infty }{\left[\displaystyle \frac{{h}_{n}({f}_{n})}{4{h}_{f}({f}_{n})}\right]}^{2}{f}_{n}{T}_{\mathrm{obs}}.\end{eqnarray} \tag{ 22 }$$

Finally, for stationary and circular sources, we have

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{\rm{S}}}{{\rm{N}}}\right)}^{2}\approx {\left[\displaystyle \frac{{h}_{2}(f)}{4{h}_{f}(f)}\right]}^{2}{{fT}}_{\mathrm{obs}}.\end{eqnarray} \tag{ 23 }$$

In the equations above, T_obs is the LISA observation time, which we take to be 4 yr. In practice, the vast majority of sources are stationary or can be approximated as stationary because their evolution over the observation time spans less than 500 LISA bins. Based on a 4 yr observation time, this amounts to a frequency change of Δf_GW ≲ 5 × 10⁻⁶ Hz and thus has a negligible effect on S/N. In the following section, we assume all sources are stationary.

We convolve the fixed population with the SFHs and spatial distributions described in Section 2 to produce a population of WD, NS, and BH binaries born in the Galaxy. Based on their birth time and the age of each Galactic component, we evolve each binary up to the present and remove all systems that either merge or fill their Roche lobes. Although accreting systems may also be resolvable GW sources, and indeed may present rich opportunities for studying various aspects of binary evolution (e.g., Kremer et al. 2017; Breivik et al. 2018), we limit this study to only detached sources for simplicity. We note that since one model is used for this example population, the uncertainty across binary evolution models is not well represented. Future studies will investigate, in detail, predicted close binary populations and the uncertainty across several models.

Table 2 shows a summary of the number of sources per population and Milky Way component, as well as the number of systems for which S/N > 7 from a Milky Way population realization. Our approach for selecting resolved sources is more simple than the approach used in Korol et al. (2017) and Lamberts et al. (2019). Instead, we apply a running median with a window of 100 frequency bins to the total root power spectral density shown in Figure 4 following Benacquista & Holley-Bockelmann (2006). This produces a synthetic foreground signal similar to that shown in Korol et al. (2017). This foreground is added to the LISA power spectral density curve of Robson et al. (2019) and used to compute the S/N as described in Section 4.1.

Zoom In Zoom Out Reset image size

The ASD of each population in our Milky Way realization, along with the ASD of the LISA noise floor, is plotted in Figure 4. The vast majority of the signal in the ASD is due to the population of DWDs. The characteristic shape, especially the falloff near 1 mHz, is due to the number of systems radiating GWs in a given frequency bin whose width is determined by LISA's observing duration as ∼T_obs⁻¹. At low frequencies, the number of sources per bin is high, both because most sources preferentially have longer orbital periods and because the number of bins is low when compared with the number of bins at higher frequencies. As a result, at frequencies in excess of 1 mHz, there is a higher likelihood that a source occupies a unique frequency bin. The features present in the ASD are often smoothed out to produce a foreground (e.g., Littenberg et al. 2013; Korol et al. 2017; Lamberts et al. 2019). Here we present the unsmoothed ASD in order to illustrate the contribution of each population to the signal. The GW signals at harmonics of the orbital frequency are seen in the populations containing an NS or BH due to eccentric sources.

We note that different binary evolution models will produce different ASDs and thus different irreducible foregrounds. This may provide an avenue for distinguishing different binary evolution models. However, we caution that in order to properly determine how well LISA can constrain binary evolution models, resolved sources should be compared to ASDs produced from the same fixed populations. This ensures that self-consistent comparisons are made between different models.

Figure 5 shows the systems resolved with S/N > 7 above the irreducible foreground created from the running median of the ASD of the full population. The total number in each Galactic component is listed in Table 2. We find that the vast majority of the population of resolved compact binaries is composed of DWD systems. Our findings are broadly consistent with most previous work that uses population synthesis to predict LISA populations (e.g., Hils et al. 1990; Nelemans et al. 2001b; Ruiter et al. 2010; Yu & Jeffery 2010; Nissanke et al. 2012; Korol et al. 2017; Lamberts et al. 2018, 2019; Lau et al. 2020). We find that the number of resolved DWDs from our simulations is in agreement within a factor of ∼2 (Nelemans et al. 2001a; Ruiter et al. 2010; Nissanke et al. 2012; Korol et al. 2017; Lamberts et al. 2019). However, we note that we find relatively more He + He than CO + He DWDs than Lamberts et al. (2019). Our resolved DWD population has a factor of ∼4 smaller DWDs than Yu & Jeffery (2010). We find an order of magnitude more NS + WD binaries than Nelemans et al. (2001a), resulting from updated physics for NS progenitors as described in Section 3.

Zoom In Zoom Out Reset image size

We find a lower rate of NS + NS binaries detected by LISA than both Lau et al. (2020) and Andrews et al. (2020) by a factor of ∼3. The assumptions for both of these studies are different from the model used here. In particular, Lau et al. (2020), following Stevenson et al. (2017) and Vigna-Gómez et al. (2018), assume that the initial distribution of separation is limited to 1000 au, while we place an upper limit of 5000 au. They also use the Fryer et al. (2012) delayed model, instead of the rapid model for compact object formation. Finally, they assume that Roche-lobe overflow from stripped stars onto NSs is stable, which produces more NS + NS binaries at relatively short orbital periods than the model used in this study. The COMPAS_α, detailed in Stevenson et al. (2017), is the model most similar to the model chosen for this study and results in ∼10 resolved LISA sources, which is in good agreement with our findings. Andrews et al. (2020) links the NS + NS population resolved by LISA directly to the observed NS + NS merger rates derived from radio observations or LIGO/Virgo. The relative agreement of Andrews et al. (2020) and Lau et al. (2019) suggests that the model chosen for this study does not reproduce the Milky Way population of NS + NS binaries. Since this study is simply a proof of concept for the capabilities of COSMIC, we leave a careful analysis of the effects of different compact object models as a topic of future study.

Similar to Lamberts et al. (2018), we find relatively few resolved BH + BH binaries even though there are >10⁶ of them in our total simulated population. This is also true for BH + NS and NS + NS binaries, which are produced at lower rates. The high number of BH binaries relative to NS + NS binaries is due to the relationship between the mass and natal kick of the BHs and NSs at birth. In the Fryer et al. (2012) rapid model, BHs can be formed with a significant amount of fallback and thus reduced natal kicks, which are less likely to unbind the binary. However, BH binaries are mostly unresolved because they tend to occupy frequencies well below 1 mHz. We further note that we use the optimistic assumption that stars that undergo a common envelope on the Hertzsprung gap survive, leading to a larger population.