The Multi-messenger Matrix: The Future of Neutron Star Merger Constraints on the Nuclear Equation of State

The discovery of gravitational waves (GWs) from the binary neutron star (BNS) merger GW170817 (Abbott et al. 2017a), and commensurately of electromagnetic (EM) counterpart emission across the frequency spectrum (Abbott et al. 2017b), has initiated a new era of multi-messenger astrophysics. GW170817 was followed by a prompt burst of gamma-rays (Abbott et al. 2017c); thermal radiation from the radioactive decay of heavy r-process nuclei (e.g., Drout et al. 2017; Kasen et al. 2017); and nonthermal X-ray/radio emission (e.g., Margutti et al. 2018), likely produced by an off-axis relativistic jet similar to those favored to produce the short-duration gamma-ray bursts (GRBs) observed on-axis at cosmological distances (e.g., Eichler et al. 1989; Berger 2014). Joint messenger observations of BNS mergers offer an unprecedented opportunity to explore the dynamical processes at work in these cataclysmic events as well as fundamental properties of neutron-rich matter near and above nuclear saturation density.

GW detectors, such as the Laser Interferometer Gravitational-Wave Observatory (LIGO)/Virgo, accurately measure several properties of the inspiraling binary, particularly the chirp mass ${{ \mathcal M }}_{{\rm{c}}}$. They can also be used to measure or constrain the tidal deformability of the inspiraling stars prior to their disruption (e.g., Read et al. 2009; Abbott et al. 2018; Annala et al. 2018; De et al. 2018; Fattoyev et al. 2018; Most et al. 2018; Raithel et al. 2018). However, the current generation of detectors are far less sensitive to the post-merger signal and thus of the ultimate fate of the merger remnant, such as whether and when a black hole is formed (Abbott et al. 2017d). Here, EM observations provide a complementary view. For example, the lifetime of the neutron star (NS) remnant prior to black hole formation may be encoded in the luminosities and colors of the kilonova as the result of differences in the neutron abundance (and thus of the synthesis of lanthanide nuclei) in the ejecta (Metzger & Fernández 2014; Lippuner et al. 2017) and in the ejecta kinetic energy due to magnetar spin-down (Metzger & Bower 2014; Fong et al. 2016; Horesh et al. 2016).

The majority of the ejecta from GW170817, as inferred from the modeling of the kilonova, likely originated in the post-merger phase from accretion disk outflows (e.g., Metzger et al. 2008; Fernández & Metzger 2013; Perego et al. 2014; Just et al. 2015; Siegel & Metzger 2017; Fernández et al. 2019) rather than being the result of a dynamical ejection process. This is important because key properties of the accretion disk, particularly its total mass and therefore its wind ejecta, are primarily functions of the total binary mass M_tot instead of the binary mass ratio q (e.g., Radice et al. 2018b; Coughlin et al. 2018b, although see Kiuchi et al. 2019). Although q strongly affects the properties of the tidal dynamical ejecta, this component was likely a sub-dominant contributor in GW170817. This is a fortunate circumstance because, although q is poorly measured by GW observations, M_tot depends only weakly on q for a precisely known value of ${{ \mathcal M }}_{{\rm{c}}}$ (Figure 1).

Zoom In Zoom Out Reset image size

The type of compact remnant created in the aftermath of a BNS coalescence event depends sensitively on its total mass, M_tot, relative to various threshold masses. These threshold masses depend on the neutron star (NS) equation of state (EOS), particularly those which determine key NS properties such as the maximal Tolman–Oppenheimer–Volkoff (TOV) mass M_TOV and the radius for a characteristic mass, e.g., R_1.6 for a 1.6 M_⊙ star. The most massive binaries undergo "prompt" collapse to a black hole on the dynamical timescale (e.g., Shibata & Taniguchi 2006; Bauswein et al. 2013a). Slightly less massive binaries produce hypermassive NS (HMNS) remnants, which are temporarily stabilized from gravitational collapse by their rapid and differential rotation (e.g., Kastaun & Galeazzi 2015; Hanauske et al. 2017). Less massive binaries leave supramassive NS (SMNS) remnants, which remain supported against collapse even once in a state of rigid body rotation after their differential rotation has been removed (e.g., by stresses due to strong internal magnetic fields). Finally, the lowest mass binaries could in principle create remnants that are below the TOV mass and therefore indefinitely stable (e.g., Giacomazzo & Perna 2013), even once their entire angular momentum has been removed (over potentially much longer timescales, e.g., by magnetic dipole spin-down).

