The Mass Distribution of Neutron Stars in Gravitational-wave Binaries

The discovery of two neutron star-black hole coalescences by LIGO and Virgo brings the total number of likely neutron stars observed in gravitational waves to six. We perform the first inference of the mass distribution of this extragalactic population of neutron stars. In contrast to the bimodal Galactic population detected primarily as radio pulsars, the masses of neutron stars in gravitational-wave binaries are thus far consistent with a uniform distribution, with a greater prevalence of high-mass neutron stars. The maximum mass in the gravitational-wave population agrees with that inferred from the neutron stars in our Galaxy and with expectations from dense matter.


INTRODUCTION
To date, the vast majority of neutron star (NS) masses have been measured through the binary interactions of Galactic radio pulsars (Özel & Freire 2016). Beginning with Thorsett & Chakrabarty (1999), population-level analyses have been performed on these Galactic pulsars, revealing a mass distribution with a double-peaked structure (Antoniadis et al. 2016;Alsing et al. 2018;Shao et al. 2020). The double NSs among these Galactic observations, however, have only been found with both component masses near 1.35 M (Özel et al. 2012;Kiziltan et al. 2013).
The first NS binary observed in gravitational waves (GWs), GW170817 (Abbott et al. 2017), fell neatly into the mass range anticipated from these radio observations. On the other hand, GW190425 (Abbott et al. 2020a), the second low-mass GW binary discovered, illustrated that the compact object population detected with GWs need not be the same as that observed in our Galaxy: its total mass is a 5σ outlier relative to known Galactic double NSs. With the recent neutron star-black hole (NSBH) discoveries GW200105 and GW200115 (Abbott et al. 2021a), LIGO (LIGO Scientific Collaboration et al. 2015) and Virgo (Acernese et al. 2015) have now recorded four confirmed signals whose progenitor likely contains at least one NS. In this Letter, we consider these to be six extragalactic NSs known through GW astronomy, and investigate the collective properties of this incipient NS population. We follow the hierarchical Bayesian population inference framework of Mandel (2010) and Mandel et al. (2019), which accounts for statistical uncertainty in the measured parameters of each GW signal and for selection effects in the observed population. This kind of population inference has a history of applications inside (e.g. Abbott et al. 2019;Roulet et al. 2020;Abbott et al. 2021b) and outside (e.g. Loredo 2004;Hogg et al. 2010;Foreman-Mackey et al. 2014) of GW astronomy.
One of the chief aims of this work is to compare the GW and Galactic NS mass distributions. We also investigate the maximum mass inferred for the extragalactic NS population, the NS pairing probability for assembling BNSs, and the effect of including or excluding an outlier observation, GW190814 (Abbott et al. 2020b); the nature of its 2.6 M secondary component is challenging to interpret, as it falls into the mass gap between known high-mass NSs (Cromartie et al. 2020) and lowmass black holes (BHs; Thompson et al. 2019).
On the whole, we find that a minimalistic model consisting of a uniform mass distribution with random pairing into BNSs is adequate for explaining the population of NSs observed so far with GWs. Moreover, these first observations already reveal interesting differences between this extragalactic NS population and that which we observe in our own Galaxy.

