#  Investigating the Lower Mass Gap with Low Mass X-ray Binary Population
Synthesis

Jared C. Siegel Department of Astronomy and Astrophysics, University of
Chicago, Chicago, IL 60637, USA  Center for Interdisciplinary Exploration and
Research in Astrophysics (CIERA), Northwestern University, 1800 Sherman Ave,
Evanston, IL 60201, USA  Department of Astrophysical Sciences, Princeton
University, 4 Ivy Lane, Princeton, NJ 08544, USA  Ilia Kiato  Department of
Physics, University of Illinois at Urbana-Champaign, 1110 W Green St, Urbana,
IL 61801, USA  Center for Interdisciplinary Exploration and Research in
Astrophysics (CIERA), Northwestern University, 1800 Sherman Ave, Evanston, IL
60201, USA  Vicky Kalogera Department of Physics and Astronomy, Northwestern
University, 2145 Sheridan Ave, Evanston, IL 60208, USA  Center for
Interdisciplinary Exploration and Research in Astrophysics (CIERA),
Northwestern University, 1800 Sherman Ave, Evanston, IL 60201, USA
Christopher P. L. Berry SUPA, School of Physics and Astronomy, University of
Glasgow, Kelvin Building, University Ave, Glasgow G12 8QQ, UK  Center for
Interdisciplinary Exploration and Research in Astrophysics (CIERA),
Northwestern University, 1800 Sherman Ave, Evanston, IL 60201, USA  Thomas J.
Maccarone  Department of Physics & Astronomy, Texas Tech University, Box
41051, Lubbock, TX 79409, USA  Katelyn Breivik Center for Computational
Astrophysics, Flatiron Institute, 162 Fifth Ave, New York, NY 10010, USA  Jeff
J. Andrews  Center for Interdisciplinary Exploration and Research in
Astrophysics (CIERA), Northwestern University, 1800 Sherman Ave, Evanston, IL
60201, USA  Department of Physics, University of Florida, 2001 Museum Rd.,
Gainesville, FL 32611  Simone S. Bavera Departement d’Astronomie, Université
de Genève, Chemin Pegasi 51, CH-1290 Versoix, Switzerland  Aaron Dotter
Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA),
Northwestern University, 1800 Sherman Ave, Evanston, IL 60201, USA  Tassos
Fragos Departement d’Astronomie, Université de Genève, Chemin Pegasi 51,
CH-1290 Versoix, Switzerland  Konstantinos Kovlakas Departement d’Astronomie,
Université de Genève, Chemin Pegasi 51, CH-1290 Versoix, Switzerland  Devina
Misra Departement d’Astronomie, Université de Genève, Chemin Pegasi 51,
CH-1290 Versoix, Switzerland  Kyle A. Rocha Department of Physics and
Astronomy, Northwestern University, 2145 Sheridan Ave, Evanston, IL 60208, USA
Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA),
Northwestern University, 1800 Sherman Ave, Evanston, IL 60201, USA  Philipp M.
Srivastava  Electrical and Computer Engineering, Northwestern University, 2145
Sheridan Road, Evanston, IL 60208, USA  Center for Interdisciplinary
Exploration and Research in Astrophysics (CIERA), Northwestern University,
1800 Sherman Ave, Evanston, IL 60201, USA  Meng Sun Center for
Interdisciplinary Exploration and Research in Astrophysics (CIERA),
Northwestern University, 1800 Sherman Ave, Evanston, IL 60201, USA  Zepei Xing
Departement d’Astronomie, Université de Genève, Chemin Pegasi 51, CH-1290
Versoix, Switzerland  Emmanouil Zapartas  IAASARS, National Observatory of
Athens, Penteli, 15236, Greece

###### Abstract

Mass measurements from low-mass black hole X-ray binaries (LMXBs) and radio
pulsars have been used to identify a gap between the most massive neutron
stars (NSs) and the least massive black holes (BHs). BH mass measurements in
LMXBs are typically only possible for transient systems: outburst periods
enable detection via all-sky X-ray monitors, while quiescent periods enable
radial-velocity measurements of the low-mass donor. We quantitatively study
selection biases due to the requirement of transient behavior for BH mass
measurements. Using rapid population synthesis simulations (COSMIC), detailed
binary stellar-evolution models (MESA), and the disk instability model of
transient behavior, we demonstrate that transient-LMXB selection effects
introduce observational biases, and can suppress mass-gap BHs in the observed
sample. However, we find a population of transient LMXBs with mass-gap BHs
form through accretion-induced collapse of a NS during the LMXB phase, which
is inconsistent with observations. These results are robust against variations
of binary evolution prescriptions. The significance of this accretion-induced
collapse population depends upon the maximum NS birth mass
$M_{\mathrm{NS,birth-max}}$. To reflect the observed dearth of low-mass BHs,
COSMIC and MESA models favor $M_{\mathrm{NS,birth-max}}\lesssim 2M_{\odot}$.
In the absence of further observational biases against LMXBs with mass-gap
BHs, our results indicate the need for additional physics connected to the
modeling of LMXB formation and evolution.

X-ray transient sources (1852); Low-mass X-ray binary stars (939); Stellar
mass black holes (1611); Stellar evolutionary models (2046)

††software: astropy (Robitaille et al., 2013; Price-Whelan et al., 2018);
COSMIC (Breivik et al., 2020); matplotlib (Hunter, 2007); MESA (Paxton et al.,
2011, 2013, 2015, 2018, 2019); numpy (van der Walt et al., 2011); pandas (Wes
McKinney, 2010; The pandas development team, 2020); POSYDON (Fragos et al.,
2023); scipy (Virtanen et al., 2020)

##  1 Introduction

Black hole X-ray binaries (XRBs) are composed of a black hole (BH) accreting
mass from a non-degenerate donor star: either through Roche lobe overflow
(RLO) or captured winds. There are currently $>50$ candidate BH XRB systems
(Remillard & McClintock, 2006; Corral-Santana et al., 2016), broadly divided
by donor mass into high-mass XRBs (HMXB) and low-mass XRBs (LMXB). HMXBs host
donor stars of masses $M_{\mathrm{donor}}\gtrsim 5M_{\odot}$, and
predominately transfer mass to the compact object via strong winds, while
LMXBs host donor stars of masses $M_{\mathrm{donor}}\lesssim 2M_{\odot}$ and
transfer mass via RLO. XRBs offer powerful insights into binary evolution,
compact objects, and accretion disks. The XRB Cygnus X-1 provided the first
evidence for the existence of a BH (Bolton, 1972), and since then X-ray and
optical-infrared observations have allowed the identification of tens of
transient BH systems (Corral-Santana et al., 2016).

In cases of RLO (most commonly in LMXBs) mass-transferred from the donor star
forms an accretion disk around the BH. Depending on the structure and
composition of the system, a thermal instability may form within the disk and
generate transient X-ray emission, cycling through periods of bright outburst
and dim quiescence (Cannizzo et al., 1982; King et al., 1996; Lasota et al.,
2008). Transient systems are prime targets for follow-up radial velocity
measurements of the donor star: periods of bright outburst make the system
discoverable in X-rays, while periods of quiescence allow for inference of the
accretor’s mass via uncontaminated radial velocity measurements of the donor
star.

Mass measurements provide key insights into the formation and evolution of
BHs, and the majority of Galactic BHs with mass estimates are transient XRBs.
Exceptions include the X-ray faint binary VFTS 243, which Shenar et al. (2022)
proposed is composed of a $25M_{\odot}$ O-type star and a $>9M_{\odot}$ BH, as
well as the proposed non-interacting low-mass binary companions to 2MASS
J05215658+4359220 (Thompson et al., 2019) and V723 Mon (Jayasinghe et al.,
2021); however, V723 Mon is consistent with a stripped-star companion instead
of a compact object (El-Badry et al., 2022). These systems are outnumbered by
the $\sim 20$ transient LMXBs currently with dynamical mass measurements
(Corral-Santana et al., 2016). Using the mass measurements of the observed
LMXB population, prior studies have attempted to constrain the underlying
stellar-mass BH distribution (Bailyn et al., 1998; Özel et al., 2010; Farr et
al., 2011). Each of these investigations identified a gap between the lowest
mass BHs ($\gtrsim 5M_{\odot}$) and the maximum mass ($2$–$3M_{\odot}$) for
non-rotating neutron stars (NSs; Rhoades & Ruffini, 1974; Kalogera & Baym,
1996; Müller & Serot, 1996; Özel & Freire, 2016; Margalit & Metzger, 2017; Ai
et al., 2020; Shao et al., 2020; Raaijmakers et al., 2021); this feature is
commonly referred to as the _lower mass gap_.

Understanding the nature of the purported lower mass gap offers key insights
into the properties of degenerate matter and supernova explosion mechanisms
(Fryer et al., 2012; Mandel & Müller, 2020; Zevin et al., 2020; Liu et al.,
2021; Patton et al., 2022). If instability growth and launch of the supernova
proceeds on rapid timescales ($\sim 10~{}\mathrm{ms}$), accretion onto proto-
NSs is suppressed, and simulations predict a paucity of compact objects born
with masses between $2$–$5M_{\odot}$ (Fryer et al., 2012; Belczynski et al.,
2012; Fryer et al., 2022); however, if instability growth is delayed
(timescales greater than $\sim 200~{}\mathrm{ms}$), simulations predict a
continuous compact object birth mass distribution. Motivated by the observed
Galactic LMXB sample, prior studies have typically favored the rapid
instability timescale, but an underlying assumption of these studies is that
BH mass measurements are unbiased.

Recent detections of BH masses in systems other than LMXBs complicate the
picture of the lower mass gap. The binary sources of gravitational waves
GW190814 and GW200210_092254 are each inferred to have a source containing a
compact object lying in the mass gap: $2.59^{+0.08}_{-0.09}M_{\odot}$ and
$2.83^{+0.47}_{-0.42}M_{\odot}$, respectively (Abbott et al., 2020, 2021b).
The mass estimates for the compact-object companions to 2MASS
J05215658+4359220 and V723 Mon (assuming a compact object interpretation) also
lie within the mass gap, $3.3_{-0.7}^{+2.8}M_{\odot}$ and $(3.04\pm
0.06)M_{\odot}$ respectively (Thompson et al., 2019; Jayasinghe et al., 2021).
These measurements are in apparent contradiction with the existence of a lower
mass gap, and must be reconciled with the Galactic LMXB BH mass distribution.

Understanding the nature of detection biases in the Galactic LMXB population
is critical to interpreting the observed BH mass distribution. Through
approximate scaling of the observed XRB sample, Özel et al. (2010) argued
detection biases due to outburst flux cannot account for the mass gap. Jonker
et al. (2021) noted that one potential selection effect comes from the fact
that objects with high extinction are more difficult to obtain follow-up
radial velocity measurements for. Assuming an anti-correlation between
supernova kick magnitude and compact object mass (e.g., Fryer et al., 2012;
Atri et al., 2019), this preference for LMXBs off the Galactic plane may bias
against mass inference of higher mass BHs. While this effect alone would not
produce a mass gap as an observational selection effect, it may still shape
the observed BH mass distribution. Kreidberg et al. (2012) also demonstrated
that assuming zero or constant emission from the accretion flow can lead to
systematic overestimates of compact object masses. The role of detection
biases in the Galactic LMXB population, particularly the requirement of
transient behavior for BH mass measurement, remains unclear.

Here we adopt a forward modeling approach to assess selection biases:
generating samples of synthetic binary systems, evolving the systems forward
in time, and inferring the detectable population using the disk instability
model of LMXB transient behavior. With these methods, we investigate the
impact of observational biases on the lower mass gap and the underlying BH
mass distribution. Through our analysis, we uncover that the observed gap has
implications for the maximum NS mass at birth ($M_{\mathrm{NS,birth-max}}$).

In Section 2.1 and Section 2.3, we describe our population synthesis (COSMIC;
Breivik et al., 2020) and stellar evolution (MESA; Paxton et al., 2011, 2013,
2015, 2018, 2019) methods, respectively; the population synthesis Milky Way
model is outlined in Section 2.2. In Section 2.4 we outline our treatment of
LMXB transient behavior, and in Section 2.5 we present the formulation of
detection probabilities. In Section 2.6 we detail how we draw together a
synthetic population. We describe our results in Section 3, discuss their
implications in Section 4, and summarize our conclusions in Section 5.

##  2 LMXB Simulations

Here we outline the methods for simulating LMXB samples with rapid binary
population synthesis (COSMIC; Section 2.1) and binary stellar-evolution models
(MESA; Section 2.3). While MESA provides detailed modeling, the computational
efficiency of COSMIC allows us to explore uncertain aspects of binary
evolution physics. Taken together, the suite of proof-of-concept COSMIC
simulations contextualize the results of improved mass-transfer modeling with
MESA, by marginalizing over uncertain aspects of binary evolution physics.