The expectation from numerical relativity that total mass is the dominant binary property that determines the outcome of a merger (a "Vogt–Russell" theorem for BNS mergers) is one that a large sample of future mergers with both M_tot measurements and EM counterparts could test. Likewise, with such a quantitative theoretical framework firmly in place, determination of whether an indefinitely stable NS, a SMNS, HMNS, or prompt occurs directly links EM and GW observation, permiting constraints on unknown NS EOS properties. This method was utilized following the first BNS merger GW170817, where the measured ${M}_{\mathrm{tot}}={2.74}_{-0.01}^{+0.04}\,{M}_{\odot }$ (Abbott et al. 2019) along with the inference of a HMNS remnant from the kilonova colors and kinetic energy was used to place an upper limit of M_TOV ≲ 2.17 M_⊙ (Margalit & Metzger 2017; Shibata et al. 2017, and later corroborated by Rezzolla et al. 2018; Ruiz et al. 2018) as well as a lower limit on the NS radius R_1.6 ≳ 10.3–10.7 km (Bauswein et al. 2017; Radice et al. 2018b).⁴

Here, we more systematically develop the multi-messenger approach to constraining the NS EOS. A primary goal of this work is to elucidate the landscape of what can be learned about the EOS by a large sample of BNS mergers with well-characterized EM counterparts. We further use this information to motivate search strategies and provide recommendations for information that LIGO/Virgo should promptly disseminate to the community of EM observers to maximize the scientific impact of its discoveries.

Figure 2 shows the current parameter space of GW and EM observables, delineated by the class of merger remnant. The total remnant mass (vertical axis), depicted by the black curve, is a function of the GW-measured chirp mass ${{ \mathcal M }}_{{\rm{c}}}$ (horizontal axis), and only weakly dependent on the mass ratio q. The EOS sets mass cuts that determine the post-merger remnant for a given ${{ \mathcal M }}_{{\rm{c}}}$. Thus, different remnant types are separated along the vertical axis. Given that different remnants can imprint different observational signatures on the EM counterpart, this vertical axis can also be read as an "EM signature" axis.

Zoom In Zoom Out Reset image size

The EM signatures listed in Figure 2 include the possible presence or absence of a relativistic GRB jet⁵ (e.g., Murguia-Berthier et al. 2017); the kilonova colors and the quantity of ejecta mass M_ej, which depend on the NS remnant lifetime (e.g., Metzger & Fernández 2014); and the ejecta energy (e.g., enhancement due to energy injection from a long-lived magnetar remnant; Metzger & Piro 2014). Figure 3 shows estimates of two measurable properties of the ejecta, particularly the total ejecta mass (measurable from the kilonova) and ejecta kinetic energy (measurable from the kilonova or synchrotron emission from interaction with the interstellar medium), as a function of the chirp mass for a range of q. The kilonova color is qualitatively shown in the bottom panel as well.

Zoom In Zoom Out Reset image size

Uncertainties in our knowledge of the NS EOS translates into an ambiguity in the exact mass cut transitions between remnant types, or an overlap between possible merger outcomes (i.e., multiple possible EM observational signatures) given a measured GW chirp mass. In the idealized approach adopted in this Letter, this overlap is entirely a consequence of our uncertainty of the EOS.⁶ In other words, if the underlying EOS were precisely known, different measured chirp masses would have a one-to-one correspondence to differing EM signatures (different remnant types). Regions in which no overlap between possible remnant types exists should therefore be interpreted as predictions—if LIGO/Virgo measure a chirp mass in this region, our current understanding of the EOS implies only one possible merger outcome, giving a clear prediction for the type of EM counterpart expected to accompany such GW detection. Conversely, chirp-mass parameter space regions in which there is an overlap in the possible merger product should be seen as learning opportunities to constrain the EOS—if an EM counterpart is observed in such an event and can be classified into one of the categories along the vertical axis, this would allow an inference of the merger outcome and break the ambiguity arising from our ignorance of the NS EOS.