NEUTRON STARS IN GRAVITATIONAL-WAVE BINARIES
The set of observations informing our population inference consists of the BNS mergers GW170817 and GW190425, and the NSBH mergers GW200105 and GW200115. We separately consider treating GW190814 as another NSBH merger, given the uncertainty in the nature of its secondary component. GWs from compact binary inspirals such as these provide accurate measurements of their source's primary and secondary masses, m 1,2 , as well as its luminosity distance D L . These mea-surements, in the form of a three-dimensional marginal likelihood P (d|m 1 , m 2 , D L ), are the observational input for our population inference. The marginal likelihood can be obtained from the full GW likelihood P (d|m 1 , m 2 , D L , θ) by marginalizing over nuisance parameters θ, such as the source's inclination and spins. Although our NS population model is independent of D L , we must retain the luminosity distance dependence in the likelihood in order to model GW selection effects.
For each of the observations, we calculate the marginal likelihood from the GW posterior samples released by the LIGO-Virgo Collaboration. 1 That is, we weight every posterior sample (m 1 , m 2 , D L ) by a factor of 1/π ∅ (m 1 , m 2 , D L ) to convert from the prior π ∅ (m 1 , m 2 , D L ) adopted for GW parameter estimation, which is uniform in redshifted masses m z 1,2 = (1+z)m 1,2 and quadratic in the luminosity distance, to the uniformin-m 1,2 -D L prior required for the likelihood. Gaussian kernel density estimates of the one-dimensional marginal component mass likelihoods constructed from these weighted samples are plotted in the upper panel of Fig. 1 to illustrate the mass measurements that drive the population inference.

MODELING THE NEUTRON STAR MASS DISTRIBUTION
Motivated by a few basic considerations, we investigate several simple models for the NS mass distri-bution. Firstly, causality arguments place an upper bound of ∼ 3 M on the maximum NS mass (Rhoades & Ruffini 1974;Kalogera & Baym 1996). Plausible supernova formation channels do not produce NSs below 1 M Woosley et al. 2020) and the standard LIGO/Virgo searches target masses greater than 1 M (Abbott et al. 2016;Usman et al. 2016;Mukherjee et al. 2021). Our population models are consequently formulated for m ∈ [1, 3) M .
Second, although there is evidence that recycled and nonrecycled pulsars in the Galactic population originate from different mass distributions (Özel et al. 2012;Farrow et al. 2019), we do not expect the current limited set of GW mass measurements to have the resolving power to discriminate between different subpopulations. We therefore assume that all NSs in BNSs and NSBHs originate from a common mass distribution π(m|λ).
Third, again because of the relatively small number of NS observations at hand, we prefer simple shapesuniform or Gaussian-for the common mass distribution. We also consider a sum of two Gaussians to facilitate a direct comparison with the bimodal Galactic population (Antoniadis et al. 2016;Alsing et al. 2018;Shao et al. 2020).
The astrophysical processes of stellar evolution and binary interaction that eventually result in a merging system containing one or two NSs have significant complexity and uncertainty (Dominik et al. 2012;Fragos et al. 2019;Broekgaarden et al. 2021;Mandel et al. 2021;Santoliquido et al. 2021). In the absence of a definite, astrophysically motivated prescription for the pairing function that connects the masses of the two components in a GW binary, we suppose that the pairing probability is a power law in the mass ratio q = m 2 /m 1 . This phenomenological model is a good descriptor of pairing in the binary BH population (Kovetz et al. 2017;Abbott et al. 2019;Fishbach & Holz 2020;Abbott et al. 2021b). The binary mass distribution characteristic of BNSs and NSBHs is thus taken to be with separate pairing parameters β and β bh for the BNSs and NSBHs. We assume definite a priori classifications for the GW sources. 2 Because there are so few observations of BHs in NSBH systems, we do not attempt to infer the BH mass distribution π bh (m); rather, we fix it to be uniform and nonoverlapping with the NS mass distribution: π bh (m) = U (m|3 M , 30 M ). In practice, the criterion in Eq.
(2) that distinguishes BNSs from NSBHs is m 1 < 3 M . The BH mass distribu-2 If we were instead to classify the sources probabilistically based on their measured parameters, Eq.
(2) would involve a mixing fraction between BNSs and NSBHs that depends on their relative detection rate.
tion's upper bound is chosen to encompass the measurement for GW190814's primary component. We also choose to fix the NSBH pairing parameter to β bh = 0, which corresponds to a minimal assumption of random pairing. These choices do not strongly impact the inferred NS mass distribution, which is driven by the NS mass model and measurements. Indeed, we have tested a different BH mass model-a power law ∝ m −2 -and a different NSBH pairing parameter-β bh = 2, a value that lies within the 90% credible region preferred by the entire compact binary population in the second LIGO-Virgo GW transient catalog (Abbott et al. 2021d,b)and found a negligible impact on the inferred NS mass distribution.
We also take Eq.
(2) to be independent of the GW binary's luminosity distance. This amounts to an assumption that the NS and BH populations do not evolve with distance, which is a reasonable approximation for nearby signals (Fishbach et al. 2018;Mapelli et al. 2019).
Hence, the binary mass model (2) is specified by a choice of basic NS mass distribution π(m|λ) and a choice of BNS pairing parameter β. In light of limited existing GW observations, we make a minimalistic assumption for the latter: we assume random pairing, i.e. β = 0. (We will revisit this assumption in Section 5.) Our three fiducial models, which we call flat, peak and bimodal in Table 1, therefore differ only by the basic NS mass distribution shape they adopt: respectively uniform, Gaussian, and double Gaussian.