###  2.1 Population models

Binary population synthesis offers valuable insights into binary evolution
channels and the expected distributions of target populations. We employ the
rapid binary population synthesis code COSMIC (Breivik et al., 2020) to
generate Milky Way LMXB populations. Originating from BSE (Hurley et al.,
2002), COSMIC approximates single-star evolution using fitting formulae to
detailed models and applies prescriptions for binary evolution to these
formulae. To complement these simulations, we also consider a three-
dimensional grid (in compact object mass $M_{\mathrm{CO}}$, donor mass
$M_{\mathrm{donor}}$, and orbital period $P$ space) of LMXBs using MESA
(Section 2.3). The computational efficiency of COSMIC allows us to explore
uncertain aspects of binary evolution physics that would not be
computationally feasible with MESA.

Independent of contextualizing MESA simulations, verifying the extent to which
binary population synthesis methods can replicate the observed LMXB sample is
of interest. For example, Ivanova (2006) proposed that some LMXBs could be fed
by pre-main-sequence donor stars given the observed relation between LMXBs’
orbital periods and the donors’ effective temperatures, as well as the fact
that only stars more massive than $\sim 2M_{\odot}$ have a pre-main-sequence
lifetime shorter than the main-sequence lifetime of a star that explodes and
forms a BH ($\lesssim 10^{7}~{}\mathrm{yr}$). However, no rapid population
synthesis codes simulate pre-main-sequence interactions. These constraints on
the binary population suppress the formation of LMXBs with massive BHs and
shorter period binaries. Any rapid population synthesis investigation of LMXBs
that does not consider pre-main-sequence interactions is impacted by these
limitations. While population synthesis methods likely cannot reproduce all
aspects of the observed LMXB sample, quantifying these discrepancies is vital
to developing improved population synthesis engines. Moreover, we are
employing population synthesis as a proof-of-concept survey of binary
evolution prescriptions: rather than replicating all aspects of the observed
LMXB distribution, we are interested in uncovering any significant
dependencies between the synthetic observed LMXB sample and uncertain aspects
of binary evolution physics.

Table 1: Milky Way Component Parameters.

Component | AgeaaThe adopted Milky Way age is $\mathcal{A}_{\mathrm{MW}}=11$ Gyr. | Metallicity | Star formation history | Stellar mass | Spatial distribution parameters
---|---|---|---|---|---
| (Gyr) | $[\mathrm{Fe}/\mathrm{H}]$ |  | ($10^{10}M_{\odot}$) |
Thin disk |  $\mathcal{A}_{\mathrm{MW}}-1$ Gyr | 0.0 | Continuous | 4.32 |  $z_{\mathrm{disk}}=0.3$ kpc, $R_{\mathrm{disk}}=2.90$ kpc
Thick disk | $\mathcal{A}_{\mathrm{MW}}$ | $-0.8$ |  $1$ Gyr burst | 1.44 |  $z_{\mathrm{disk}}=0.9$ kpc, $R_{\mathrm{disk}}=3.31$kpc
Bulge |  $\mathcal{A}_{\mathrm{MW}}-1$ Gyr | 0.0 |  $1$ Gyr burst | 0.89 |  $a=1.8$, $r_{0}=0.075$ kpc, $r_{\mathrm{cut}}=2.1$ kpc, $q=0.5$

Note. — Component masses and spatial distributions follow McMillan (2011).
Star formation histories and metallicities are adopted from Robin et al.
(2003).

To generate Galactic LMXB population synthesis samples, we simulate Milky Way
stellar populations for a range of binary evolution prescriptions. For each
COSMIC population:

  1. 1.

The maximum NS mass is set to $3M_{\odot}$ (Section 1; Rhoades & Ruffini,
1974; Kalogera & Baym, 1996).

  2. 2.

Compact object birth masses from core-collapse supernovae are assigned
following the _Delayed_ prescription of Fryer et al. (2012); this differs from
the alternative _Rapid_ prescription in the assumed timescale of instability
growth and launch of the supernova. The Delayed remnant prescription produces
a continuous mass distribution, while the Rapid prescription produces a gap
between the most massive NSs and the least massive BHs. The Delayed
prescription is adopted to explore whether LMXB selection biases can hide an
underlying population of low-mass BHs.

  3. 3.

Compact object natal kicks are drawn from a bimodal distribution: standard
iron core-collapse supernova kicks are drawn from a Maxwellian distribution
with a dispersion of $265$ $\mathrm{km\,s^{-1}}$ (Hobbs et al., 2005), while
kicks for electron-capture supernovae and ultra-stripped supernovae are drawn
according to a Maxwellian distribution with a dispersion of
$20~{}\mathrm{km\,s^{-1}}$ (Giacobbo & Mapelli, 2019). Natal kicks are reduced
by a factor of $1-f_{\mathrm{fb}}$, where $f_{\mathrm{fb}}$ is the fraction of
the ejected supernova mass that will fall back onto the proto-compact object
(Fryer et al., 2012).

  4. 4.

Initial conditions (e.g., eccentricity, orbital period, primary mass, and mass
ratio) are independently drawn following Sana et al. (2012).

We investigate model dependence on six uncertain factors of binary evolution:

  1. 1.

Stars without a distinct core–envelope boundary (e.g., stars on the
Hertzsprung gap) that instigate a common envelope (CE) event are either
assumed to survive the CE (_Optimistic_) or are assumed to merge
(_Pessimistic_ ; Belczynski et al., 2008).

  2. 2.

The CE efficiency $\alpha$ parameterizes the transfer of orbital energy to the
envelope, i.e., how easily the CE is unbound from the system (Webbink, 1984;
de Kool, 1990); higher values of $\alpha$ lead to longer period binaries post-
CE. Prior studies using one-dimensional simulations (Fragos et al., 2019) and
comparison to gravitation-wave observations (Giacobbo & Mapelli, 2018;
Santoliquido et al., 2021; Zevin et al., 2021) have favored higher CE
efficiencies ($\alpha\sim 5$). To explore these higher values, $\alpha$ is set
to $\alpha=1.0$ or $\alpha=5.0$.

  3. 3.

The minimum ZAMS mass ratio between binary members is set to $q>0.01$, where
$q\equiv M_{\mathrm{donor}}/M_{\mathrm{primary}}$, or is defined as a function
of primary mass (_Lifetime limited_). For the Lifetime limited case, $q$ is
restricted such that the pre-main-sequence lifetime of the secondary is
shorter than the main-sequence lifetime of the primary.

  4. 4.

In addition to gravitational radiation, close-in binaries (e.g., cataclysmic
variables and LMXBs) can efficiently lose orbital angular momentum via
magnetic braking of a tidally coupled magnetic wind. Magnetic braking is
implemented following Hurley et al. (2002, hereafter H02) or Ivanova & Taam
(2003, hereafter IT03).

  5. 5.

The amount of mass accreted during Roche-lobe overflow is either fixed to
$50\%$ efficiency (Belczynski et al., 2008) or is a function of the accretor’s
type (_State dependent_). For the State dependent prescription, the accretion
efficiency during Roche-lobe overflow is either (i) limited to ten-times the
thermal rate if the accretor is a main-sequence, Hertzsprung gap or core
helium burning star, or (ii) unlimited if the accretor is a giant branch,
early asymptotic giant branch or asymptotic giant branch star. Compact objects
also experience conservative mass transfer, provided the mass-transfer rate is
sub-Eddington. Accretion onto a NS is treated independently of whether the
system is transient. While outburst events during accretion onto a white dwarf
can lead to a loss of mass, only a small of fraction of accreting NSs show
evidence of mass loss during bursts (Degenaar et al., 2018); this difference
between white dwarf and NS accretors is likely attributable to the deeper
potential wells of NSs.

  6. 6.

The maximum NS birth mass from core-collapse supernova $M_{\mathrm{NS,birth-
max}}$ is of particular interest to the formation and evolution of LMXBs;
lowering $M_{\mathrm{NS,birth-max}}$ suppresses the formation of LMXBs through
accretion-induced collapse (AIC) of a NS. To vary $M_{\mathrm{NS,birth-max}}$
while preserving computational feasibility, $M_{\mathrm{NS,birth-max}}$ is
fixed to $3M_{\odot}$ within COSMIC and is lowered in post-processing by
removing systems with a NS birth mass above a given threshold. This method of
imposing a maximum NS birth mass is a non-physical toy model. The goal is to
demonstrate how a truncated NS birth mass distribution impacts LMXB
demographics.

To simulate stellar populations analogous to the Milky Way, we assign birth
times, spatial positions, and metallicities following a three-component Milky
Way model (Section 2.2).

###  2.2 Milky Way model

For a given set of binary evolution prescriptions, we generate a synthetic
Galactic sample of LMXBs. To do so, we resample a given COSMIC population and
assign system properties in accordance with the Milky Way (e.g., star-
formation history and spatial distribution). The Milky Way is approximated as
a composite of three sub-populations: thin disk, thick disk, and bulge.

The star-formation history and metallicity of each component is adopted from
Robin et al. (2003). For the thick disk and bulge, star formation is modeled
as a $1~{}\mathrm{Gyr}$ continuous burst, beginning at the birth of the Milky
Way and $1~{}\mathrm{Gyr}$ after the birth of the Milky Way, respectively.
Star formation for the thin disk is continuous from $1~{}\mathrm{Gyr}$ after
the birth of the Milky Way to the present. Stars in the thin disk and bulge
have solar metallicity, while stars in the thick disk have $15\%$ of solar
metallicity. System birth times are relative to a Milky Way age of
$\mathcal{A}_{\mathrm{MW}}=11~{}\mathrm{Gyr}$. LMXBs’ positions within the
galaxy are fixed to their birthplaces; we do not include the effects of SN
kicks on systemic velocity.

The component masses and spatial distributions for the Milky Way model follow
McMillan (2011). The thin and thick disk spatial distributions are treated as
an exponential,

$\displaystyle\rho_{\mathrm{disk}}(R,z)\propto\exp\left(-\frac{\mid
z\mid}{z_{\mathrm{disk}}}-\frac{R}{R_{\mathrm{disk}}}\right),$ (1)

where $z_{\mathrm{disk}}$ is the scale-height and $R_{\mathrm{disk}}$ is the
scale length; for the thin disk, $z_{\mathrm{disk}}=0.3~{}\mathrm{kpc}$ and
$R_{\mathrm{disk}}=2.9~{}\mathrm{kpc}$; for the thick disk,
$z_{\mathrm{disk}}=0.9~{}\mathrm{kpc}$ and
$R_{\mathrm{disk}}=3.31~{}\mathrm{kpc}$ (McMillan, 2011). The bulge spatial
distribution is

$\displaystyle\rho_{\mathrm{bulge}}(R,z)$
$\displaystyle\propto\frac{\exp\left[-(r^{\prime}/r_{\mathrm{cut}})^{2}\right]}{(1+r^{\prime}/r_{0})^{a}},$
(2) $\displaystyle r^{\prime}$ $\displaystyle=\sqrt{R^{2}+(z/q)^{2}},$ (3)

where $a=1.8$, $r_{0}=0.075~{}\mathrm{kpc}$,
$r_{\mathrm{cut}}=2.1~{}\mathrm{kpc}$, and the axial ratio $q=0.5$. Each
component is treated as axisymmetric.

The parameters of the Milky Way model are summarized in Table 1.

###  2.3 MESA models

Rapid binary population synthesis enables the efficient modeling of
statistically significant stellar samples, but such methods have drawbacks.
Rapid binary population synthesis codes, including COSMIC, assume binary
members have the same properties as a single star with matching mass and
metallicity at thermal equilibrium; this can lead to systematic biases in the
binaries’ evolution (Gallegos-Garcia et al., 2021; Fragos et al., 2023). As
outlined above, the COSMIC simulations serve as a proof-of-concept survey over
uncertain aspects of binary evolution physics and contextualize a smaller grid
of binaries modeled in greater detail with MESA. Here we outline the LMXB grid
evolved using MESA.

We consider a three-dimensional grid of binaries consisting of a compact
object and a hydrogen-rich main-sequence star. Adopted from Fragos et al.
(2023), the grid consists of $1375$ binaries, initialized uniformly in
$\log_{10}(M_{\mathrm{CO}}/M_{\odot})$,
$\log_{10}(M_{\mathrm{donor}}/M_{\odot})$, and $\log_{10}(P/\mathrm{day})$.
Compact object masses span the range $1.1$–$10M_{\odot}$, donor star masses
span $0.5$–$3M_{\odot}$, and orbital periods span $1.26$–$3162$ days. The grid
includes LMXBs with mass-gap BH, which allows us to investigate the
observability of such sources. The maximum NS mass is set to $3M_{\odot}$ (as
in the COSMIC simulations). All systems are initialized at solar metallicity,
and the donor stars are seeded with a rotational period equal to their orbital
period (Fragos et al., 2023).