2.1. Merger Outcome Mass Cuts

The cuts in chirp mass shown in Figure 2 are calculated from several criteria, mostly based on the remnant mass implied by a given chirp mass. The remnant's baryonic mass depends on the well-measured chirp mass ${{ \mathcal M }}_{{\rm{c}}}$, and weakly on the binary mass ratio q, but is also a function of the EOS through the translation between gravitational and baryonic masses:

$$\begin{eqnarray}&&{M}_{\mathrm{rem}}^{{\rm{b}}}=f({m}_{1})+f({m}_{2})-{M}_{\mathrm{ej}},\end{eqnarray} \tag{ 1 }$$

where M_ej is the ejecta mass lost from the system,

$$\begin{eqnarray}&&{m}_{1}={{ \mathcal M }}_{{\rm{c}}}{q}^{-3/5}{(1+q)}^{1/5},\,\,\,{m}_{2}={m}_{1}q,\end{eqnarray} \tag{ 2 }$$

are the masses of the two merging NSs, and

$$\begin{eqnarray}&&{m}^{{\rm{b}}}\equiv f(m)\approx m+0.058{m}^{2}+0.013{m}^{3}\end{eqnarray} \tag{ 3 }$$

is an approximate conversion between gravitational mass m and baryonic mass m^b in units of solar mass accurate to within 10% in binding energy (${m}^{{\rm{b}}}-m$) for a winde range of piecewise polytropic EOSs (J. Lattimer 2019, private communication; see also Lattimer & Prakash 2016; Zhao & Lattimer 2018). Although, in detail, the function f(m) is not universal and depends on the underlying EOS, we use this approximate expression in our current analysis for simplicity. This is a reasonable approximation given that a ±10% deviation in binding energy induces only <1% variation in the chirp-mass thresholds derived below (see Table 1).

Table 1.  Merger Remnant Scenarios, Based on EOS Constraints Available after GW170817

${{ \mathcal M }}_{{\rm{c}}}$ a (M⊙)	Remnant Typeb	EOS Constraintc	Predicted Fractiond
$0.948\times [1+0.390{M}_{\mathrm{ej}}-\delta (q)$ e]	NS (i-A)	⋯	≪1%
$1.032\times \left[1+0.353{M}_{\mathrm{ej}}-\delta (q)\right]$	NS (ii-A)	MTOV lower bound ≥ 2.01 M⊙	3%
	SMNS (ii-B)	MTOV upper bound ≤ 2.17 M⊙
$1.103\times \left[1+0.325{M}_{\mathrm{ej}}-\delta (q)\right]$	SMNS (iii-B)	⋯	18%
1.188	SMNS (iv-B)	MTOV lower bound ≥ 2.01 M⊙	47%
	HMNS (iv-C)	MTOV upper bound ≤ 2.17 M⊙
$1.433\times {\left[1+{2}^{-1/5}\delta (q)\right]}^{-1}$	HMNS (vi-C)	R1.6 lower bound ≥ 10.3 km	32%
	prompt collapse (vi-D)	R1.6 upper bound
$1.750\times {q}^{3/5}{\left[(1+q)/2\right]}^{-1/5}$	prompt collapse (vii-D)	⋯	≪1%
$1.889\times {q}^{3/5}{\left[(1+q)/2\right]}^{-1/5}$	prompt collapse (viii-D)	⋯	⋯
	BH–BH, NS–BH (viii-E)	MTOV upper bound ≤ 2.17 M⊙f
	BH–BH, NS–BH (ix-E)	⋯	⋯