GRAVITATIONAL-WAVE POPULATION INFERENCE
Bayesian population inference supplies a prescription for evaluating the likelihood of our set d of GW observations, given a population model π(m 1 , m 2 |λ) that depends on unknown population parameters λ (Mandel 2010;Mandel et al. 2019). Formally, where the three-dimensional marginal likelihood P (d|m 1 , m 2 , D L ) is the observational input from the GWs. The detection fraction which acts as a normalization in Eq.
(3), captures GW selection effects. The selection function Table 1. Inferred Population Parameters for the NS Mass Distribution According to Each of the Models Considered. Medians and symmetric 90% credible intervals of the marginal one-dimensional posterior distributions are reported. The minimum mass is fixed to 1 M for the bimodal and bimodal q models. The AIC statistic, used for ranking the models, is also reported relative to the flat model.
gives the probability of detection for a GW signal as a function of its source properties. We adopt a simple analytic model of a hard signal-to-noise (SNR) ratio threshold following , which produces a selection function that scales like M z 5/2 , where M z = (1+z)(m 1 m 2 ) 3/5 /(m 1 +m 2 ) 1/5 is the redshifted chirp mass. In Eqs. (3)-(4), the distance prior P (D L ) ∝ D L 2 encodes the quadratic increase in a GW survey's ratio of sensitive volume to distance in the local universe.
The posterior probability of a set of population parameters λ within a given population model follows from Eq. (3) via Bayes' theorem: P (λ|d) ∝ P (d|λ)P (λ), where P (λ) is their prior distribution. The model's prediction for the intrinsic NS mass distribution is obtained by marginalizing its binary mass distribution over the posterior uncertainty in λ and averaging over the individual m 1 and m 2 distributions: The corresponding prediction for the observed NS mass distribution P obs λ (m|d)-i.e., the one filtered through selection effects-can be obtained by inserting factors of ζ(λ) within the integrals.