To facilitate comparison between the population synthesis and MESA methods, we
initialize and evolve grids of binaries with COSMIC at the same initial
conditions as the MESA grid. We consider the same combinations of binary
evolution prescriptions as Section 2.1, with the exception of CE efficiency
and minimum ZAMS mass ratio (which do not apply to this phase of evolution).

Figure 1: X-ray outburst light curves for a typical LMXB system for a range of
BH masses and transitions to RIA. On the _left_ , we assume no transition to
RIA and vary the BH mass with fixed orbital period and companion mass. On the
_right_ , we present the sharp and smooth transitions to RIA for fixed BH
mass. Both panels assume $M_{\mathrm{donor}}=0.7M_{\odot}$,
$P=10~{}\mathrm{hour}$, and $d=8~{}\mathrm{kpc}$. The shaded gray region
highlights were the X-ray luminosity drops below the detectable limit
($10~{}\mathrm{mCrab}$). All luminosities have been converted to the
$2$–$10~{}\mathrm{keV}$ range by dividing the bolometric luminosity by
$f_{\mathrm{corr}}=4$.

###  2.4 Transient LMXB Behavior

####  2.4.1 Instability Cycle

LMXB accretion disks are complex rotating, turbulent, irradiated systems.
Despite this complexity, the outburst cycles of LMXBs are reasonably well
approximated by the disk instability model (King et al., 1996; Dubus et al.,
1999; Lasota, 2001, 2016). This simplified model of accretion physics assumes
(i) a constant mass-transfer rate, determined by the secular mass-transfer
rate from the binary evolution calculations, (ii) that each outburst is
identical, a consequence of the assumed accretion of the entire disk during an
outburst, and (iii) that the mass transfer is conservative when sub-Eddington.
Observational work shows these assumptions are imperfect (e.g., the observed
bimodality of outburst durations; Lasota, 2001, Section 4.1). The disk
instability model is adopted for this work given that it is quantitatively
tractable, and efforts to improve the model are beyond the scope of this
paper.

Under the disk instability model, a disk is thermally stable provided
radiative cooling varies faster with temperature than viscous heating (Frank
et al., 1992; Dubus et al., 1999). As a result, a rapid change in opacity
instigates a thermal instability in the disk (e.g., if a disk is composed of
ionized hydrogen, and the temperature falls low enough for hydrogen to
recombine). If thermally unstable, a disk will cycle through hot outburst
states and cold quiescent states.

For a given LMXB, the conditions for thermal equilibrium of the accretion disk
(typically presented in the disk surface density $\Sigma$ and effective
temperature $T_{\mathrm{eff}}$ plane) can be numerically calculated as a
function of disk radius. Equilibrium solutions in the
$\Sigma$–$T_{\mathrm{eff}}$ plane (_S-curves_) are composed of a stable hot
branch and a stable cold branch connected by an instability strip (Lasota,
2001, Section 3). LMXBs where the mass-transfer rate results in a
$T_{\mathrm{eff}}$ on the instability strip will be thermally unstable,
leading the disk to cycle between the hot and cold branches. The inflection
points of the S-curve correspond to critical accretion rates that bound the
instability strip. By calculating S-curves for a range of disk radii and
binary system parameters (e.g., orbital period, donor mass, and accretor
mass), the cold branch and hot branch bounding critical accretion rates can be
interpolated for a given LMXB.

In principle, stable LMXB accretion disks can occupy the cold branch or the
hot branch of the S-curve. However, binary evolution favors mass-transfer
rates above the cold branch critical mass-transfer rate (Lasota, 2001).
Systems may occupy the cold branch late in their evolution, when their donor
stars become degenerate around the point that orbital periods begin to
increase (Knigge et al., 2011); in such cases, the outbursts would not
contribute significantly to the detectable population of sources, because the
outbursts would be short, faint, and infrequent. In this study, a LMXB is
considered a transient source if its mass-transfer rate falls below the hot
branch critical accretion rate
$\dot{M}_{\mathrm{donor}}<\dot{M}_{\mathrm{crit}}$.

For systems with hydrogen-dominated donor stars, the hot branch critical mass-
transfer rate for X-ray irradiated disks is (Dubus et al., 1999)

$\displaystyle\dot{M}_{\mathrm{crit}}=1.5\times 10^{15}$
$\displaystyle\left(\frac{M_{\mathrm{BH}}}{M_{\odot}}\right)^{-0.4}\left(\frac{R_{\mathrm{disk}}}{10^{10}~{}\mathrm{cm}}\right)^{2.1}$
$\displaystyle\times\left(\frac{C}{5\times
10^{-4}}\right)~{}\mathrm{g~{}s^{-1}},$ (4)

where $M_{\mathrm{BH}}$ is the BH mass, $R_{\mathrm{disk}}$ is the outer
radius of the accretion disk, and $C$ is a disk structure constant (which we
fix to $C=5\times 10^{-4}$; Dubus et al., 1999). For helium-dominated donors,
where the opacity of the inflow and accretion disk are markedly different from
the hydrogen-dominated case, the hot branch irradiated disk critical mass-
transfer rate is (Lasota et al., 2008)

$\displaystyle\dot{M}_{\mathrm{crit}}$ $\displaystyle=2.1\times
10^{16}\left(\frac{C_{\mathrm{He}}}{10^{-3}}\right)^{-0.22}$
$\displaystyle\times\left(\frac{\alpha}{0.1}\right)^{-0.03-0.01\log(C_{\mathrm{He}}/10^{-3})}$
$\displaystyle\times\left(\frac{R_{\mathrm{disk}}}{10^{10}~{}\mathrm{cm}}\right)^{2.51-0.05\log(C_{\mathrm{He}}/10^{-3})}$
$\displaystyle\times\left(\frac{M_{\mathrm{BH}}}{M_{\odot}}\right)^{-0.74+0.05\log(C_{\mathrm{He}}/10^{-3})}~{}\mathrm{g\,s}^{-1},$
(5)

where $C_{\mathrm{He}}$ is a disk structure constant for helium-dominated
donors and $\alpha$ is a viscosity parameter (we fix $C_{\mathrm{He}}=10^{-3}$
and $\alpha=0.1$; Lasota et al., 2008). We adopt the standard formulation of
disk size: $R_{\mathrm{disk}}=70\%$ of the BH’s Roche lobe radius (King et
al., 1997).

If the mass-transfer rate is below the critical rate, the system will undergo
a cycle of outburst and quiescent phases of duration $T_{\mathrm{outburst}}$
and $T_{\mathrm{quiescence}}$, respectively. If the mass-transfer rate is
above the critical rate, the system will be a persistent X-ray source. For
transient LMXBs, the period of the instability cycle is
$T_{\mathrm{cycle}}\equiv T_{\mathrm{quiescence}}+T_{\mathrm{outburst}}$.
Assuming the entire disk is accreted onto the BH during an outburst event, the
duration of the quiescent phase is

$T_{\mathrm{quiescence}}=\frac{M_{\mathrm{disk-
max}}}{\dot{M}_{\mathrm{donor}}},$ (6)

where $\dot{M}_{\mathrm{donor}}$ is the mass-transfer rate from the donor
star. Adopting the critical surface density $\Sigma_{\mathrm{crit}}$ of
Cannizzo et al. (1988), and assuming the surface density of the quiescent disk
prior to in-fall is comparable to $\Sigma_{\mathrm{crit}}$ (Dubus et al.,
2001), the maximum accretion disk mass is

${M}_{\mathrm{disk-
max}}\approx\int_{R_{\mathrm{ISCO}}}^{R_{\mathrm{disk}}}2\pi\Sigma_{\mathrm{crit}}(R)R\,\mathrm{d}R,$
(7)

where $R_{\mathrm{ISCO}}$ is the radius of the innermost stable circular
orbit.

The duration of the outburst $T_{\mathrm{outburst}}$ is the total time
necessary for the disk to accrete onto the BH. Since the mass in-fall rate
during the instability $\dot{M}_{\mathrm{ins}}(t)$ varies during the outburst,
$T_{\mathrm{outburst}}$ cannot be approximated analogously to
$T_{\mathrm{quiescence}}$ in Eq. (6). Following King & Ritter (1998), the mass
in-fall rate $\dot{M}_{\mathrm{ins}}(t)$ is treated as an exponential decay
followed by a linear decay.

The bolometric outburst light curve is related to the in fall rate by

$L_{\mathrm{bol}}(t)=\eta(t)c^{2}\dot{M}_{\mathrm{ins.}}(t),$ (8)

where $\eta(t)$ is the radiative efficiency and $c$ is the speed of light.
Adopting the methods of Knevitt et al. (2014), we consider two potential
formulations of $\eta(t)$, corresponding to two different transitions to
radiatively inefficient accretion (RIA): a _sharp_ transition,

$\eta=\begin{cases}0,&L_{\mathrm{bol}}<\nu L_{\mathrm{Edd}}\\\
0.1,&L_{\mathrm{bol}}\geq\nu L_{\mathrm{Edd}}\end{cases},$ (9)

or a _smooth_ transition,

$\eta=\begin{cases}0.1\left(\dfrac{\dot{M}_{\mathrm{ins}}}{\nu\dot{M}_{\mathrm{Edd}}}\right),&L_{\mathrm{bol}}<\nu
L_{\mathrm{Edd}}\\\ 0.1,&L_{\mathrm{bol}}\geq\nu L_{\mathrm{Edd}}\end{cases},$
(10)

where we either fix $\nu$ to $0$ (no transition to RIA) or $0.05$; here
$\dot{M}_{\mathrm{Edd}}$ is the Eddington mass-transfer rate, and
$L_{\mathrm{Edd}}$ is the corresponding Eddington luminosity.

In this approximate outburst model, including a transition to RIA shortens the
X-ray outburst lifetime without changing the peak luminosity of the outburst.
To highlight this dependence, outburst light curves for a typical LMXB system
under varying BH masses and transitions to RIA are presented in Figure 1.

Since the outburst luminosity may drop below a detectable level before the
full $T_{\mathrm{outburst}}$ has elapsed, we introduce the effective outburst
time $T^{\prime}_{\mathrm{outburst}}$, the total amount of time that the X-ray
outburst flux is above an observable threshold. For a limiting X-ray flux of
$F_{\mathrm{fiducial}}$, the outburst is observable if

$F_{\mathrm{fiducial}}<\frac{L_{\mathrm{bol}}(t)}{4\pi
d^{2}f_{\mathrm{corr}}},$ (11)

where $f_{\mathrm{corr}}$ converts the bolometric flux into the observable
range and $d$ is the distance to the system. We adopt $f_{\mathrm{corr}}=4$
and $F_{\mathrm{fiducial}}=10~{}\mathrm{mCrab}$ (Section 2.5). For COSMIC, $d$
is drawn for each system from the spatial distributions described in Section
2.2; for the MESA grid, $d$ is fixed to $8$ kpc (Section 2.6).

In Section 2.5, LMXBs’ effective outburst times
$T^{\prime}_{\mathrm{outburst}}$ and cycle lifetimes $T_{\mathrm{cycle}}$ are
mapped to probabilities of detection with an all-sky X-ray monitor.

###  2.5 Detection Probabilities

Transient X-ray binaries have predominately been discovered by all-sky
monitors (ASMs). Although the observed sample of transient LMXBs has been
built up by a myriad of instruments, we treat the ASM aboard the Rossi X-ray
Timing Explorer (RXTE; Bradt et al., 1993) as our benchmark detector. Active
from 1996 to 2012, the RXTE ASM observed in the $1.5$–$12~{}\mathrm{keV}$ band
with a highly stochastic pointing pattern that typically scanned a given
source 5–10 times per day (Levine et al., 1996). RXTE discovered a significant
portion of the transient LMXB sample (Wen et al., 2006). The detection
requirements of the RXTE ASM can be approximated as a flux limit of
$F_{\mathrm{fiducial}}=10~{}\mathrm{mCrab}$ for a minimum emission time of
$T_{\mathrm{fiducial}}=1$ day.

For a transient LMXB source to be identified by an all-sky survey and be
eligible for mass inference, it must (i) be bright enough during outburst to
be discovered, and (ii) transition from an outburst to a quiescent state
during the survey. The latter requirement enables mass inference via radial-
velocity follow-up of the low-mass donor. We neglect inclination effects and
potential biases induced by the position of the Sun.

If an outburst occurs during an ASM survey, the probability the outburst is
detected by the ASM is (Knevitt et al., 2014)

$\displaystyle p_{\mathrm{discover}}=$
$\displaystyle\min\left\\{1,\frac{T^{\prime}_{\mathrm{outburst}}}{T_{\mathrm{fiducial}}}\right\\},$
(12)

where $T^{\prime}_{\mathrm{outburst}}$ is the effective outburst time: the
total amount of time that the X-ray outburst is brighter than
$F_{\mathrm{fiducial}}$ (Section 2.4).