Notes.

^aThreshold chirp masses for different merger outcomes. Systematic uncertainty of ±10% in the approximate binding energy relation (Equation (3)) induces <1% variation in the first three ${{ \mathcal M }}_{{\rm{c}}}$ cuts (the rest are unaffected). ^bType of compact remnant formed from the merger, labeled following the scheme from Figure 2. ^cEOS constraint obtained by ascertaining the type of merger remnant from EM observations. ^dCalculated assuming a cosmological population of BNS mergers with a mass distribution following that of the known Galactic BNS population. ^eHere $\delta (q)\equiv \left[{q}^{-3/5}{\left(1+q\right)}^{6/5}-{2}^{6/5}\right]/2$ is a small correction factor such that $\delta (q)\lesssim 0.02$ for 0.7 ≤ q ≤ 1, and M_ej is the ejecta mass in units of M_⊙. ^fOr evidence of a formation channel for BHs with M < M_TOV.

Download table as:  ASCIITypeset image

The maximal mass of a non-rotating NS (the TOV mass) is one of the poorly constrained properties of NSs that depends sensitively on the pressure at the highest densities achieved in the NS core. Extant constraints on its value, ${M}_{\mathrm{TOV},\min }\lesssim {M}_{\mathrm{TOV}}\lesssim {M}_{\mathrm{TOV},\max }$ are nevertheless sufficient information for defining the different mass cuts illustrated in Figure 2.

If the mass of the merger remnant is less than the minimum allowed TOV mass, ${M}_{\mathrm{TOV},\min }$, then the remnant is itself guaranteed to be indefinitely stable. This criterion can be translated into a cut on chirp mass, ${M}_{\mathrm{rem}}^{{\rm{b}}}({{ \mathcal M }}_{{\rm{c}}})\lt {M}_{\mathrm{TOV},\min }^{{\rm{b}}}$ (region (i) in Figure 2). For slightly larger chirp/remnant masses we have ${M}_{\mathrm{TOV},\min }\lt {M}_{\mathrm{rem}}\lt {M}_{\mathrm{TOV},\max }$, making it unclear whether M_rem is above or below the stable threshold mass, M_TOV. If EM observations can distinguish a stable remnant (${M}_{\mathrm{rem}}\lt {M}_{\mathrm{TOV}}$) from one forming an unstable SMNS (${M}_{\mathrm{rem}}\gt {M}_{\mathrm{TOV}}$), then a new upper (lower) limit on M_TOV is obtained. Such a distinction could be made, for instance, based on the inferred kinetic energy of the merger ejecta (e.g., Metzger & Bower 2014), or from evidence for BH formation arising from an abrupt drop in the X-ray light curve (e.g., Rowlinson et al. 2010).

For ${M}_{\mathrm{rem}}\gt {M}_{\mathrm{TOV},\max }$ a stable NS cannot form. The transition between an SMNS and HMNS remnant is well approximated by ${M}_{\mathrm{rem}}=\xi {M}_{\mathrm{TOV}}$ with the constant of proportionality ξ ≃ 1.18 largely independent on the EOS (Cook et al. 1994; Lasota et al. 1996). Thus, for ${M}_{\mathrm{rem}}\lt \xi {M}_{\mathrm{TOV},\min }$ we do not expect a HMNS as the final pre-collapse merger remnant type. As a result, the region in which ${M}_{\mathrm{TOV},\max }\lt {M}_{\mathrm{rem}}\lt \xi {M}_{\mathrm{TOV},\min }$ is predicted to have only one possible merger outcome—an SMNS (region (iii) in Figure 2).

The remnants of mergers with slightly larger inferred masses satisfying $\xi {M}_{\mathrm{TOV},\min }\lt {M}_{\mathrm{rem}}\lt \xi {M}_{\mathrm{TOV},\max }$ could form either an SMNS or HMNS remnant (region (iv)). Again, if a BNS merger detected in this mass range with a well-characterized EM counterpart could disentangle these two possibilities, then this would provide a more stringent constraint on M_TOV. This was the strategy employed following the inference that GW170817 likely formed an HMNS remnant to constrain ${M}_{\mathrm{TOV},\max }$ (e.g., Margalit & Metzger 2017). This inference was made based on the quantity, kinetic energy, and composition of the kilonova ejecta (Figure 3).