RESULTS AND DISCUSSION
For each of the models described in Sec. 3, we follow the previous section's formalism and evaluate P (λ|d) with a Markov Chain Monte Carlo algorithm implemented with emcee (Foreman-Mackey et al. 2013). We use the marginal likelihoods from GW170817, GW190425, GW200105, and GW200115 as the observational input. We sample from uniform priors over the population parameters, subject to the constraints listed after Eq. (1) (5) is computed from the posterior distribution P (λ|d).
Inferred neutron star mass distribution -We first examine the inferred mass distribution for the three fiducial models described above to discern the key characteristics of the NS population in GW binaries. The median and symmetric 90% confidence interval of the mass distribution inferred for the flat, bimodal and peak models are plotted in Figs. 1, 2 and 3, respectively. These plots also feature traces of individual realizations of the mass distribution drawn from the population posterior. Table 1 reports the median and symmetric 90% confidence interval of the one-dimensional marginal posterior distribution for each population parameter.
The inferred flat mass distribution extends from approximately 1-2 M , with considerable uncertainty in the minimum and maximum mass cutoffs, as expected from only a handful of GW observations;  estimate that 50 BNS signals are needed to determine both cutoffs to within 0.2 M . Nonetheless, the observations constrain the m min and m max parameters away from the prior, disfavoring in particular the largest maximum mass values allowed a priori. The lower bound on m max is driven by the precise mass measurement for GW200105's 1.9 M secondary, which must be accommodated by the mass distribution; this causes the pileup in traces that corresponds to the precipitous fall in the median at 2 M in Fig. 1. The overlap of the less sharply peaked likelihoods in the 1.1-1.6 M range is responsible for the slightly more gradual initial rise in the median.
The mass distributions inferred with the bimodal and peak models are very similar. Both are unimodal Figure 2. Comparison between the inferred mass distributions for NSs in GW binaries and NSs in the Galaxy. Top: posterior predictive check of the bimodal and Galactic population models. One hundred realizations of the observed NS masses are plotted against 100 samples from the predicted bimodal and Galactic mass distributions P obs λ (m|d). The black line with unit slope indicates perfect agreement between the model and the observations; a model that overpredicts (underpredicts) the number of NS observations of a given mass will systematically produce scatter points below (above) the line. Bottom: same as Fig. 1, but for the bimodal model as compared to the median of the Farr & Chatziioannou (2020) Galactic mass distribution. on average, with a broad peak at ∼ 1.6 M ; although many bimodal traces have two peaks, the peak locations are relatively unconstrained by the observational data, which washes out the bimodality in the median. Similarly, despite many individual realizations of the mass distribution having a sharp high-mass cutoff, the GW observations primarily bound m max 2 M for these models, allowing the median to taper off smoothly at high masses. At the low-mass end, we observe a similar tapering in the peak mass distribution, while the bimodal mass distribution is dominated by the imposed sharp cutoff at m min = 1 M .
Besides examining the NS mass distribution predicted by each of the three models, we can attempt to determine which model is preferred by the observations. We adopt the Akaike information criterion (AIC; Akaike 1981), AIC = 2n − 2 log P (d|λ), as the metric for this comparison. The AIC is simply the (logarithm of the) maximal likelihood P (d|λ) for a given population model, whereλ are the maximum-likelihood parameters, plus a correction term that penalizes the model according to its number of free parameters n. The smaller the AIC, the greater the data's preference for the model. The AIC is used for model selection in the astrophysics literature (e.g. Shi et al. 2012;Krishak & Desai 2020), with exp (∆AIC/2) corresponding to the relative likelihood of model B compared to model A for ∆AIC = AIC A − AIC B . Values of ∆AIC > 6 are considered significant for model selection (Liddle 2009).
To rank our population models, for each one we identify the mass distribution realization with the largest likelihood among our discrete population posterior and calculate ∆AIC relative to the flat model. Based on the results listed in Table 1, we find that the peak and bimodal models are moderately disfavored by the data, in part from their larger number of free parameters (n = 4 and n = 6, respectively). The more parsimonious (n = 2) flat model is preferred.
Thus, a uniform NS mass distribution is sufficient to fit current GW observations. There is no clear evidence of bimodality in this GW population of NSs, nor is there unambiguous evidence of a sharp minimum or maximum mass cutoff; these latter features only appear in the inferred flat mass distribution, which includes them by construction.