In addition to the requirement on the system’s outburst flux, a LMXB must also
transition from outburst to quiescence during the survey in order to be
classified as a transient. The probability of such a transition is

$p_{\mathrm{transition}}=\min\left\\{1,\frac{T_{\mathrm{survey}}}{T_{\mathrm{cycle}}}\right\\},$
(13)

where $T_{\mathrm{cycle}}$ is the period of the instability cycle,
$T_{\mathrm{cycle}}\equiv T_{\mathrm{quiescence}}+T_{\mathrm{outburst}}$, and
$T_{\mathrm{survey}}$ is the duration of the X-ray monitor’s survey.

Considering both the probability of being seen in outburst and the probability
of observing a transition, it follows that the probability of a transient LMXB
source being identified is

$\displaystyle
p_{\mathrm{detect}}=p_{\mathrm{discover}}(T^{\prime}_{\mathrm{outburst}})\times
p_{\mathrm{transition}}(T_{\mathrm{cycle}}).$ (14)

Motivated by RXTE, we adopt $T_{\mathrm{survey}}=15$ yr. For persistent LMXBs,
we set $p_{\mathrm{detect}}=0$.

Using Eq. (14), the probability of identifying a LMXB as a transient source
can be inferred from the outburst time $T^{\prime}_{\mathrm{outburst}}$ and
cycle lifetime $T_{\mathrm{cycle}}$; $T^{\prime}_{\mathrm{outburst}}$ and
$T_{\mathrm{cycle}}$ can in-turn be inferred from the evolutionary state of
the LMXB, e.g., mass-transfer rate, BH mass, donor mass, and orbital period.
In Section 2.6, we leverage Eq. (14) to generate synthetic samples of observed
LMXBs for both COSMIC and MESA.

Table 2: Population Models

CE survival | $\alpha$ |  $q$ limit | Magnetic braking | Accretion limit |  $f_{\mathrm{low}}(M_{\mathrm{NS,birth}}<1.5M_{\odot})$aaThe fraction $f_{\mathrm{low}}(p_{\mathrm{detect}})$ is defined relative to $M_{\mathrm{bound}}=4.5M_{\odot}$.  | $f_{\mathrm{low}}(M_{\mathrm{NS,birth}}<3.0M_{\odot})$
---|---|---|---|---|---|---
Optimistic | 1.0 | Lifetime limited | IT03 | State dependent | $0.39^{+0.04}_{-0.04}$ | $0.89^{+0.02}_{-0.03}$
Optimistic | 1.0 | Lifetime limited | IT03 |  $50\%$ efficiency | $0.59^{+0.02}_{-0.04}$ | $0.67^{+0.02}_{-0.04}$
Optimistic | 1.0 | Lifetime limited | H02 | State dependent | $0.79^{+0.02}_{-0.02}$ | $0.955^{+0.005}_{-0.006}$
Optimistic | 1.0 | Lifetime limited | H02 |  $50\%$ efficiency | $0.58^{+0.03}_{-0.03}$ | $0.72^{+0.02}_{-0.02}$
Optimistic | 1.0 | 0.01 | IT03 | State dependent | $0.37^{+0.02}_{-0.02}$ | $0.80^{+0.01}_{-0.02}$
Optimistic | 1.0 | 0.01 | IT03 |  $50\%$ efficiency | $0.58^{+0.03}_{-0.02}$ | $0.65^{+0.04}_{-0.02}$
Optimistic | 1.0 | 0.01 | H02 | State dependent | $0.603^{+0.013}_{-0.009}$ | $0.880^{+0.006}_{-0.007}$
Optimistic | 1.0 | 0.01 | H02 |  $50\%$ efficiency | $0.56^{+0.01}_{-0.02}$ | $0.67^{+0.01}_{-0.01}$
Optimistic | 5.0 | Lifetime limited | IT03 | State dependent | $0.44^{+0.02}_{-0.02}$ | $0.943^{+0.004}_{-0.004}$
Optimistic | 5.0 | Lifetime limited | IT03 |  $50\%$ efficiency | $0.65^{+0.02}_{-0.03}$ | $0.80^{+0.01}_{-0.01}$
Optimistic | 5.0 | Lifetime limited | H02 | State dependent | $0.56^{+0.03}_{-0.03}$ | $0.965^{+0.002}_{-0.005}$
Optimistic | 5.0 | Lifetime limited | H02 |  $50\%$ efficiency | $0.60^{+0.03}_{-0.04}$ | $0.84^{+0.01}_{-0.01}$
Optimistic | 5.0 | 0.01 | IT03 | State dependent | $0.39^{+0.02}_{-0.01}$ | $0.893^{+0.012}_{-0.008}$
Optimistic | 5.0 | 0.01 | IT03 |  $50\%$ efficiency | $0.53^{+0.01}_{-0.02}$ | $0.70^{+0.02}_{-0.01}$
Optimistic | 5.0 | 0.01 | H02 | State dependent | $0.43^{+0.01}_{-0.02}$ | $0.914^{+0.005}_{-0.004}$
Optimistic | 5.0 | 0.01 | H02 |  $50\%$ efficiency | $0.50^{+0.01}_{-0.02}$ | $0.741^{+0.009}_{-0.012}$
Pessimistic | 1.0 | Lifetime limited | IT03 | State dependent | $0.36^{+0.04}_{-0.05}$ | $0.88^{+0.01}_{-0.01}$
Pessimistic | 1.0 | Lifetime limited | IT03 |  $50\%$ efficiency | $0.57^{+0.04}_{-0.05}$ | $0.64^{+0.03}_{-0.04}$
Pessimistic | 1.0 | Lifetime limited | H02 | State dependent | $0.63^{+0.03}_{-0.03}$ | $0.948^{+0.005}_{-0.004}$
Pessimistic | 1.0 | Lifetime limited | H02 |  $50\%$ efficiency | $0.56^{+0.02}_{-0.03}$ | $0.71^{+0.04}_{-0.02}$
Pessimistic | 1.0 | 0.01 | IT03 | State dependent | $0.47^{+0.01}_{-0.04}$ | $0.837^{+0.008}_{-0.034}$
Pessimistic | 1.0 | 0.01 | IT03 |  $50\%$ efficiency | $0.56^{+0.02}_{-0.02}$ | $0.64^{+0.04}_{-0.02}$
Pessimistic | 1.0 | 0.01 | H02 | State dependent | $0.48^{+0.01}_{-0.01}$ | $0.856^{+0.003}_{-0.003}$
Pessimistic | 1.0 | 0.01 | H02 |  $50\%$ efficiency | $0.597^{+0.014}_{-0.009}$ | $0.69^{+0.02}_{-0.01}$
Pessimistic | 5.0 | Lifetime limited | IT03 | State dependent | $0.35^{+0.06}_{-0.05}$ | $0.988^{+0.002}_{-0.001}$
Pessimistic | 5.0 | Lifetime limited | IT03 |  $50\%$ efficiency | $0.63^{+0.04}_{-0.05}$ | $0.83^{+0.02}_{-0.04}$
Pessimistic | 5.0 | Lifetime limited | H02 | State dependent | $0.70^{+0.02}_{-0.02}$ | $0.984^{+0.001}_{-0.002}$
Pessimistic | 5.0 | Lifetime limited | H02 |  $50\%$ efficiency | $0.59^{+0.03}_{-0.02}$ | $0.930^{+0.006}_{-0.008}$
Pessimistic | 5.0 | 0.01 | IT03 | State dependent | $0.28^{+0.03}_{-0.03}$ | $0.928^{+0.014}_{-0.003}$
Pessimistic | 5.0 | 0.01 | IT03 |  $50\%$ efficiency | $0.455^{+0.01}_{-0.0090}$ | $0.643^{+0.009}_{-0.013}$
Pessimistic | 5.0 | 0.01 | H02 | State dependent | $0.351^{+0.012}_{-0.007}$ | $0.911^{+0.002}_{-0.002}$
Pessimistic | 5.0 | 0.01 | H02 |  $50\%$ efficiency | $0.40^{+0.03}_{-0.03}$ | $0.72^{+0.02}_{-0.01}$

Note. — Using COSMIC, a set of population synthesis simulations spanning
uncertain aspects of binary evolution physics reveals that rapid population
synthesis consistently fills the mass gap. The detection-weighted fraction of
LMXBs with a BH in the mass gap $f_{\mathrm{low}}(p_{\mathrm{detect}})$ is
presented for thirty-two exploratory COSMIC Milky Way populations (assuming no
transition to RIA and either $M_{\mathrm{NS,birth-max}}=1.5M_{\odot}$ or
$3.0M_{\odot}$). For a given population, the average and the $68\%$
uncertainty is reported from $5000$ snapshots of the evolution tracks between
$10.5~{}\mathrm{Gyr}$ and $11.5~{}\mathrm{Gyr}$, to reflect the uncertainty in
the age of the Milky Way. Each population adopts a different combination of
binary evolution prescriptions, and each row of the table corresponds to a
different COSMIC population. As outlined in Section 2.1, the prescriptions for
uncertain aspects of binary evolution physics are varied between the COSMIC
populations, including the conditions for merging during a CE event
(Optimistic versus Pessimistic), the efficiency of CE ejection ($\alpha=1$ or
$5$), the minimum mass ratio $q$ at ZAMS ($q>0.01$ or Lifetime limited, where
$q$ is restricted such that the pre-main-sequence lifetime of the secondary is
shorter than the main-sequence lifetime of the primary), the magnetic braking
implementation (following either H02 or IT03), the efficiency of accretion
during RLO (either $50\%$ efficiency or State dependent efficiency, where
compact objects have conservative mass accretion up to the Eddington limit).
For each population, the compact object birth mass distribution is fixed to
the Delayed prescription of Fryer et al. (2012), the maximum NS mass is fixed
to $3M_{\odot}$, and the binaries’ initial conditions are independently drawn
from Sana et al. (2012).

Figure 2:  The BH mass (_upper panels_) and orbital period (_lower panels_)
distributions for a selection of synthetic Galactic LMXB samples. Detection
probabilities are calculated with no transition to RIA. Each panel corresponds
to a different set of binary evolution prescriptions. For each presented
population $\alpha=5$, the minimum $q$ at ZAMS follows the Lifetime limited
prescription, and $M_{\mathrm{NS,birth-max}}=3.0M_{\odot}$. The solid lines
represent the average distribution across all $5000$ snapshots of the
evolution tracks, and the shaded regions demarcate the $16$th–$84$th
percentiles. AIC LMXBs are highlighted in purple. The detection-weighted BH
mass and orbital period distributions of LMXBs from population synthesis are
broadly insensitive to changes in the adopted binary evolution prescriptions.

###  2.6 Synthetic observed samples

To facilitate comparison between the observed LMXB sample and synthetic
populations of LMXBs, a Monte Carlo observing procedure is employed. Samples
are drawn from the evolution tracks of COSMIC and MESA based on the assumed
age of the Milky Way, and are then compared with the Galactic LMXB population
to investigate the nature of the lower mass gap.

While similar in philosophy, the synthetic observing process differs between
the COSMIC Milky Way populations and the MESA grid; the differences are
motivated by the far lower number of binaries simulated in the MESA grid. For
a given COSMIC Milky Way population, we collect $5000$ snapshots of the
population within a time window of $10.5$–$11.5~{}\mathrm{Gyr}$, informed by
our adopted Milky Way age of $\mathcal{A}_{\mathrm{MW}}=(11\pm
0.5)~{}\mathrm{Gyr}$. The time of each snapshot is independently and uniformly
drawn within the window. At every snapshot, we infer the properties of the
population via linear interpolation of the COSMIC time steps. Statistical
properties of the synthetic observed sample can then be inferred by
considering all $5000$ snapshots.

As described in Section 2.3, we also consider a three-dimensional grid of
compact object–main sequence binaries with MESA, and initialize COSMIC
populations with matching initial conditions for a range of binary evolution
prescriptions. Given the discrete distribution of initial conditions and the
lower number of binaries, synthetic observed samples are drawn from the grids
using a modified version of the Monte Carlo procedure applied to the COSMIC
populations. For a given trial, we draw a random Milky Way age
$\mathcal{A}_{\mathrm{MW}}$ from $\mathcal{U}(10.5,11.5)~{}\mathrm{Gyr}$ and
draw each binary a birth time independently from
$\mathcal{U}(0,\mathcal{A}_{\mathrm{MW}})$; this assignment of binary birth
times is analogous to continuous star formation and is based on the thin disk
component of the Milky Way model (Table 1). Linear interpolation is then
applied to infer each binary’s properties at $\mathcal{A}_{\mathrm{MW}}$.
While distances are assigned via spatial distributions for the COSMIC Milky
Way populations, a nominal distance of $8~{}\mathrm{kpc}$ is adopted for each
binary in the MESA grid, due to the low number of binaries. In total, $5000$
trails are drawn for each grid.