For a more massive merging NS system, the remnant may undergo a prompt collapse to a BH, which will be accompanied by systematically less massive and more neutron-rich ejecta, resulting in a less luminous and redder kilonova than for the HMNS case. The threshold for such collapse depends on the dense-matter EOS parameterized through M_TOV, but also on a compactness parameter $\propto {M}_{\mathrm{TOV}}/{R}_{1.6}$. Numerical and semi-analytic work have investigated this dependence (Bauswein et al. 2013b; Bauswein & Stergioulas 2017) and found that the remnant will promptly collapse if

$$\begin{eqnarray}&&{M}_{\mathrm{tot}}\gt {M}_{\mathrm{thresh}}\approx \left(2.38-3.606\displaystyle \frac{{{GM}}_{\mathrm{TOV}}}{{c}^{2}{R}_{1.6}}\right){M}_{\mathrm{TOV}},\end{eqnarray} \tag{ 4 }$$

where R_1.6 is the radius of a 1.6 M_⊙ NS. The threshold limit above which the merger undergoes a prompt-collapse is therefore an increasing function of both M_TOV and R_1.6. An HMNS remnant can be ruled out for ${M}_{\mathrm{tot}}({{ \mathcal M }}_{{\rm{c}}})\gt \max {M}_{\mathrm{thresh}}$, and similarly a prompt collapse is ruled out if ${M}_{\mathrm{tot}}({{ \mathcal M }}_{{\rm{c}}})\,\lt \min {M}_{\mathrm{thresh}}$, where the maximum (minimum) is determined by the maximal (minimal) NS radius and TOV mass allowed based on extant constraints, or by the lowest (highest) chirp-mass merger observed to undergo prompt collapse (form an HMNS). The sharp drop-off in the predicted ejecta mass approaching the prompt collapse threshold (Figure 3) should produce a kilonova signature readily distinguishable from the HMNS case. This has been used to constrain R_1.6 following GW170817 (Bauswein et al. 2017), and can further be applied to future detections (Bauswein et al. 2017; see also Table 1).

Finally, under standard astrophysical formation scenarios (e.g., excluding primordial BHs), BHs are not expected to form with masses less than M_TOV, and hence the region where ${m}_{1}\lt {M}_{\mathrm{TOV},\min }$ is ruled out for BH–BH or BH–NS mergers.⁷

Table 1 summarizes the ranges in chirp mass that delineate different predicted remnant types, as well as the EOS constraints that could be obtained from distinguishing the EM counterpart for a GW event in each category. This table is calculated based on extant constraints on the NS EOS following GW170817: ${M}_{\mathrm{TOV},\min }=2.01\,{M}_{\odot }$ (Antoniadis et al. 2013); ${M}_{\mathrm{TOV},\max }=2.17\,{M}_{\odot }$ (Margalit & Metzger 2017); ${R}_{1.6,\min }=10.3$ km (Bauswein et al. 2017); ${R}_{1.6,\max }=13.5$ km (De et al. 2018), but can be updated following future GW/EM events. Except for the two highest chirp-mass cuts related to the BH/NS mass transition, the chirp-mass boundaries are only weakly dependant on binary mass ratio and ejecta mass. This dependence is also given (to first order) in Table 1.

The previous section outlined the extant parameter space and illustrated the predicted outcomes for the next binary NS merger given our current knowledge and uncertainty regarding the NS EOS. We now address the implications of these results applied to a large population of merging binary NSs.

3.1. Merger Outcome Statistics

In order to model the population of extragalactic merging binary NS mergers, we assume the latter population is similar to that of known Galactic double NSs. The inferred NS masses for GW170817 are consistent with being drawn from this population (e.g., Abbott et al. 2019; see also Figure 4), which—albeit for a sample size of N = 1—supports this naive assumption.