Comparison with the Galactic neutron star population -
The Galactic population of NSs, observed via electromagnetic astronomy, has been studied extensively (Özel & Freire 2016). The double NS subpopulation appears to conform to a sharply peaked Gaussian mass distribution with µ = 1.33 M and σ ≈ 0.11 M (Kiziltan et al. 2013), at least when assuming a common distribution for recycled and nonrecycled pulsars (Özel et al. 2012;Farrow et al. 2019). The full NS population is characterized by a bimodal mass distribution (Antoniadis et al. 2016;Alsing et al. 2018;Shao et al. 2020). A recent study  found Gaussian peaks at µ = 1.35 +0.04 −0.03 M and µ = 1.8 +0.6 −0.2 M , with respective widths σ = 0.08 +0.03 −0.03 M and σ = 0.3 +0.3 −0.1 M , a maximum mass cutoff at m max = 2.3 +0.8 −0.3 M , and w = 0.7 +0.1 −0.2 . 3 Comparing these two Galactic models with our inferred peak and bimodal mass distributions, which are directly analogous, allows us to investigate the differences between the Galactic and GW populations of NSs.
In Fig. 2, we compare the double-Gaussian mass distribution from  with the bimodal model. Unlike the Galactic mass distribution, the one inferred from GWs is unimodal, predicting far fewer low-mass and moderately more high-mass NSs in the population. The bimodal mass distribution may also support a gentler maximum mass cutoff. The apparent difference at the low-mass end of the Galactic and GW mass distributions is an artifact of our choice to fix m min = 1 M in the bimodal model-consistent with our prior bounds-whereas Similarly, in Fig. 3, we compare the Gaussian double NS mass distribution from Kiziltan et al. (2013) against the peak model. The peak in the latter mass distribution is much broader than the Galactic one and occurs at a higher mass scale. Hence, the peak model predicts many more heavy double NSs in the population than are observed in our Galaxy. This is consistent with the discovery of GW190425, whose total mass of approximately 3.4 M is already recognized as an outlier from the observed Galactic double NS population (Abbott et al. 2020a).
To reinforce the conclusion that the observed populations of NSs in the Galaxy and in GW binaries are distinct, we perform a posterior predictive check of the models. Essentially, we test whether the Galactic mod-els provide an equally good fit to the observed GW population of NSs as our bimodal and peak models. To do so, we draw one mass sample from each of the GW marginal mass likelihoods to constitute a realization of the observed NS population, and we draw the same number of samples from either the GW or Galactic mass distribution P obs λ (m|d). This constitutes the prediction for the detected NS population, accounting for selection effects. Iterating over 100 instances of this procedure, we pair up the observed and predicted masses in sequence of increasing mass, and plot the resulting scatter points in the upper panels of Figs. 2 and 3. A model that fits the observations perfectly would have its scatter points distributed along the line of unit slope. Despite the considerable statistical uncertainty in the mass likelihoods and inferred mass distribution for the GW population models, we see that their scatter points are generally oriented along this line. On the other hand, the scatter points for both Galactic models skew above the line, indicating that they underpredict the number of highmass NSs observed with GWs. This is especially apparent for the double NS model from Kiziltan et al. (2013) in Fig. 3.
Thus, a mass model distinct from the existing Galactic ones is indeed required to accurately describe the detected GW population of NSs. We account for the GW selection bias toward detecting heavy masses, so this may reflect a difference in the observed properties of Galactic pulsars vs. the intrinsic properties of merging compact binaries as a result of different binary evolution pathways. For example, Galaudage et al. (2021) suggest that the rapid merger of high-mass systems reduces their radio visibility, which could explain the relative lack of such systems seen in our Galaxy. More generally, the properties of the observed Galactic population are affected by radio selection effects (Chattopadhyay et al. 2020), which we do not explore in this work. It may be possible to reconcile the two observed populations by building additional parameters such as the spins or the orbital period into the model (Kruckow 2020).
Maximum mass -The best-constrained population parameter in each model is the maximum mass, m max . We investigate whether the maximum mass in the GW population of NSs agrees with (1) the maximum mass in the Galactic population, and (2) the maximum Tolman-Oppenheimer-Volkoff (TOV) mass supported by the NS equation of state (EOS). A discrepancy with the former could be informative about differences between GW and electromagnetic selection effects; a discrepancy with the latter could indicate that the NS mass spectrum is limited by the astrophysical formation channel rather than the EOS.
In Fig. 4, we plot the m max posterior distribution inferred with the flat and bimodal models. Both posteriors peak around 2 M , but the bimodal model has support for m max as large as 3 M because its allows for tapering instead of a sharp cutoff. The flat model predicts a maximum mass of 2.0 +0.4 −0.3 M , while the bimodal model predicts 2.1 +0.7 −0.3 M . The lower bound on m max is driven by the well-resolved mass of the secondary in GW200105 in both cases. The m max posterior for the peak model is almost identical to the bimodal model.
Our estimate of the Galactic NS population's maximum mass comes from the double-Gaussian model of , which predicts m max = 2.3 +0.8 −0.3 M . The corresponding posterior distribution is displayed in Fig. 4. As can be seen, this Galactic maximum mass is completely consistent with the flat and bimodal estimates of m max within current uncertainties: there is no discrepancy in the maximum masses of Galactic and extragalactic NSs.
For the maximum TOV mass-the maximum mass of a nonrotating NS-predicted by the dense-matter EOS, we take the value of M TOV = 2.2 +0.4 −0.2 M obtained by Landry et al. (2020) from nonparametric inference based on radio pulsar, X-ray and GW observations of NSs, an analysis that did not interpret GW190814's secondary as a NS. We find good agreement between the EOSinformed maximum NS mass posterior, also shown in Fig. 4, and that inferred from the GW population with both the flat and bimodal models: the maximum mass for extragalactic NSs is consistent with M TOV , suggesting that stellar and binary evolution are able to produce GW binaries containing the heaviest NSs. This conclusion persists even when GW190814 is incorporated into the population inference as an NSBH observation, as long we also account for it in the analysis of the EOS; as a high-mass outlier, it dictates M TOV .
GW190814 -The results presented up to this point have relied on four GW observations-GW170817, GW190425, GW200105, and GW200115-that are readily interpretable as BNS or NSBH mergers. In this subsection, we also consider GW190814, whose secondary component has a mass of 2.6 +0.1 −0.1 M (Abbott et al. 2020b). Although its nature is uncertain (e.g. Most et al. 2020;Biswas et al. 2021;Tews et al. 2021), and comparisons with various maximum NS mass estimates suggest its secondary may be a low-mass BH (Abbott et al. 2020b;, we now entertain the scenario in which it is a NS and revisit our population inference under this assumption. Because GW190814 is such an outlier relative to the other GW observations, Figure 4. Inferred maximum mass for the population of NSs in GW binaries. Posterior distributions P (mmax|d) for the maximum mass parameter in the flat and bimodal models are compared against the maximum mass inferred from the Galactic NS population in  and the maximum TOV mass from the EOS inference in Landry et al. (2020). The dotted curves show the maximum mass posterior when GW190814's secondary is treated as a NS; we simulate the impact of this assumption on MTOV by discarding samples with MTOV < 2.5 M , the 90% credible lower bound on GW190814's secondary mass, from the Landry et al. (2020) posterior.
its inclusion in the population has a significant impact on the inferred mass distribution.
Its primary effect is to flatten and extend the NS mass distribution beyond 2.5 M , as can be seen in Figs. 1, 2 and 3 when it is included in the flat, bimodal and peak analyses, respectively. Indeed, the m max posterior is particularly sensitive to the inclusion of GW190814. For the flat model, the inferred maximum mass in the GW population of NSs shifts from 2.0 +0.4 −0.3 M to 2.7 +0.2 −0.2 M , as illustrated in Fig. 4. For the bimodal model, a comparable shift from 2.1 +0.7 −0.3 M to 2.7 +0.3 −0.1 M occurs (the m max values for the peak model are very similar). If the population of NSs in GW binaries is taken to include GW190814's secondary component, there is likely a difference between the maximum masses observed in the Galactic and GW populations, although the uncertainties in m max are large enough that the posteriors still overlap. This difference could be indicative of a radio selection effect like the one discussed in Galaudage et al. (2021). Of course, the simplest reconciliation of this discrepancy between the Galactic and GW populations is that GW190814 does not contain a NS, as is suggested by a detailed comparison of its secondary mass with an EOS-informed estimate of M TOV based on other astrophysical observations .
Pairing function -Thus far, our population models have assumed random pairing of NSs into BNSs. We now revise that assumption and investigate mass-ratiodependent pairing as an alternative, using a power-law pairing function as in Fishbach & Holz (2020). We choose β = 2 as our fixed BNS pairing parameter to match our aforementioned alternative NSBH pairing parameter choice β bh = 2. We introduce flat q, peak q and bimodal q models that differ from the fiducial ones only by this choice of β in Eq. (2). For each of these new models, we repeat the analysis described at the beginning of this section and examine how the inferred mass distributions and population parameter constraints change. We then compute AICs to determine whether the data prefer random or q-dependent pairing.
The medians of the inferred mass distributions for the flat q, bimodal q and peak q models are plotted in Figs. 1, 2 and 3, respectively, alongside their randompairing counterparts. In all cases there is little difference between the mass distributions inferred with and without q-dependent pairing. This is also apparent from the population parameter constraints in Table 1.
We compute AICs for the new models, listing their ∆AIC values relative to the flat model in Table 1. The AIC is marginally larger for the flat q, peak q and bimodal q models than their respective random-pairing counterparts. Overall, we find no systematic evidence that a q 2 pairing function is preferred over random pairing for BNSs. This conclusion also holds when we test a more extreme pairing power law (β = 7); GW observations cannot yet constrain the BNS pairing function.

CONCLUSIONS
With this first look at NSs in GW binaries as a population, we find a broad spread of masses over the range compatible with NSs, consistent with a uniform mass distribution. There is no clear evidence in the GW population of the double-peaked structure inferred from Galactic NS mass measurements, nor of the narrow low-mass peak associated specifically with Galactic double NSs. Indeed, we find that high-mass NSs are relatively more prevalent in GW binaries than in the those surveyed with electromagnetic astronomy. However, the maximum NS mass we infer from the GW population agrees with the sharp high-mass cutoff in the Galactic NS mass distribution, provided we do not interpret the outlier event GW190814 as an NSBH. Regardless of GW190814's classification, this maximum mass is consistent with the maximum nonrotating NS mass supported by the dense-matter EOS. The current set of extragalactic NS observations does not allow us to distinguish whether NSs pair randomly or preferentially in equal-mass BNSs. Given the small number of observations currently at hand, near-future NS mass measurements with GWs will have a significant impact on this early picture of the population.
The authors thank Katerina Chatziioannou, Reed Essick, Maya Fishbach, Shanika Galaudage, Richard O'Shaughnessy, Simon Stevenson and the LIGO-Virgo Rates & Populations working group for useful feedback about the manuscript. P.L. is supported by National Science Foundation award PHY-1836734 and by a gift from the Dan Black Family Foundation to the Nicholas & Lee Begovich Center for Gravitational-Wave Physics & Astronomy. J.R. thanks NSF PHY-1836734 and PHY-1806962 and the Carnegie Observatories for support. This research makes use of data, software and/or web tools obtained from the Gravitational Wave Open Science Center (https://www. gw-openscience.org), a service of LIGO Laboratory, the LIGO Scientific Collaboration and the Virgo Collaboration. LIGO is funded by the U.S. National Science Foundation. Virgo is funded by the French Centre National de Recherche Scientifique (CNRS), the Italian Istituto Nazionale della Fisica Nucleare (INFN) and the Dutch Nikhef, with contributions by Polish and Hungarian institutes. The authors are grateful for computational resources provided by the LIGO Laboratory and supported by National Science Foundation grants PHY-0757058 and PHY-0823459. This material is based upon work supported by NSF's LIGO Laboratory, which is a major facility fully funded by the National Science Foundation.