For a single snapshot, the number of observed transient LMXB systems is

$N_{\mathrm{transient}}=\sum_{i=1}^{N_{\mathrm{LMXB}}}p_{\mathrm{detect},i},$
(15)

where $N_{\mathrm{LMXB}}$ is the total number of LMXB systems in the
population during the snapshot. Systems composed of a BH accreting mass from a
non-degenerate donor star via RLO are considered LMXBs, which excludes
degenerate donors, to mirror the observed sample of LMXBs with BH mass
estimates.

To quantify the density of transient LMXB systems within the lower mass gap,
we define the fraction

$f_{\mathrm{low}}(p_{\mathrm{detect}})=\frac{1}{N_{\mathrm{transient}}}\sum_{i=1}^{N_{\mathrm{LMXB}}}\begin{cases}p_{\mathrm{detect},i},&M_{\mathrm{BH}}<M_{\mathrm{bound}}\\\
0,&M_{\mathrm{BH}}\geq M_{\mathrm{bound}}\end{cases}.$ (16)

Here $M_{\mathrm{bound}}$ is nominally fixed to $4.5M_{\odot}$.

The fraction $f_{\mathrm{low-intrinsic}}$ is then defined as the unweighted
fraction of LMXB systems with a BH mass less than $M_{\mathrm{bound}}$,

$f_{\mathrm{low-
intrinsic}}=\frac{1}{N_{\mathrm{LMXB}}}\sum_{i=1}^{N_{\mathrm{LMXB}}}\begin{cases}1,&M_{\mathrm{BH}}<M_{\mathrm{bound}}\\\
0,&M_{\mathrm{BH}}\geq M_{\mathrm{bound}}\end{cases}.$ (17)

This $f_{\mathrm{low-intrinsic}}$ fraction characterizes the _intrinsic_
density of LMXBs with BHs within the lower mass gap, while
$f_{\mathrm{low}}(p_{\mathrm{detect}})$ characterizes the _observed_ density
of LMXBs with a mass-gap BH. Considering
$f_{\mathrm{low}}(p_{\mathrm{detect}})$ and $f_{\mathrm{low-intrinsic}}$
together highlights the influence of transient detection effects on the
observed Galactic LMXB sample.

Following the procedures above, we generate synthetic LMXB samples for each
COSMIC population and MESA grid. We discuss these samples and their
implications for the lower mass gap in Section 3.

##  3 Results

Using rapid binary population synthesis (COSMIC) coupled with detailed
stellar-structure and mass-transfer models (MESA), we investigate the impact
of observational biases on the Galactic LMXB sample. We demonstrate that
transient LMXB selection effects do introduce some biases into the observed
LMXB sample, and that the observed lower mass gap has implications for the
maximum NS birth-mass $M_{\mathrm{NS,birth-max}}$.

To survey how well population synthesis methods match the observed LMXB
sample, and to explore dependencies between the LMXB sample and uncertain
aspects of binary evolution physics, our investigation begins with a suite of
COSMIC rapid binary population synthesis simulations. The COSMIC results are
discussed in Section 3.1. For thirty-two combinations of binary evolution
prescriptions, Galactic LMXB samples from COSMIC host mass-gap BHs, with a
significant fraction of mass-gap BH forming through AIC of a NS. Regardless of
binary evolution prescriptions, the COSMIC populations show minimal to
moderate dependence on the adopted X-ray outburst selection effects and favor
$M_{\mathrm{NS,birth-max}}\lesssim 2M_{\odot}$ (to suppress the formation of
AIC LMXB).

Given the consistency of the COSMIC results, we next consider a three-
dimensional grid of LMXBs evolved using MESA. These results are discussed in
Section 3.2. While the MESA grid samples a narrower range of binary evolution
prescriptions than the COSMIC populations, MESA employs a more detailed
treatment of XRB mass-transfer. For $M_{\mathrm{NS,birth-max}}=3M_{\odot}$,
MESA and COSMIC yield similar BH mass distributions, but if
$M_{\mathrm{NS,birth-max}}\lesssim 2M_{\odot}$ MESA produces a dearth of LMXBs
with mass-gap BHs and reveals a greater dependence on the adopted transition
to RIA in the disk instability model.

###  3.1 COSMIC results

Using COSMIC, we consider proof-of-concept Milky Way populations over a range
of binary evolution prescriptions. In addition to evaluating the ability of
population synthesis methods to replicate the observed LMXB sample, we
leverage the simulation suite to investigate dependencies between the LMXB
sample and uncertain aspects of binary evolution physics. Given the
substantial level of uncertainty in several binary evolution prescriptions and
the order-of-magnitude nature of our study, treating the COSMIC simulations in
aggregate is preferred. Here, we first outline the principal results of the
simulation suite, before investigating dependencies between BH mass-gap
occupancy and different prescriptions of binary evolution physics.

For a given Milky Way COSMIC population, a synthetic observed sample of LMXBs
is extracted from the evolution tracks through a Monte Carlo process (Section
2.6). In Table 2, the $p_{\mathrm{detect}}$-weighted fraction of LMXBs in the
synthetic observed sample hosting a BH in the mass gap
$f_{\mathrm{low}}(p_{\mathrm{detect}})$ is reported for each population.
Detection probabilities are assigned under the assumption that there is no
transition to RIA. For a selection of representative populations, the mass and
period distributions of the synthetic observed LMXBs are presented in Figure
2. Every model predicts $f_{\mathrm{low}}(p_{\mathrm{detect}})$ is greater
than $40\%$ (using a benchmark $M_{\mathrm{NS,birth-max}}=3M_{\odot}$ and
$M_{\mathrm{bound}}=4.5M_{\odot}$).

For LMXBs hosting a BH within the mass gap, the BHs formed through one of two
channels: AIC of a NS or intrinsic formation of a BH (directly post core-
collapse of a star). While both channels are dependent on the maximum NS mass
allowed by the equation of state (because it sets the transition between NS
and BH in the compact object birth mass distribution), the AIC LMXBs also rely
on $M_{\mathrm{NS,birth-max}}$: if $M_{\mathrm{NS,birth-max}}$ lies
significantly below the maximum NS mass, the formation of AIC LMXBs will be
suppressed relative to a population where $M_{\mathrm{NS,birth-max}}$ lies
near the maximum NS mass.

Figure 3:  From population synthesis, each combination of binary evolution
prescription and disk instability model yields an observed LMXB mass
distribution with a substantial fraction of mass-gap BHs (for
$M_{\mathrm{NS,birth-max}}=3M_{\odot}$). _Left_ : The detection probability-
weighted fraction of LMXBs with mass gap BHs,
$f_{\mathrm{low}}(p_{\mathrm{detect}})$, is presented against the unweighted
fraction, $f_{\mathrm{low-intrinsic}}$, for each model (using either no
transition to RIA, a sharp transition, or a smooth transition). For points
that fall above the dashed-black line, detection effects lead to a dearth of
mass-gap BH detections relative to the intrinsic population. _Right_ : The
$p_{\mathrm{discover}}$ detection-weighted fraction of LMXBs with a mass-gap
BH is presented against the $p_{\mathrm{transition}}$-weighted fraction. This
panel breaks-down the two components of the LMXB detection probability: (i)
being bright enough during outburst for discovery $p_{\mathrm{discover}}$, and
(ii) transitioning from outburst to quiescence during the survey
$p_{\mathrm{transition}}$. Relative to $p_{\mathrm{transition}}$, weighting by
$p_{\mathrm{discover}}$ suppresses observation of mass-gap BHs.  Figure 4:
The detection-weighted fraction of LMXBs with a mass gap BH
$f_{\mathrm{low}}(p_{\mathrm{detect}})$ as a function of $M_{\mathrm{NS,birth-
max}}$ for a selection of populations. Fractions for each population are
calculated by averaging over all $5000$ snapshots and the uncertainties are
taken from the $68\%$ uncertainty across the snapshots. Lowering the maximum
NS birth mass $M_{\mathrm{NS,birth-max}}$ reduces the occurrence of LMXB
systems with $M_{\mathrm{BH}}<4.5M_{\odot}$. The magnitude of this trend
depends on the adopted accretion prescription; changing the magnetic braking
prescription moderately effects the fraction of mass-gap BH detections, while
modifying the adopted accretion physics significantly affects
$f_{\mathrm{low}}(p_{\mathrm{detect}})$.

AIC LMXBs account for a significant fraction ($>15\%$ and reaching as high as
$90\%$, for $M_{\mathrm{NS,birth-max}}=3M_{\odot}$ and assuming no transition
to RIA) of the detection-weighted synthetic observed populations. The AIC
LMXBs typically host BHs with $M_{\mathrm{BH}}<5M_{\odot}$. The populations
where the fraction of LMXBs formed via AIC is below $30\%$ all adopted a mass-
accretion efficiency of $50\%$, while the populations where the fraction of
LMXBs formed via AIC is above $90\%$ all adopted the State dependent mass-
accretion prescription. The fraction of LMXBs formed via AIC is lowest for the
population with $\alpha=1.0$, Lifetime limited $q$, Pessimistic CE survival,
magnetic braking following IT03, and $50\%$ accretion efficiency; the fraction
of LMXBs formed via AIC is maximized for the population with $\alpha=5.0$,
Lifetime limited minimum $q$, Pessimistic CE survival, magnetic braking
following IT03, and State dependent accretion efficiency.

AIC LMXBs also dominate ($>20\%$ and reaching as high as $99\%$) the sub-
population of systems with orbital periods shorter than $10$ hours. Consistent
with the total fraction of LMXBs formed via AIC, the fraction of AIC LMXBs
among systems with periods shorter than $10$ hours is reduced by adopting
$50\%$ accretion efficiency and maximized by adopting State dependent
accretion efficiency.

We next consider how the principal results of the COSMIC simulation suite
(i.e., population synthesis fails to replicate the observed LMXB population
and consistently produces a population of AIC LMXBs) depends on the adopted
binary evolution prescriptions, beginning with the details of the disk
instability physics. From population synthesis, the detection of LMXBs hosting
mass-gap BHs is largely insensitive to the adopted transition to RIA. In
Figure 3, $f_{\mathrm{low}}(p_{\mathrm{detect}})$ is presented against
$f_{\mathrm{low-intrinsic}}$ for every combination of population synthesis
prescription and X-ray outburst prescription. Since $f_{\mathrm{low-
intrinsic}}$ characterizes the _intrinsic_ density of BHs within the lower
mass gap, and $f_{\mathrm{low}}(p_{\mathrm{detect}})$ characterizes the
_observed_ density of LMXBs with a mass-gap BH, considering them together
highlights the influence of transient detection effects on the observed
Galactic LMXB sample.

As shown in Figure 3, modeling the X-ray outburst with no RIA transition
yields $f_{\mathrm{low}}(p_{\mathrm{detect}})$ consistent with modeling the
X-ray outburst with a smooth RIA transition. For both of these prescriptions,
the detection-weighted mass distribution is only marginally different from the
unweighted distribution: $f_{\mathrm{low}}(p_{\mathrm{detect}})$ is $\sim
10\%$ higher than $f_{\mathrm{low-intrinsic}}$. Adopting a sharp transition to
RIA lowers $f_{\mathrm{low}}(p_{\mathrm{detect}})$ relative to the other RIA
transitions; however, each population still yields an
$f_{\mathrm{low}}(p_{\mathrm{detect}})$ greater than $40\%$ (using a benchmark
$M_{\mathrm{NS,birth-max}}=3M_{\odot}$ and $M_{\mathrm{bound}}=4.5M_{\odot}$).
Since the outburst duration $T_{\mathrm{outburst}}$ is positively correlated
with orbital period, and the AIC LMXBs typically have orbital periods shorter
than $10$ hours, the sensitivity of $f_{\mathrm{low}}(p_{\mathrm{detect}})$ to
the RIA transition prescription is expected. The sharp transition dramatically
decreases the outburst duration (Figure 1), making the short-period LMXBs,
which tend to be AIC LMXBs, no longer detectable.

Discovering a transient LMXB source with an ASM relies on two independent
probabilities: (i) the probability that the source is bright enough in
outburst to be discovered ($p_{\mathrm{discover}}$), and (ii) the probability
that the source transitions from outburst to quiescence during the ASM’s
observing window ($p_{\mathrm{transition}}$). In Figure 3, the
$p_{\mathrm{discover}}$-weighted fraction of LMXB with
$M_{\mathrm{BH}}<4.5M_{\odot}$ is presented against the
$p_{\mathrm{transition}}$-weighted fraction. Similar to $f_{\mathrm{low-
intrinsic}}$, these fractions are not directly comparable to the observed LMXB
sample; however, they outline the relative effects of $p_{\mathrm{discover}}$
and $p_{\mathrm{transition}}$ on the mass gap. As shown in Figure 3, weighting
by $p_{\mathrm{transition}}$ fills the mass gap relative to
$p_{\mathrm{discover}}$ weighting. Compared to short-period LMXBs, long-period
LMXBs are disfavored by $p_{\mathrm{transition}}$, because the quiescence time
is positively related with orbital period. Since AIC LMXBs tend to have
shorter periods than LMXBs with intrinsic BHs, $p_{\mathrm{transition}}$
favors AIC LMXBs and raises the density of BHs in the mass gap relative to
$p_{\mathrm{discover}}$ weighting.