Zoom In Zoom Out Reset image size

The Galactic double NS population is narrowly peaked around ${m}_{\mathrm{1,2}}\simeq 1.32\,{M}_{\odot }$ with q ≈ 1 (e.g., Zhao & Lattimer 2018). We draw NS masses from the Gaussian Galactic distribution found by Kiziltan et al. (2013),⁸ weighted by ${{ \mathcal M }}_{{\rm{c}}}^{5/2}$. We include this chirp mass weighting for completeness, as it accounts for the GW detectors' bias of detecting higher-mass systems with larger GW-strain amplitudes (however, in practice we find negligible effect of this weighting on our final results). Our analysis is similar to that of Piro et al. (2017), with the important exception that Piro et al. (2017) examined specific example EOSs, thus fixing M_TOV and R_1.6. Here we instead calculate the statistics based on our current knowledge and uncertainty in the EOS. This gives realistic bounds on the fraction of different merger remnant types, and additionally allows us to estimate the number of mergers that will occur in regions of ambiguity in the remnant fate and thus can be used to further constrain the EOS.

Figure 4 and and final column of Table 1 provide the expected fraction of mergers that will produce a given remnant. We predict that the vast majority of mergers will result in either an HMNS or SMNS remnant (≈18%–68% and ≈0%–79%, respectively). By contrast, the fraction that undergo prompt collapse could range from ≈0% to 32%. Likewise, the fraction that end up as indefinitely stable remnants could range from ≈0% to 3%. A potentially high fraction of SMNS remnants is consistent with findings based on analysis of short GRB X-ray plateaus (Lasky et al. 2014; Gao et al. 2016) or temporally extended prompt X-ray emission (Norris & Bonnell 2006; Perley et al. 2009), assuming that the latter are indicative of the presence of long-lived magnetar SMNS remnants (e.g., Metzger et al. 2008; Rowlinson et al. 2013; Zhang 2013). The rarity of indefinitely stable remnants, which may inject enormous rotational energy ∼10⁵³ erg into the merger environment (Figure 3), is consistent with the tight limits available from radio transient surveys (Metzger et al. 2015).

3.2. Future EOS Constraints

We now consider what improvement could be made in EOS constraints with a sample of future joint messenger detections. Figure 5 shows the expected constraints on M_TOV and R_1.6 as a function of the number of binary NS merger detections with sufficiently well-observed EM counterparts to accurately ascertain the remnant outcome. For this calculation we have assumed a particular EOS as being the "correct" one, for which M_TOV = 2.1 M_⊙ and R_1.6 = 11 km (horizontal dashed gray curves in Figure 5). We then draw a sample of N − 1 binary NS masses from the Galactic binary distribution, as described in the previous subsection. This simulates a set of N successive NS merger detections (the first of which is always taken to be GW170817). Given the chirp mass of the merging system and assuming that the observed EM counterpart is able to break the degeneracy between different possible merger remnants (see Figures 2, 3), we use the methods of Margalit & Metzger (2017) and Bauswein et al. (2017) to update the constraints on M_TOV and R_1.6 following each detection. We then repeat this process for a newly drawn set of N − 1 binary masses, creating a large ensemble of such detection sets. The median upper (lower) bounds on M_TOV and R_1.6 obtained following N detections are then shown by the dark red (blue) curves, with 68% and 95% percentiles around this median depicted by dark and light shadowed regions, respectively (where the median and percentiles are calculated over the ensemble).

Zoom In Zoom Out Reset image size

We find that after N ≳ 10 mergers, one could in principle constrain the value of M_TOV and R_1.6 to within a few percent, close to the level where systematic uncertainties will set in (see the discussion below). The convergence on R_1.6 is sensitive to the true NS radius (e.g., significantly slower for R_1.6 = 12 km) because it is only constrained by the prompt collapse condition (Bauswein et al. 2017), and thus depends on whether the true M_thresh (Equation (4)) is low enough to generate a large population of prompt-collapse events.

In Figure 5 we have accounted only for EOS constraints obtained using the BNS merger multi-messenger approach described above (e.g., the R_1.6 constraint does not account for potential tidal-deformability limits placed by future BNS mergers). Additional NS observations can and should be folded in to even further constrain the EOS. For example, well-measured masses of massive NSs (Demorest et al. 2010; Antoniadis et al. 2013; Thankful Cromartie et al. 2019), experimental laboratory data, NS radius measurements, and the multi-messenger BNS merger constraints described in this work can be combined using the statistical framework described by Miller et al. (2019), providing more stringent constraints than any single method used alone.

Following the groundbreaking discovery, follow-up, and multi-messenger detection of the first binary NS merger GW170817, a flurry of theoretical works have used this watershed event to place constraints on fundamental physics, including the NS EOS. Several of the methods and important results in this field came from analyzing the multi-messenger nature of this event and not, as previously envisioned, solely through GW waveform analysis. These newly applicable methods have illustrated that important scientific results can be obtained from even a single well-observed event such as GW170817. This is in contrast to the post-merger or inspiral GW constraints on the EOS, which are typically thought to require multiple events to build-up significant signal-to-noise and place meaningful constraints on the EOS (e.g., Chatziioannou et al. 2017; Torres-Rivas et al. 2019).

Given the success of the multi-messenger approach for GW170817, now is an opportune time to explore the path ahead for this budding field. The ongoing LIGO/Virgo O3 run may detect between ∼1 and 10 additional BNS mergers, while such discoveries could occur as frequently as once per week following the planned LIGO A+ upgrade. We have here summarized the current methodology for placing constraints on the NS EOS via such multi-messenger observations and laid out the post-GW170817 parameter space, highlighting regions of particular importance for EOS constraints, as well as providing explicit falsifiable predictions of our methods.

Our envisioned program relies on extensive EM follow-up of GW-detected mergers, particularly those with chirp masses indicating that they will provide the most leverage on the EOS if the outcome can be constrained ("learning opportunity"). Given the finite number of observing hours that can be dedicated to EM searches of LIGO/Virgo localization regions, and the large number of near-future GW detections expected, it is inevitable that not all mergers will be followed-up with the most sensitive wide-field optical/near-IR facilities. It is therefore crucial to identify them as early as possible and, before a follow-up campaign is launched or aborted, also identify the possible scientific impact of obtaining multi-messenger observations for each event.

An underappreciated inference from GW170817 is that ejecta from the disk wind, rather than dynamically ejected material, is likely the most important source of the observed kilonova emission (Figure 3). The properties of the kilonova emission is therefore likely to be more sensitive to the total binary mass (a proxy for stability of the remnant) than to the binary mass ratio q. This is important because, while q is poorly constrained by GW data, M_tot is tightly related to the well-measured chirp mass (Figure 1). As shown in Table 1, if the GW-measured chirp mass is indeed the key indicator of the expected EM counterpart, then certain ranges of chirp mass values are more important than others in obtaining new and meaningful constraints on the NS EOS through the multi-messenger analysis. We therefore strongly advocate in favor of releasing chirp mass estimates as early as possible to the scientific community, ideally already with the initial LIGO/Virgo triggers, so that the observational community will be able to allocate its resources in the most efficient way to maximize the scientific gain from multi-messenger observations.

Focusing on an idealized scenario, in which the binary mass is measured to infinite precision and the EM counterparts provide an accurate discriminator between different merger remnants, we find that the radius and maximum mass of NSs could be constrained to within a few percent after ≳10 events (Figure 5). However, a number of systematic uncertainties will enter at the percent level, if not earlier, which limit the precision of the constraints. First, although the chirp mass is generally measured to effectively infinite precision, percent-level uncertainty in the translation of ${{ \mathcal M }}_{{\rm{c}}}$ to M_tot results because of the much larger uncertainty in q (Figure 1). Determining M_rem accurately also requires subtracting off M_ej (which is at best constrained from the kilonova emission to a factor of ≲2; e.g., Wu et al. 2019), and accounting for the NS binding energy (Equations (1), (3)) which implicitly depends on the true EOS. Furthermore, the intrinsic chirp mass is degenerate with the redshift of the event; therefore, as merger events begin to probe cosmological-scale volumes, uncertainties in the cosmological distance scale will ultimately propagate to those in the EOS.