The prevalence of LMXBs hosting a mass-gap BH in the COSMIC populations is
discrepant with the Galactic LMXB sample. While AIC LMXBs comprise a
substantial fraction of the COSMIC populations for a variety of disk
instability models and binary evolution prescriptions, these systems are
dependent on $M_{\mathrm{NS,birth-max}}$. Lowering $M_{\mathrm{NS,birth-max}}$
suppresses the number of AIC LMXBs, so tuning $M_{\mathrm{NS,birth-max}}$
could bring the population synthesis models into better agreement with
observations. Below, we investigate the dependence between the AIC LMXB
population and $M_{\mathrm{NS,birth-max}}$.

In Figure 4, $f_{\mathrm{low}}(p_{\mathrm{detect}})$ is presented as a
function of $M_{\mathrm{NS,birth-max}}$ for a selection of COSMIC populations.
While $f_{\mathrm{low}}(p_{\mathrm{detect}})$ decreases for each population as
$M_{\mathrm{NS,birth-max}}$ is lowered, the magnitude of the effect depends
strongly on the accretion-limit prescription. This dependence reflects the
higher occurrence rate of AIC LMXBs in populations with State dependent
accretion limits. For the State dependent prescription, compact objects have
conservative mass-transfer, provided the mass-transfer rate is sub-Eddington;
accretors therefore retain a greater fraction of the transferred mass, which
corresponds to more NSs undergoing AIC. The increased accretion rates from the
State dependent prescription also translate to greater mass gain for the
intrinsic BH LMXBs, which in turn reduces
$f_{\mathrm{low}}(p_{\mathrm{detect}})$. Populations with $50\%$ accretion
efficiency have both higher $f_{\mathrm{low}}(p_{\mathrm{detect}})$ for
intrinsic BHs and weaker dependence on $M_{\mathrm{NS,birth-max}}$.

Over the range of uncertain binary evolution physics that we considered, LMXB
populations modeled with COSMIC yield a substantial fraction of LMXB BHs with
$M_{\mathrm{BH}}<4.5M_{\odot}$ (over $40\%$, assuming no transition to RIA)
for $M_{\mathrm{NS,birth-max}}=3M_{\odot}$. The simulations are in tension
with the observed sample of LMXBs, where there is a dearth of BHs with masses
below $\approx 4.5M_{\odot}$ (Bailyn et al., 1998; Özel et al., 2010; Farr et
al., 2011). The fraction of mass-gap BHs in the synthetic Galactic sample is
minimized by adopting a CE efficiency of $\alpha=5$, $q>0.01$, magnetic
braking following IT03, State dependent accretion efficiency, pessimistic CE
survival, no transition to RIA, and $M_{\mathrm{NS,birth-max}}=1.5M_{\odot}$;
for this population, $28\%$ of LMXBs host a BH in the mass gap. While lowering
$M_{\mathrm{NS,birth-max}}$ lowers the density of LMXB systems in the mass
gap, $f_{\mathrm{low}}(p_{\mathrm{detect}})$ remains greater than $25\%$ for
every COSMIC population and is substantially higher for most, assuming no
transition to RIA.

Table 3: COSMIC and MESA grid comparison

| Magnetic braking | Accretion limit | RIA transition |  $f_{\mathrm{low}}(M_{\mathrm{NS,birth}}<1.5M_{\odot})$aaThe fraction $f_{\mathrm{low}}(p_{\mathrm{detect}})$ is defined relative to $M_{\mathrm{bound}}=4.5M_{\odot}$.  | $f_{\mathrm{low}}(M_{\mathrm{NS,birth}}<3.0M_{\odot})$
---|---|---|---|---|---
COSMIC | H02 |  $50\%$ efficiency | No transition | 0.345 | 0.406
| H02 |  $50\%$ efficiency |  $\nu=0.05$, sharp | 0.422 | 0.492
| H02 |  $50\%$ efficiency |  $\nu=0.05$, smooth | 0.347 | 0.407
| IT03 |  $50\%$ efficiency | No transition | 0.345 | 0.406
| IT03 |  $50\%$ efficiency |  $\nu=0.05$, sharp | 0.418 | 0.490
| IT03 |  $50\%$ efficiency |  $\nu=0.05$, smooth | 0.344 | 0.405
| H02 | State dependent | No transition | 0.354 | 0.437
| H02 | State dependent |  $\nu=0.05$, sharp | 0.439 | 0.516
| H02 | State dependent |  $\nu=0.05$, smooth | 0.351 | 0.435
| IT03 | State dependent | No transition | 0.352 | 0.435
| IT03 | State dependent |  $\nu=0.05$, sharp | 0.439 | 0.515
| IT03 | State dependent |  $\nu=0.05$, smooth | 0.353 | 0.437
MESA |  |  | No transition | 0.105 | 0.428
|  |  |  $\nu=0.05$, sharp | 0.242 | 0.507
|  |  |  $\nu=0.05$, smooth | 0.105 | 0.429

Note. — Relative to rapid population synthesis, LMXBs evolved using MESA
under-populate the mass gap if $M_{\mathrm{NS,birth-max}}\lesssim 2M_{\odot}$,
and show greater dependence on the adopted RIA transition. The detection-
weighted fraction of LMXBs with a BH in the mass gap
$f_{\mathrm{low}}(p_{\mathrm{detect}})$ is presented for a three-dimensional
grid of binaries evolved with MESA, as well as for grids evolved using COSMIC;
the grids evolved with COSMIC are initialized at the same initial conditions
as the grid evolved using MESA. Values of the detection-weighted fraction
$f_{\mathrm{low}}(p_{\mathrm{detect}})$ are reported for each combination of
$M_{\mathrm{NS,birth-max}}=1.5M_{\odot}$ or $3.0M_{\odot}$ and smooth, sharp,
or non-existent transition to RIA. To marginalize over uncertain aspects of
binary evolution physics, different combinations of binary evolution
prescriptions are adopted for the COSMIC grids, including the magnetic braking
implementation (following either H02 or IT03) and the efficiency of accretion
during RLO (either $50\%$ efficiency or State dependent efficiency, where
compact objects have conservative mass accretion up to the Eddington limit).

Figure 5:  _Top_ : The detection-weighted fraction of LMXBs with a mass gap BH
$f_{\mathrm{low}}(p_{\mathrm{detect}})$ as a function of $M_{\mathrm{NS,birth-
max}}$. Each color corresponds to a different RIA transition. For each panel,
the results for no transition to RIA overlap with the results for a smooth
transition to RIA. _Bottom_ : The average number of detections of AIC LMXBs
(over $N_{\mathrm{trials}}=5000$ snapshots of the evolutionary tracks) with
$M_{\mathrm{BH}}<4.5M_{\odot}$ (as a solid line); the dashed horizontal line
is the average number of detections of LMXB that formed with a BH and have
$M_{\mathrm{BH}}<4.5M_{\odot}$; the dotted horizontal line is the average
number of detections of LMXB (from any formation channel) with
$M_{\mathrm{BH}}>4.5M_{\odot}$. The _left column_ presents the results from
the MESA grid, while the _right column_ presents the results from the COSMIC
grid. While MESA and COSMIC predict similar LMXB BH mass distributions when
$M_{\mathrm{NS,birth-max}}=3M_{\odot}$, MESA produces a dearth of LMXBs with
mass gap BHs relative to COSMIC if $M_{\mathrm{NS,birth-max}}\lesssim
2M_{\odot}$.  Figure 6:  The LMXB orbital period distribution for the MESA
grid is presented in the _left_ column, while the orbital period distributions
of the COSMIC grids (which are generated with matching initial conditions) are
presented on the _right_. _Top_ : Monte Carlo sampled orbital period
distributions are presented for three sub-populations of LMXBs: AIC LMXBs
(purple), LMXBs with intrinsic BHs in the mass gap (gold), and LMXBs with
intrinsic BHs above the mass gap (gray). _Bottom_ : Monte Carlo sampled and
detection-weighted orbital period distributions for the entire population of
LMXBs are presented for either no transition to RIA (blue) or a sharp
transition (teal); here, applying a smooth transition to RIA is nearly
indistinguishable from applying no transition. LMXBs evolved using MESA
typically occupy shorter orbital periods than LMXBs evolved with COSMIC.
Figure 7: The differences in the orbital period and mass distributions between
COSMIC and MESA both originate from the treatment of RLO mass accretion.
Relative to COSMIC, evolving LMXBs with MESA leads to greater mass accretion
by the BHs, and by extension, shorter orbital periods. Here the evolution
tracks of a single LMXB evolved with MESA (dashed lines) are presented against
the evolution tracks when COSMIC is used (solid lines). The binary is
initialized with $M_{\mathrm{BH}}=2.5M_{\odot}$,
$M_{\mathrm{donor}}=1.16M_{\odot}$, and $P=0.413$ days. For the COSMIC
evolution tracks, the magnetic braking formulation of H02 and the State
dependent accretion efficiency were adopted.

###  3.2  MESA results

Our proof-of-concept rapid population synthesis Milky Way populations
consistently populate the lower mass gap. These simulations also reveal a
dependence between $M_{\mathrm{NS,birth-max}}$ and the existence of the mass
gap. Informed by the COSMIC simulations, here we investigate a grid of binary
MESA models. We adopt a three-dimensional grid of binaries (in compact object
mass, donor mass, and orbital period space), consisting of a compact object
and a hydrogen-rich main-sequence star from Fragos et al. (2023). A series of
COSMIC populations are initialized with the same initial conditions as the
MESA grid, to facilitate comparison.

Using the Monte Carlo sampling procedure described in Section 2.6,
$f_{\mathrm{low}}(p_{\mathrm{detect}})$ is reported for the MESA grid and each
COSMIC grid in Table 3, for a range of accretion disk treatments and
$M_{\mathrm{NS,birth-max}}$. For $M_{\mathrm{NS,birth-max}}=3M_{\odot}$, the
MESA and COSMIC grids yield generally consistent
$f_{\mathrm{low}}(p_{\mathrm{detect}})$. However, when $M_{\mathrm{NS,birth-
max}}$ is lowered, MESA predicts a consistently lower
$f_{\mathrm{low}}(p_{\mathrm{detect}})$ than COSMIC, and the dependence on the
assumed transition to RIA is exacerbated. When $M_{\mathrm{NS,birth-
max}}\lesssim 2M_{\odot}$, $f_{\mathrm{low}}(p_{\mathrm{detect}})$ for a LMXB
population modeled with MESA is approximately half that of a LMXB population
modeled with COSMIC.

In Figure 5, $f_{\mathrm{low}}(p_{\mathrm{detect}})$ is presented as a
function of $M_{\mathrm{NS,birth-max}}$ for the MESA grid and a representative
COSMIC grid; for a more comprehensive view, the average number of LMXB
detections (separated by formation channel and BH mass) is presented as well.
For $M_{\mathrm{NS,birth-max}}\lesssim 2M_{\odot}$, the MESA grid yields
$f_{\mathrm{low}}(p_{\mathrm{detect}})$ a factor of two lower than the COSMIC
grids. This disparity is enhanced when either a smooth or nonexistent
transition to RIA is adopted. The sharp RIA transition actually reduces the
number of detected AIC LMXBs for the MESA grid. However, the behavior of
$f_{\mathrm{low}}(p_{\mathrm{detect}})$ is dominated by the increased
detection of LMXBs with $M_{\mathrm{BH}}>4.5M_{\odot}$ when either a smooth or
nonexistent transition to RIA is adopted for the MESA grid and the less
dramatic increase for the COSMIC grid (i.e., a smooth or nonexistent
transition to RIA leads to slightly more mass-gap BH detections but a far
greater increase in the number of higher mass LMXB detections).

Similar to the BH mass distributions, the orbital period distributions from
MESA and COSMIC also diverge. In Figure 6, orbital period distributions are
presented for three sub-populations of LMXBs (AIC LMXBs, LMXBs with intrinsic
BH in the mass gap, and LMXBs with intrinsic BH above the mass gap) from the
MESA grid and each COSMIC grid. To highlight the influence of
$p_{\mathrm{detect}}$, Figure 6 includes both the underlying orbital period
distributions (i.e., only Monte Carlo sampled) and the synthetic observed
orbital period distributions (i.e., Monte Carlo sampled and detection
weighted).