The largest source of uncertainty in our outlined procedure likely come from our incomplete understanding of the merger physics. We thus note the obvious need for further theoretical investigations into the late-time evolution of NS merger remnants. These efforts should be two-fold—focusing on both the post-merger dynamics and remnant formation (this has implications for the precise mass cuts separating different remnant classes) and an enhanced characterization of the EM signatures that may identify each remnant type (including a better understanding of natural sources of variability and possible contaminants).

Theoretical maps like Figure 3 are based on simplified assumptions that do not account for some physical effects that could introduce quantitative uncertainties in the EOS constraints. These include thermal effects on the EOS over the neutrino cooling timescale ∼1 s (e.g., Kaplan et al. 2014), and the precise way in which angular momentum is redistributed throughout the merger remnant. For instance, Fujibayashi et al. (2018) and Radice et al. (2018a) included the effects of an artificial α-viscosity into GR hydro simulations of the post-merger evolution to mimic physical angular momentum transport; they found that the presence of additional ejecta sources was not included in the usual "dynamical + disk wind" dichotomy. In a recent study, Kiuchi et al. (2019) show that mass ratio can affect the prompt-collapse threshold and disk mass estimates, potentially motivating a revised examination of the prompt-collapse criteria.

The "Vogt–Russell" theorem for BNS mergers, which underpins the strategy advanced here, is motivated by numerical relativity simulations, but must itself be tested with a large ensemble of events. This highlights the importance of characterizing the EM counterparts of not just mergers in the mass range that tighten EOS constraints, but those in mass ranges for which the merger outcome is "predicted." For instance, our approach is calibrated on the results of numerical simulations that largely use a standard hadronic equation of state, which ignores the possibility of a first-order phase transition at high densities or temperatures to deconfined quarks (e.g., Drago & Pagliara 2018; Han & Steiner 2019; Paschalidis et al. 2018; Bauswein et al. 2019; Most et al. 2019). Deviations in an ensemble of mergers from the predictions of a "well-behaved" EOS might provide evidence for such a phase transition.

Finally, we note that the multi-messenger methods described here rely solely on a qualitative categorization of which remnant type was produced by the merger. Even stronger EOS constraints can in principle be obtained by quantitatively fitting inferred kilonova properties to numerical relativity simulations (Radice et al. 2018b; Coughlin et al. 2018a, 2018b; Radice & Dai 2019; Lazzati & Perna 2019). Such methods are somewhat more susceptible to uncertainties in the kilonova modeling and numerical relativity simulations, but are a promising way forward.

On the opposite end, GW-waveform analysis alone can provide an extremely clean model-independent probe of the NS radius. However, this method requires very high signal-to-noise events, especially at ∼kHz frequencies. Events as loud as GW170817 should be rare, implying that such constraints may not improve quickly. In contrast, the multi-messenger methods discussed above require only information about the chirp mass, which should be well measured even for lower signal-to-noise events, although much will hinge on the success rate of detecting and characterizing the EM counterparts to the requisite fidelity. Dedicated follow-up programs on large survey instruments, such as the Large Synoptic Survey Telescope, could play a major role in the success of this endeavor (Cowperthwaite et al. 2019).

We thank Jim Lattimer, Andreas Bauswein, and Kenta Kiuchi for helpful comments on the manuscript. B.M. is supported by NASA through the NASA Hubble Fellowship grant #HST-HF2-51412.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. B.D.M. is supported in part by NASA through the Astrophysics Theory Program (grant No. NNX16AB30G). The final stages of this work benefited from conversations at KITP, supported in part by the National Science Foundation under grant No. NSF PHY-1748958.