While the observed orbital period distributions of the COSMIC grids are
generally consistent with MESA, the underlying distributions show significant
differences. Relative to the COSMIC grids, the underlying orbital period
distribution of AIC LMXBs is skewed towards shorter orbital periods in the
MESA grid; without $p_{\mathrm{detect}}$ weighting, $80\%$ of AIC LMXBs in the
MESA grid have $P<10$ hour, but $<35\%$ of AIC LMXBs in the COSMIC grids have
$P<10$ hour. For LMXBs with intrinsic BHs and $P<10$ hour, the density of
LMXBs with mass-gap BHs is lower in the MESA grid, relative to the COSMIC
grids; in the MESA grid, $f_{\mathrm{low-intrinsic}}=7\%$ for LMXBs with
intrinsic BHs and $P<10$ hour, while the fraction is $\sim 30\%$ in the COSMIC
grids. For $P>10$ hour, this trend continues, but the magnitude is lessened;
in the MESA grid, $f_{\mathrm{low-intrinsic}}=15\%$ for LMXBs with intrinsic
BH and $P>10$ hour, while the fraction is $\sim 20\%$ for the COSMIC grids.

The differences in the orbital period and mass distributions between COSMIC
and MESA are both consequences of how RLO mass transfer is treated in each
code. COSMIC treats RLO mass transfer following Hurley et al. (2002), in which
stars are (nonphysically) approximated as in thermal equilibrium throughout
mass loss. For main-sequence donors in the MESA models, mass is removed for
stars overfilling their Roche lobes at the beginning of each time-step, such
that the star’s radius remains within its Roche lobe, while for giant stars,
the methods of Kolb & Ritter (1990) are adopted; stars’ radii are allowed to
extend beyond their Roche lobes, and mass transfer is self-consistently
calculated based on the local fluid conditions (Fragos et al., 2023).
Differences in mass-transfer rate naturally affect a LMXB’s BH mass and
orbital evolution.

Within the first $100~{}\mathrm{Myr}$ of evolution, every compact object in
the grid accretes more mass when MESA is used than when COSMIC used.
Similarly, every donor star in the grid loses more mass when MESA is used than
when COSMIC is used. Of LMXBs in the grid that do not merge by
$100~{}\mathrm{Myr}$, the compact objects of LMXBs evolved using MESA accrete
on average $0.7M_{\odot}$ more mass by $100~{}\mathrm{Myr}$ than compact
objects in LMXBs evolved using COSMIC; weighting by detection probability
(assuming no transition to RIA), the compact objects of LMXBs evolved using
MESA accrete on average $0.4M_{\odot}$ more mass by $100~{}\mathrm{Myr}$ than
the compact objects of LMXBs evolved using COSMIC. The greater mass accretion
by compact objects in MESA results in both more massive compact objects and
shorter orbital periods.

To highlight these differences, the evolution tracks of a single LMXB evolved
using MESA are presented against the evolution tracks from COSMIC (for a LMXB
with the same initial conditions) in Figure 7. The majority of both the
orbital period evolution and mass gain for the compact object occurs in the
first few hundred million years in the MESA grid. During this period, the
mass-loss rate from the donor star is larger in MESA than in COSMIC by at
least two orders of magnitude. Since the donor mass loss rate is still sub-
Eddington, the mass exchange is conservative for both codes. In contrast, the
timescales for orbital period evolution and mass gain in the COSMIC grid
extend to several billion years. This leads to more dramatic orbital period
shrinking and more accretion onto the compact object and thus results in the
overall higher-mass BHs and shorter periods in the MESA grid when compared to
COSMIC.

In summary, a COSMIC population study of uncertain aspects of binary evolution
and disk instability models consistently found a substantial fraction of LMXBs
hosting BHs with $M_{\mathrm{BH}}<4.5M_{\odot}$ (over $40\%$) for
$M_{\mathrm{NS,birth-max}}=3M_{\odot}$. After lowering $M_{\mathrm{NS,birth-
max}}$ to $1.5M_{\odot}$, $f_{\mathrm{low}}(p_{\mathrm{detect}})>25\%$ for
every COSMIC population. For $M_{\mathrm{NS,birth-max}}=3M_{\odot}$, MESA and
COSMIC yield similar BH mass distributions and both fill the lower mass gap.
However, if $M_{\mathrm{NS,birth-max}}$ is lowered to below $2M_{\odot}$, the
LMXBs produced with MESA have at least a factor of two fewer BHs with
$M_{\mathrm{BH}}<4.5M_{\odot}$. This difference is exacerbated when a smooth
or nonexistent transition to RIA is adopted. The results highlight the
importance of accurately modeling mass transfer in understanding the formation
of compact objects.

##  4 Discussion

While Galactic LMXBs have previously been used to constrain the underlying
stellar-mass BH distribution (Bailyn et al., 1998; Özel et al., 2010; Farr et
al., 2011), the impact of detection biases on the observed LMXB sample has
remained unclear (Özel et al., 2010; Kreidberg et al., 2012; Jonker et al.,
2021). Here, we presented a quantitative study of selection biases, motivated
by the fact that BH mass measurement is possible only for _transient_ XRB
systems. We examined whether the requirement for transient behavior can hide
an underlying population of low-mass BHs and create an _apparent_ mass gap in
the observed sample.

We followed a forward-modeling approach. Using both rapid binary population
synthesis (COSMIC) and detailed stellar evolution models (MESA), we generated
initial synthetic binary systems, evolved the systems forward in time, and
inferred the detectable populations using the disk instability model of LMXB
transient behavior. To reflect the criteria for LMXB mass inference, synthetic
samples were weighted by the probability a RXTE-like ASM (Levine et al., 1996)
would detect the source as transient. From these results, we demonstrated that
transient LMXB selection effects do introduce some biases into the observed
sample, but this is not sufficient to completely hide mass-gap BHs.
Additionally, observations of the lower mass gap have implications for
$M_{\mathrm{NS,birth-max}}$, as this impacts the number of BHs (or lack
thereof) formed through AIC. We conclude that while selection effects are
important in understanding XRB observations, the lack of LMXB BH seen in the
lower mass gap places constraints upon viable models of binary stellar
evolution.

Using COSMIC, we generated Milky Way binary populations for an array of binary
evolution prescriptions and initial conditions. While current rapid population
synthesis methods likely cannot reproduce all properties of the observed LMXB
sample, they enable exploratory studies to uncover any significant
dependencies between the synthetic observed LMXB sample and the binary
evolution prescriptions. In concert with the suite of population synthesis
simulations, we considered a grid of MESA binary models. We adopted a three-
dimensional grid of binaries consisting of a compact object and a hydrogen-
rich main-sequence star from Fragos et al. (2023). While the MESA grid does
not consider the full range of binary evolution prescriptions investigated
with COSMIC, here we explored an improved physical treatment of XRB mass-
transfer.

Our key results are:

  1. 1.

Our rapid population synthesis results cannot replicate the observed mass
distribution of Galactic LMXBs. The COSMIC results indicate that detection of
LMXBs containing mass-gap BHs is a robust prediction of the Delayed explosion
mechanism. For thirty-two combinations of binary evolution prescriptions and
three disk instability models (smooth, sharp or nonexistent transition to RIA)
over $40\%$ of the LMXBs in the synthetic Galactic samples have
$M_{\mathrm{BH}}<4.5M_{\odot}$ (with $M_{\mathrm{NS,birth-max}}=3M_{\odot}$).

  2. 2.

LMXBs with a mass-gap BH form through one of two channels: AIC of a NS or
directly post core-collapse of a star. Unlike LMXBs with intrinsic BHs, the
formation of AIC LMXBs relies upon $M_{\mathrm{NS,birth-max}}$ (i.e., lowering
$M_{\mathrm{NS,birth-max}}$ suppresses the formation of AIC LMXBs). Since
$M_{\mathrm{BH}}\lesssim 5M_{\odot}$ for the AIC LMXBs, the existence of a
lower mass gap in the Galactic LMXB sample can potentially constrain
$M_{\mathrm{NS,birth-max}}$.

  3. 3.

While rapid population synthesis and MESA models are consistent for high
values of the maximum NS birth mass (e.g., $M_{\mathrm{NS,birth-
max}}=3M_{\odot}$), LMXBs evolved with MESA have a lower fraction of BHs in
the mass gap relative to the COSMIC models, if $M_{\mathrm{NS,birth-
max}}\lesssim 2M_{\odot}$. Among the synthetic LMXB samples, the MESA models
produce a factor of two fewer LMXBs with a mass-gap BH than the analogous
COSMIC systems. This trend is exacerbated when a smooth or nonexistent
transition to RIA is adopted.

  4. 4.

Relative to our rapid population synthesis models, BHs in LMXBs evolved using
MESA generally undergo greater mass accretion. This results in both more
massive BHs and shorter period LMXBs in the MESA models.

Drawing population-level conclusions from the MESA grid alone presents
considerable challenges. Unlike the COSMIC Milky Way models, the initial
conditions for the MESA grid are not drawn from informed distributions, and
the MESA grid only considers the compact object–main sequence phase of
evolution. While we find that the COSMIC simulated LMXB samples host
significantly more BHs in the mass gap than samples evolved using MESA, we
cannot predict the occurrence rate of mass-gap BH for the Milky Way using the
MESA grid.

To suppress the formation of AIC LMXBs, both the COSMIC Milky Way populations
and the MESA grid favor $M_{\mathrm{NS,birth-max}}\lesssim 2M_{\odot}$. This
constraint on $M_{\mathrm{NS,birth-max}}$ is consistent with the observed NS
mass distribution. Masses have been inferred for $>70$ NSs via pulsar timing,
eclipsing X-ray binaries, and optical radial velocity measurements of a NS’s
companion (Alsing et al., 2018; Lattimer, 2019). From these measurements, the
NS mass distribution displays a strong peak near $1.4M_{\odot}$ with a few
objects heavier than $\sim 1.6M_{\odot}$ and a maximum around $\sim
2.3M_{\odot}$ (Shao et al., 2020). However, the most massive NSs are in
binaries with low-mass companions (typically white dwarfs) and likely
underwent accretion growth. With the exceptions of GW190814 and
GW200210_092254, which are consistent with containing a compact object lying
in the mass gap (Abbott et al., 2020, 2021c, 2021b), there are few
gravitational-wave sources with significant support for a mass-gap component,
and observations are consistent with a dip in the NS population above $\sim
2.1M_{\odot}$ (Abbott et al., 2021a). The preference for $M_{\mathrm{NS,birth-
max}}\lesssim 2M_{\odot}$ in both the COSMIC and MESA simulations is
approximately consistent with these observed samples.

Over a wide range of binary evolution prescriptions, $M_{\mathrm{NS,birth-
max}}$, and disk instability models, COSMIC Milky Way populations fill the
mass gap. In isolation, the robustness of this result (and its contradiction
with observation) would favor a supernova explosion mechanism that inherently
creates a mass gap (e.g., the Rapid prescription of Fryer et al., 2012).
However, gravitational-wave observations still favor mechanisms that (at least
partially) populate the gap (e.g., the Delayed prescription of Fryer et al.,
2012). Therefore, there is still physics regarding the formation of mass-gap
BHs (either from AIC of a NS or core collapse) that remains to be understood.
Given the overproduction of mass-gap BHs with COSMIC relative to MESA, greater
study of LMXB formation and observably using detailed stellar evolution models
(e.g., POSYDON; Fragos et al., 2023) including up-to-date supernova and
stellar-collapse prescriptions (e.g., Ertl et al., 2020; Fryer et al., 2022)
is vital to untangling the mystery of the lower mass gap. Additionally, the
presence of observational selection biases should motivate the building of
more sensitive ASMs than RXTE, such as those proposed for missions like eXTP
and STROBE-X (Zhang et al., 2017; Ray et al., 2019), as well as the use of
sufficient optical and infrared resources to make mass estimates for the
fainter X-ray transients which will tend to have fainter optical counterparts.

##  5 Conclusions

For a supernova engine that fills the lower BH mass gap, transient LMXB
selection effects do introduce significant biases into the observed LMXB
sample. However, unless there are further (unaccounted for) observational
biases against finding LMXBs with mass-gap BHs, population synthesis models
fail to reproduce this aspect of the observed LMXB population. This result is
robust against variations of uncertain aspects of binary evolution physics,
e.g., CE efficiency, CE survivability, minimum ZAMS mass ratio, magnetic
braking, and accretion efficiency. This points to the need for additional
physics not currently included in our COSMIC and MESA simulations, such as a
supernova mechanism that suppresses formation of mass-gap objects.

Regardless of whether the low-mass BHs form from core collapse, the results of
our COSMIC and MESA models lead to the robust implication that the NS birth
masses must be suppressed above $\sim 2M_{\odot}$. Otherwise the mass gap
would be filled by low-mass BHs formed through NS AIC. This constraint on the
maximum NS birth mass is independent of whether the supernova engine forms
mass-gap BHs. This result alone motivates a reexamination of the physics
included in the COSMIC and MESA simulations, such as a supernova engine that
limits the maximum birth mass of NS.

## Acknowledgments

The authors thank Michael Zevin and Simon Stevenson for useful discussions,
and the referee for useful suggestions. JS and IK were supported as CIERA REU
students by the National Science Foundation (NSF) under Grant Number 1757792.
Any opinions, findings, and conclusions or recommendations expressed in this
material are those of the author(s) and do not necessarily reflect the views
of the NSF. JS was also partially supported by CIFAR (through VK’s Senior
Fellowship). VK was partially supported through a CIFAR Senior Fellowship, a
Guggenheim Fellowship, and the Gordon and Betty Moore Foundation (grant award
GBMF8477). CPLB was supported by the CIERA Board of Visitors Research
Professorship. JJA was supported by Northwestern University through a CIERA
Postdoctoral Fellowship. KR was supported by the Gordon and Betty Moore
Foundation (PI Kalogera, grant award GBMF8477) and the Riedel Family Graduate
Fellowship in CIERA. AD, PS, MS were supported by the Gordon and Betty Moore
Foundation (grant award GBMF8477). SBB, TF, KK, DM, and ZX were supported from
the Swiss National Science Foundation Professorship Grant (PP00P2_176868; PI
Fragos). KK acknowledges support from the Federal Commission for Scholarships
for Foreign Students for the Swiss Government Excellence Scholarship (ESKAS
No. 2021.0277). ZX acknowledges support from the Chinese Scholarship Council
(CSC). EZ acknowledges funding support from the European Research Council
(ERC) under the European Union’s Horizon 2020 research and innovation program
(Grant agreement No. 772086). The Flatiron Institute is funded by the Simons
Foundation. This research was supported in part through the computational
resources and staff contributions provided for the Quest high performance
computing facility at Northwestern University which is jointly supported by
the Office of the Provost, the Office for Research, and Northwestern
University Information Technology.

The data behind Table 2 is openly available from the Zenodo repository at
10.5281/zenodo.8155601.

## References

  * Abbott et al. (2021a) Abbott, R., Abbott, T. D., Acernese, F., et al. 2021a, arXiv e-prints, arXiv:2111.03634
  * Abbott et al. (2020) Abbott, R., Abbott, T. D., Abraham, S., et al. 2020, ApJ, 896, L44
  * Abbott et al. (2021b) Abbott, R., Abbott, T. D., Acernese, F., et al. 2021b, arXiv e-prints, arXiv:2111.03606
  * Abbott et al. (2021c) —. 2021c, arXiv e-prints, arXiv:2108.01045
  * Ai et al. (2020) Ai, S., Gao, H., & Zhang, B. 2020, ApJ, 893, 146
  * Alsing et al. (2018) Alsing, J., Silva, H. O., & Berti, E. 2018, MNRAS, 478, 1377
  * Atri et al. (2019) Atri, P., Miller-Jones, J. C. A., Bahramian, A., et al. 2019, MNRAS, 489, 3116
  * Bailyn et al. (1998) Bailyn, C. D., Jain, R. K., Coppi, P., & Orosz, J. A. 1998, ApJ, 499, 367
  * Belczynski et al. (2008) Belczynski, K., Kalogera, V., Rasio, F. A., et al. 2008, ApJS, 174, 223
  * Belczynski et al. (2012) Belczynski, K., Wiktorowicz, G., Fryer, C. L., Holz, D. E., & Kalogera, V. 2012, ApJ, 757, 91
  * Bolton (1972) Bolton, C. T. 1972, Nature, 235, 271
  * Bradt et al. (1993) Bradt, H. V., Rothschild, R. E., & Swank, J. H. 1993, A&AS, 97, 355
  * Breivik et al. (2020) Breivik, K., Coughlin, S., Zevin, M., et al. 2020, ApJ, 898, 71
  * Cannizzo et al. (1982) Cannizzo, J. K., Ghosh, P., & Wheeler, J. C. 1982, ApJ, 260, L83
  * Cannizzo et al. (1988) Cannizzo, J. K., Shafter, A. W., & Wheeler, J. C. 1988, ApJ, 333, 227
  * Corral-Santana et al. (2016) Corral-Santana, J. M., Casares, J., Muñoz-Darias, T., et al. 2016, A&A, 587, A61
  * de Kool (1990) de Kool, M. 1990, ApJ, 358, 189
  * Degenaar et al. (2018) Degenaar, N., Ballantyne, D. R., Belloni, T., et al. 2018, Space Sci. Rev., 214, 15
  * Dubus et al. (2001) Dubus, G., Hameury, J. M., & Lasota, J. P. 2001, A&A, 373, 251
  * Dubus et al. (1999) Dubus, G., Lasota, J.-P., Hameury, J.-M., & Charles, P. 1999, MNRAS, 303, 139
  * El-Badry et al. (2022) El-Badry, K., Seeburger, R., Jayasinghe, T., et al. 2022, MNRAS, 512, 5620
  * Ertl et al. (2020) Ertl, T., Woosley, S. E., Sukhbold, T., & Janka, H. T. 2020, ApJ, 890, 51
  * Farr et al. (2011) Farr, W. M., Sravan, N., Cantrell, A., et al. 2011, ApJ, 741, 103
  * Fragos et al. (2019) Fragos, T., Andrews, J. J., Ramirez-Ruiz, E., et al. 2019, ApJ, 883, L45
  * Fragos et al. (2023) Fragos, T., Andrews, J. J., Bavera, S. S., et al. 2023, ApJS, 264, 45
  * Frank et al. (1992) Frank, J., King, A., & Raine, D. 1992, Accretion power in astrophysics., Vol. 21
  * Fryer et al. (2012) Fryer, C. L., Belczynski, K., Wiktorowicz, G., et al. 2012, ApJ, 749, 91
  * Fryer et al. (2022) Fryer, C. L., Olejak, A., & Belczynski, K. 2022, ApJ, 931, 94
  * Gallegos-Garcia et al. (2021) Gallegos-Garcia, M., Berry, C. P. L., Marchant, P., & Kalogera, V. 2021, ApJ, 922, 110
  * Giacobbo & Mapelli (2018) Giacobbo, N., & Mapelli, M. 2018, MNRAS, 480, 2011
  * Giacobbo & Mapelli (2019) —. 2019, MNRAS, 482, 2234
  * Hobbs et al. (2005) Hobbs, G., Lorimer, D. R., Lyne, A. G., & Kramer, M. 2005, MNRAS, 360, 974
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90
  * Hurley et al. (2002) Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897
  * Ivanova (2006) Ivanova, N. 2006, ApJ, 653, L137
  * Ivanova & Taam (2003) Ivanova, N., & Taam, R. E. 2003, ApJ, 599, 516
  * Jayasinghe et al. (2021) Jayasinghe, T., Stanek, K. Z., Thompson, T. A., et al. 2021, MNRAS, 504, 2577
  * Jonker et al. (2021) Jonker, P. G., Kaur, K., Stone, N., & Torres, M. A. P. 2021, ApJ, 921, 131
  * Kalogera & Baym (1996) Kalogera, V., & Baym, G. 1996, ApJ, 470, L61
  * King et al. (1996) King, A. R., Kolb, U., & Burderi, L. 1996, ApJ, 464, L127
  * King et al. (1997) King, A. R., Kolb, U., & Szuszkiewicz, E. 1997, ApJ, 488, 89
  * King & Ritter (1998) King, A. R., & Ritter, H. 1998, MNRAS, 293, L42
  * Knevitt et al. (2014) Knevitt, G., Wynn, G. A., Vaughan, S., & Watson, M. G. 2014, MNRAS, 437, 3087
  * Knigge et al. (2011) Knigge, C., Baraffe, I., & Patterson, J. 2011, ApJS, 194, 28
  * Kolb & Ritter (1990) Kolb, U., & Ritter, H. 1990, A&A, 236, 385
  * Kreidberg et al. (2012) Kreidberg, L., Bailyn, C. D., Farr, W. M., & Kalogera, V. 2012, ApJ, 757, 36
  * Lasota (2001) Lasota, J.-P. 2001, New A Rev., 45, 449
  * Lasota (2016) —. 2016, Astrophysics and Space Science Library, Vol. 440, Black Hole Accretion Discs, ed. C. Bambi, 1
  * Lasota et al. (2008) Lasota, J. P., Dubus, G., & Kruk, K. 2008, A&A, 486, 523
  * Lattimer (2019) Lattimer, J. M. 2019, Universe, 5, 159
  * Levine et al. (1996) Levine, A. M., Bradt, H., Cui, W., et al. 1996, ApJ, 469, L33
  * Liu et al. (2021) Liu, T., Wei, Y.-F., Xue, L., & Sun, M.-Y. 2021, ApJ, 908, 106
  * Mandel & Müller (2020) Mandel, I., & Müller, B. 2020, MNRAS, 499, 3214
  * Margalit & Metzger (2017) Margalit, B., & Metzger, B. D. 2017, ApJ, 850, L19
  * McMillan (2011) McMillan, P. J. 2011, MNRAS, 414, 2446
  * Müller & Serot (1996) Müller, H., & Serot, B. D. 1996, Nucl. Phys. A, 606, 508
  * Özel & Freire (2016) Özel, F., & Freire, P. 2016, ARA&A, 54, 401
  * Özel et al. (2010) Özel, F., Psaltis, D., Narayan, R., & McClintock, J. E. 2010, ApJ, 725, 1918
  * Patton et al. (2022) Patton, R. A., Sukhbold, T., & Eldridge, J. J. 2022, MNRAS, 511, 903
  * Paxton et al. (2011) Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3
  * Paxton et al. (2013) Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4
  * Paxton et al. (2015) Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15
  * Paxton et al. (2018) Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34
  * Paxton et al. (2019) Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10
  * Price-Whelan et al. (2018) Price-Whelan, A. M., Sipőcz, B. M., Günther, H. M., et al. 2018, AJ, 156, 123
  * Raaijmakers et al. (2021) Raaijmakers, G., Greif, S. K., Hebeler, K., et al. 2021, ApJ, 918, L29
  * Ray et al. (2019) Ray, P., Arzoumanian, Z., Ballantyne, D., et al. 2019, in Bulletin of the American Astronomical Society, Vol. 51, 231
  * Remillard & McClintock (2006) Remillard, R. A., & McClintock, J. E. 2006, ARA&A, 44, 49
  * Rhoades & Ruffini (1974) Rhoades, C. E., & Ruffini, R. 1974, Phys. Rev. Lett., 32, 324
  * Robin et al. (2003) Robin, A. C., Reylé, C., Derrière, S., & Picaud, S. 2003, A&A, 409, 523
  * Robitaille et al. (2013) Robitaille, T. P., Tollerud, E. J., Greenfield, P., et al. 2013, A&A, 558, A33
  * Sana et al. (2012) Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Science, 337, 444
  * Santoliquido et al. (2021) Santoliquido, F., Mapelli, M., Giacobbo, N., Bouffanais, Y., & Artale, M. C. 2021, MNRAS, 502, 4877
  * Shao et al. (2020) Shao, D.-S., Tang, S.-P., Jiang, J.-L., & Fan, Y.-Z. 2020, Phys. Rev. D, 102, 063006
  * Shenar et al. (2022) Shenar, T., Sana, H., Mahy, L., et al. 2022, Nature Astronomy, 6, 1085
  * The pandas development team (2020) The pandas development team. 2020, pandas-dev/pandas: Pandas, v1.1.1, Zenodo, doi:10.5281/zenodo.3509134.  https://doi.org/10.5281/zenodo.3509134
  * Thompson et al. (2019) Thompson, T. A., Kochanek, C. S., Stanek, K. Z., et al. 2019, Science, 366, 637
  * van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science and Engineering, 13, 22
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261
  * Webbink (1984) Webbink, R. F. 1984, ApJ, 277, 355
  * Wen et al. (2006) Wen, L., Levine, A. M., Corbet, R. H. D., & Bradt, H. V. 2006, ApJS, 163, 372
  * Wes McKinney (2010) Wes McKinney. 2010, in Proceedings of the 9th Python in Science Conference, ed. Stéfan van der Walt & Jarrod Millman, 56–61
  * Zevin et al. (2020) Zevin, M., Spera, M., Berry, C. P. L., & Kalogera, V. 2020, ApJ, 899, L1
  * Zevin et al. (2021) Zevin, M., Bavera, S. S., Berry, C. P. L., et al. 2021, ApJ, 910, 152
  * Zhang et al. (2017) Zhang, S. N., Feroci, M., Santangelo, A., et al. 2017, in Space Telescopes and Instrumentation 2016: Ultraviolet to Gamma Ray, ed. J.-W. A. den Herder, T. Takahashi, & M. Bautz, Vol. 9905, International Society for Optics and Photonics (SPIE), 99051Q.  https://doi.org/10.1117/12.2232034