Neutron Star Mass and Radius Measurements

: Constraints on neutron star masses and radii now come from a variety of sources: theoretical and experimental nuclear physics, astrophysical observations including pulsar timing, thermal and bursting X-ray sources, and gravitational waves, and the assumptions inherent to general relativity and causality of the equation of state. These measurements and assumptions also result in restrictions on the dense matter equation of state. The two most important structural parameters of neutron stars are their typical radii, which impacts intermediate densities in the range of one to two times the nuclear saturation density, and the maximum mass, which impacts the densities beyond about three times the saturation density. Especially intriguing has been the multi-messenger event GW170817, the ﬁrst observed binary neutron star merger, which provided direct estimates of both stellar masses and radii as well as an upper bound to the maximum mass.


Introduction
The study of neutron stars represents our best chance to study matter under conditions of high density, extreme isospin asymmetry, and relatively cold temperatures which cannot be examined through heavy ion collisions. As explained below, neutron star matter is in strong-and weak-interaction equilibrium, which for densities larger than the nuclear saturation density, n s 0.16 fm −3 , which is the normal density found inside atomic nuclei, results in very neutron-rich compositions in which the neutron/proton ratio n n /n p is 10 to 20. In laboratory nuclei, and in heavy-ion experiments, on the other hand, one generally has more symmetric matter, n p 0.9n n − n n . Limited information (see Ref. [1] for a review and references) is available through nuclear structure studies of neutron-rich nuclei such as mass, neutron skin, and giant monopole and dipole resonances that can probe cold matter up to n s , but under relatively symmetric conditions; see Section 4. Recent advances in theoretical neutron matter studies [2] complement experiments and probe extremely neutron-rich matter, but are expansions and therefore limited to densities below about 2n s . Measurements of the radii of neutron stars are excellent probes of neutron star matter from n s − 3n s , as the pressure of such matter in this range is highly correlated with the radius [3], while the maximum mass neutron star mass probes densities in the range from 3n s − 5n s . Lower and upper bounds on neutron star radii are found through the assumptions of general relativity and causality [4], as discussed in Sections 2 and 5. Radii can be estimated through studies of thermal X-ray emission from cooling neutron stars, either from young stars or from old, but transiently accreting, sources in quiescence as discussed in Section 6, from pulsar moments of inertia, and from tidal deformabilities and masses inferred from gravitational waves from neutron star-neutron star or black hole-neutron star mergers, as summarized in Sections 7 and 8. Section 2 reviews the lower limit to the neutron star maximum mass that comes from pulsar timing measurements of the most massive pulsars, and an upper limit may be determined from multi-messenger merger events in which both gravitational waves and electromagnetic signatures are seen, another topic of Section 8. Observational data are expected to grow rapidly in the near future due to both ground-based (SKA [5], FAST [6], and LVC) and space-based (HST, NICER [7], ATHENA [8], and the proposed eXTP [9]) observatories.

General Mass and Radius Limits from First Principles
It is useful to examine the most extreme neutron star configurations. For example, Koranda, Stergioulas and Friedman [10] suggested that the most compact stars result from an equation of state (EOS) that is as soft as posssible at low (subsaturation) densities but as stiff as possible at the highest densities. Here, soft and hard refer to the relative pressures p at the corresponding energy densities ε. The extreme case can be modeled with the maximum compactness EOS where the pressure vanishes at low density but is causal (i.e., having a sound speed equal to light speed) at high densities: where ε 0 represents a transition density which is the single parameter of this EOS. When employed in the relativistic stellar structure equations dp dr = − G(m + 4π pr 3 /c 2 )(ε + p) r(rc 2 − 2Gm) , where m is the mass interior to the radius r, one obtains solutions that depend only on the central energy density ε c and the parameter ε 0 . In particular, the solution with the maximum possible mass has central values ε c /ε 0 = 3.034 and p c /ε 0 = 2.034, as well as GM max /(R max c 2 ) = 0.354, where M max and R max are the total mass and radius [11]. Moreover, the maximum mass configuration has M max = 0.08513c 8 / G 3 ε 0 , which shows that ε 0 = 16.74(M max /M ) −2 ε s where M max is the neutron star maximum mass and ε s 150 MeV fm −3 is the energy density at n s . At present, the most massive, accurately-measured pulsars are PSR J0348 + 0432 [12], with M = 2.01 ± 0.04, and PSR J0740 + 6620 [13], with M = 2.17 +0.11 −0.10 M . While the former is technically a more accurate measurement, it is subject to uncertainties in modeling thermal evolutions of white dwarfs, while the latter is a Shapiro delay measurement not requiring theoretical modeling. With time, the uncertainties of Shapiro delay measurements should decrease. Conservatively, we adopt 2.0M for the minimum value of M max , giving ε 0 < ∼ 4.19ε s and ε c < ∼ 12.7ε s . These represent the highest possible values for the energy density in a neutron star; a larger measured M max will lower them.
For masses smaller than the maximum mass, stars become less compact, i.e., they have smaller values of β = GM/(Rc 2 ). From integration of Equation (2), one can determine the dimensionless bounds on R/R max as a function of M/M max using Equation (1). Rescaling with assumed values for M max , the region in the M − R diagram excluded by causality is displayed in Figure 1 where v s is the sound speed. It is commonly believed that the realistic high-density limit of the EOS, in which quarks are asymptotically free, is s = 1/3. In this case, ε c = 4.826ε 0 , p c = 1.275ε 0 , and GM max /(R max c 2 ) = 0.271 in the maximum mass case. Thus, M max = 0.05169c 8 / G 3 ε 0 and ε 0 = 6.172(M max /M ) −2 ε s . Thus, ε 0 = 1.53ε s and ε c = 7.37ε s ; the latter is clearly smaller (by about a factor √ 3) than for s = 1, as expected.  The most rapidly-spinning pulsar, PSR J1748-2446ad, can be used to estimate the maximum radii of neutron stars. It can be shown [16] in general relativity, under the assumption of uniform rotation that the maximum (Keplerian) spin rate, where mass-shedding from the equator occurs, is, for a star of non-rotating mass M and radius R, The spin frequency of PSR J1748-2446ad, 716 Hz [17], implies that R < 13.2(M/M ) 1/3 km for this pulsar, and automatically applies to all stars with masses greater than this pulsar, in particular for the maximum mass configuration. However, although the mass of this pulsar is unknown, if it is much less than M max , this limit should apply to most neutron stars.
Measured neutron star masses via pulsar timing are displayed in Figure 2 [14,18]. The figure shows that the error-weighted average masses (dashed lines) are close to 1.4M for all categories of measurements, but the dispersion in masses is especially small in the case of double neutron star systems. This figure does not include the recently discovered pulsar J0740 + 6620 with M = 2.17 +0.11 −0.10 M [13]. Some other massive estimates are noted, such as the black-widow pulsars B19757 + 20 and PSR J1311-3430, but these have quite large uncertainties due to the unknown sizes and shapes of their companions. These pulsars tend to be massive as they appear to be 'consuming' their companions. Several additional systems are now known, and it is hoped further observations and models will refine these mass estimates. 1 This formula is accurate to a few percent, despite the fact that the equatorial radius increases by about 50% for a maximally-rotating object.

Nuclear Physics Constraints on Neutron Star Radii
Lattimer and Prakash [3] discovered that neutron star radii and neutron star matter pressure around n s − 2n s are highly correlated. Approximately, where p s (p 2s ) is the neutron star matter pressure at n s (2n s ). 2 One can define the symmetry energy as the energy difference between pure neutron matter and symmetric nuclear matter at a given density. Denoting E x (n) as the baryonic energy per baryon as a function of density and proton fraction x, the symmetry energy is S(n) = E 0 (n) − E 1/2 (n). It is usually assumed that the energy of matter at arbitrary proton fraction 0 < x < 1/2 may be found with a quadratic interpolation E x (n) = E 1/2 (n) + S(n)(1 − 2x) 2 . This is equivalent to expanding E x in a power series expansion of 2 These formula were updated by Ref. [1] to comply with the M max > ∼ 2M constraint.
(1 − 2x) and neglecting terms higher than quadratic order, for which there is theoretical justification [19]. With this assumption, the symmetry energy may be expanded about n s : where S V = S(n s ) ∼ 32 MeV, L ∼ 50 MeV and K sym ∼ −100 MeV are symmetry energy parameters. Neutron star matter has a composition determined by energy minimization ∂(E + E e )/∂x = 0, where E e = (3/4)hcx(3π 2 nx) 1/3 is the energy due to electrons. This minimization is equivalent to beta-equilibrium, µ n − µ p = µ e , where the µs are chemical potentials, which leads to x 0.04 at n s . Thus, one can approximate neutron star matter by pure neutron matter, and it follows that the neutron matter energy and pressure at n s are E 0 (n s ) Here, B 16 MeV is the bulk binding energy of nuclear matter at n s . A valuable experimental and theoretical lower limit to the energy of neutron matter and the symmetry energy parameters exists from unitary gas considerations which sets a larger lower limit to radii than does causality. A unitary gas refers to the idealized state of fermions interacting via a pairwise short-range s-wave interaction having an infinite scattering length and zero range [20]. The unitary gas conjecture posits that the lower limit to the neutron matter energy at any density is that of a unitary gas [15]. Pure neutron matter has an s-wave scattering length of a 0 = −18.9 fm, giving (a 0 k F ) −1 = −0.03 at n s , where k F = (3π 2 n) 1/3 is the Fermi wave number, whereas a unitary gas has (a 0 k F ) −1 = 0. Neutron matter also has an effective range r e f f ∼ 2.7 fm, compared to a unitary gas which has r e f f = 0. However, both distinctions, as well as its higher-order spin terms and tensor interactions, make neutron matter more repulsive, with therefore a greater energy, than a unitary gas A unitary gas has the universal behavior E UG = 3ξ 0 E F /5 12.6(n/n s ) 2/3 MeV, P UG = n 2 dE UG /dn 1.35(n/n s ) 5/3 MeV fm −3 , (8) where ξ 0 0.37 is the Bertsch parameter, well-measured in cold atom experiments, and E F = (h 2 k 2 F /2m n ) is the Fermi energy of non-interacting non-relativistic neutrons. The unitary gas conjecture automatically sets the lower limits S V ≥ 28.6 MeV, L ≥ 25.3 MeV and K sym ≥ −275 MeV [15]. This rules out large portions of the L − S V plane (right panel of Figure 1), and establishes that many commonly-used dense matter EOSs have poor choices for the symmetry parameters. The lower limit to the neutron pressure at n s , Ln s /3 ≥ 1.35 MeV fm −3 , also establishes that R 1.4 > ∼ 9.7 km, 1.6 km larger than the causal limit.

Constraints Based on Neutron Matter Theory and Nuclear Experiments
There has been much progress in refining the energy of neutron matter from chiral many-body approaches (see Ref. [15] for a summary of recent work). Some examples of their predictions for the symmetry parameters S V and L are shown in Figure 1, which are seen to be consistent with the unitary gas conjecture, and predict relatively small values of L, in the range 30 MeV< L < 70 MeV, and therefore radii, 10 km < ∼ R 1.4 < ∼ 13 km. Although uncertainties in these calculations increase as the density is increased, error analyses indicate that they are useful up to about 2n s .
The symmetry energy parameters can also be constrained from nuclear measurements. The most accurate measurements are those of nuclear masses. According to the liquid drop model, the nuclear mass can be described primarily as the sum of volume, surface and Coulomb contributions, which scale as A, A 2/3 and (Z/A) 2 A 5/3 , respectively. The volume and surface energies both decrease as the nuclear asymmetry is increased, so that the symmetry energy of a nucleus is positive and proportional to (1 − 2Z/A) 2 : where the subscript i refers to the ith nucleus. Additionally, there are shell and pairing contributions to the total energy, but the symmetry energies can be accurately found, modulo the Coulomb term, by To optimize the symmetry parameters S v and S s , one minimizes the quantity where N is the number of nuclei and σ 1 MeV is a standard error. This results in a narrow correlation ellipse in the S s − S V plane with orientation 10 • with respect to the S s axis, and widths σ v 2.3σ and σ s = 13.2σ [14]. To a reasonable approximation, L/S V = (S s /S V ) 0.7 can be used to propagate this to the L − S V plane, as seen in Figure 3. In a similar way, constraints from neutron skin, dipole polarizability, giant dipole resonance and heavy ion collisions can be overlaid to obtain a relatively small region of overall consistency ( Figure 3) suggesting that 30 MeV < ∼ S v < ∼ 32 MeV and 40 MeV < ∼ L < ∼ 60 MeV.

Extrapolations of the EOS to Higher Densities
Having observed that reasonably strict constraints on the energy density and pressure exist up to about n s from nuclear experiments and up to about 2n s from neutron matter theory, it is interesting to couple these with the global constraints of causality and M max > ∼ 2M to limit the EOS at higher densities. A useful parameterization consists of piecewise polytropes. Read et al. [21] demonstrated that three polytropes in addition to a low-density crust are sufficient to adequately describe most existing models of dense matter from Skyrme-like non-relativistic and Walecka-like relativistic field-theoretical models. In principle, continuity of the EOS at the boundaries between these regions requires eight parameters, for example, the boundary densities and pressures n i , p i for i ∈ 0 − 3. We assume the crust boundary is n 0 = n s /2.7 and p 0 = 0.2177 MeV fm −3 to match the SLy4 EOS [22] that we use at lower densities. 3 Furthermore, we assume n i = 1.85(2 i−1 )n s for i ∈ 1 − 3 following Ref. [21] who found these were expedient choices to match a variety of published EOSs. Therefore, only three parameters remain unspecified. The range of p 1 , 8.4 MeV fm −3 < p 1 < 20 MeV fm −3 can be inferred from neutron matter studies, since n 1 = 1.85n s is small enough that these calculations are still valid. We note that a considerably smaller lower limit to p 1 , 3.74 MeV fm −3 , would follow from the unitary gas conjecture. For the studies reported here, the upper limit to p 1 was arbitrarily increased by 50%. 4 The ranges of p 2 and p 3 are restricted from above by causality and from below by the M max > ∼ 2M constraint. The parameter ranges and the details of this scheme are summarized in Ref. [4].
The regions of M − R and p − ε space that result from sampling these parameters within their allowed ranges are displayed in Figure 4. These ranges are sensitive to the assumed lower limit to M max ; increasing M max raises the minimum R as a function of M and the minimum p as a function of ε. Had the lower limit to p 1 been taken from the unitary gas conjecture, the lower limits to R for M 1.4M would have been reduced by about 0.5 km for M max < 2.1M but effectively unchanged for larger values of M max . Nevertheless, the lower limits to R for all M are substantially larger than those implied by the maximum compactness conjecture (Figure 1).

Astrophysical X-Ray Constraints on Mass and Radius
In general, thermal emissions from neutron stars can be used to infer masses and radii. If neutron stars were non-relativistic blackbodies, measurements of total fluxes, peak spectral energies (which are in the X-ray band), and distances D would suffice to predict the source radius. However, thermal measurements need to be corrected for the surface redshift, their uncertain atmospheres and magnetic field structures, and interstellar absorption that may extinguish a large fraction of their emission before it reaches the Earth. Due to the gravitational redshift, the observed flux is reduced by two redshift factors, one for the energy and one for the time, and the observed temperature T e f f is reduced by a single redshift factor. Since the flux is proportional to (R/D) 2 T 4 e f f , where D is the source distance, the inferred radius is increased by one redshift factor, R ∞ = R/ 1 − 2GM/(Rc 2 ). The atmosphere/magnetic field effectively changes the location of the spectral peak from which the color temperature T c is inferred; this is described by the introduction of a color-correction factor With predictions of f c , one can thus infer the angular area of the star (R ∞ /D) 2 , and with a distance estimate, the apparent radius R ∞ . If one or more spectral lines can be identified, the star's redshift z = (1 − 2β) −1/2 − 1 could be found, allowing both M and R to be inferred. However, it has not yet been possible to observe any such spectral features. An additional complication is that a large fraction of the lower-energy X-rays may be absorbed by the interstellar medium before reaching Earth. Either the amount of absorption, usually described in terms of the column density of hydrogen, N H , is treated as a fitting parameter along with M and R, which increases the uncertainty in a radius estimate, or must be taken from pulsar dispersion measures or observed extinctions along the line-of-sight to each source. To summarize, uncertainties in D, f c and N H dominate the uncertainties in X-ray inferences of neutron star radii.
Focus has therefore been on four classes of sources: (i) isolated neutron stars nearby whose distance can be inferred by parallax, (ii) sources in globular clusters for which estimates of D and N H are available, (iii) certain bursting sources, which have the advantage of providing an additional observable, related to the Eddington luminosity that depends on M and R, and (iv) phase-resolved spectroscopy of periodic flux variations from rotating neutron stars.
The classic case of a nearby isolated source with a measured parallax is RX J1856-4754. The distance is D = 115 ± 8 pc [23], which, coupled with a magnetic hydrogen atmosphere model of Ho et al. [24], gives R ∞ = 14.3 ± 1.0 km. For M = 1.4M , this gives R 1.4 = 11.4 ± 1.1 km. For this value of R ∞ , note that dR/dM −3 km M −1 , so even if the mass is different by 0.25M , the radius would change by only 0.75 km. This source has a featureless spectrum, so its atmospheric composition cannot be confirmed. No other sources have yielded reliable radius estimates.
Globular cluster source measurements are much more prolific. Guillot and Rutledge [25] measured five quiescent (non-bursting) sources, Özel and Freire [18] summarized results for eight sources, and recently, Bogdanov et al. [26] studied an additional source. These sources were found to have best-fit values of R ∞ ranging from 8.6 to 23.0 km, but this wide range can be reduced significantly if pulsar dispersion measures [27] or reddening measurements [28] are used to find N H . In this case, it is found that 9.5 km < R ∞ < 13.7 km result [29]. The average, for a 1.4M star, implies a radius of about 9.6 km. However, some sources may have He rather than H atmospheres, which would increase R ∞ by 3-5 km. Since these sources are likely to have similar masses, the large variations in their estimated radii point to unresolved problems with systematic uncertainties involving interstellar absorption, atmospheric composition, and distance.
In the case of photospheric radius expansion bursts, the radiation pressure from the thermonuclear ignition of surface layers temporarily levitates the photosphere, and it is possible to infer the corresponding Eddington flux necessary for radiation pressure to balance gravity. With surface redshift corrections, where κ is the atmospheric opacity, F Edd can be inferred from the maximum observed flux during a burst. The flux falls rapidly from the peak, and within seconds the apparent emitting angular 2 can be inferred from simultaneous flux and temperature measurements on the 'tail' of the burst. Together with F Edd , A provides the additional information needed to separately estimate M and R, assuming that D, κ and f c estimates are valid. Using this method for six photospheric radius expansion bursts, Özel & Freire [18] obtained a mean radius of 11.1 ± 0.4 km. While the stated uncertainty is gratifyingly small, it is probably an artifact of the model being oversimplified, as one can easily demonstrate.
Observations of many bursts from each source yield values for F Edd and A, and their standard deviations δF Edd and δA. One can construct two parameters α and γ: From α, the stellar compactness can be found: β = (1/4) ± (1/4) √ 1 − 8α, which requires α ≤ 1/8 for a real solution. However, in every burst source observed so far, α exceeds 1/8 by more than δα and often by 2δα. This inconsistency may be caused by one or more of the model's simplifying assumptions: spherical source geometry, constant opacity κ and color-correction factor f c during the burst. Uncertainties in a source's composition and distance also contribute to systematic uncertainties. Furthermore, the redshift factor in Equation (13) may be diminished if F Edd is estimated while the photosphere is extended above the surface. Larger radii, and more realistic uncertainties, result if this redshift factor is eliminated [30].
Related to this burst method is the so-called cooling-tail method [31] which focuses on the late-time evolutions. In this method, one fits both the observed spectra and model spectra with blackbody spectra. In addition to the color-correction factor f c , one introduces the dilution factor w defined so that the bolometric luminosity is conserved where B ν ( f c T e f f ) is the Planck function. It is convenient to define T Edd,∞ = (F Edd,∞ D 2 /(σ B R 2 ∞ )) 1/4 as a parameter that does not depend on distance. Note that T Edd,∞ = (c 3 /(κγ)) 1/4 in terms of the previously-used γ. Another fitting factor introduced is = L/L Edd = (T e f f /T Edd ) 1/4 .
In this method, one defines the blackbody normalization K BB = wR 2 ∞ /D 2 . Note that, if K BB was equal to the normalization A, it would follow that w f 4 c = 1, but this would lead to non-conservation of bolometric flux. All changes to K BB occurring after the burst's peak are due to changes in w. The observed evolution of K BB − F BB on the tail of the burst is equivalent to the evolution of w − w f 4 c , allowing the determination of R 2 ∞ /D 2 and F Edd , and therefore T Edd,∞ . This gives a trajectory in the M − R plane, which still depends on the assumed opacity and gravity used in the model atmosphere models. Fortunately, for the composition there are only two realistic cases: a solar H/He atmosphere or a pure helium atmosphere. They give results so different that only one choice produces reasonable radii. For a particular composition, an assumed distance then determines another trajectory in the M − R plane. The resulting radii determinations are, of course, limited by the distance uncertainties as in the photospheric radius expansion technique. Results with the cooling-tail method indicate radii larger than in the photospheric radius expansion technique [32][33][34], but their similarly small uncertainties will be systematically affected by imprecisely known distances.
M − R measurements can also be deduced from periodic flux oscillations stemming from thermonuclear bursts on the surfaces of neutron stars, which introduce temperature inhomogeneities [35]. The amplitude of the bursts, their deviation from sinusoids, and their dependence on X-ray energy, are sensitive to gravitational light-bending and can be used to probe M − R space. Recently, the Neutron Star Interior Composition and ExploreR (NICER) was attached to the International Space Station with a primary mission to measure neutron star masses and radii to 5% accuracy through phase-resolved spectroscopy [7]. Although the redshift will be the best-determined quantity, oscillation amplitudes vary with X-ray energies in a way that probes a function of M and R orthogonal to redshift. The two pulsars under primary initial consideration are PSR J0437-4715, whose mass (1.44 ± 0.07M ) is already known from pulsar timing, and PSR J0030+0451.
A mid-2020s Chinese-European Enhanced X-ray Timing and Polarimetry (eXTP) space telescope [9] is also planned. This mission will contribute to spin measurements, burst spectra and properties of accretion flows onto neutron stars. The effective surface area will be approximately five times larger than NICER with significantly lower background. A variety of pulsars, including those with known masses in a wide range from 1.20M to 1.93M , can be studied.

Neutron Star Moments of Inertia, Binding Energies, and Deformabilities
The discovery of the system PSR J0737-3039 containing two pulsars will allow the measurement of the moment of inertia I of the more massive pulsar [36] through the additional periastron advance induced by spin-orbit coupling. A measurement of I accurate to about 10% is ultimately expected. Since I ∝ MR 2 , and M is known to better than 0.1%, a tight constraint on the radius should follow. Piecewise polytropes, with constraints of causality and M max > 2M , yield [37]: with +(−) signs indicating upper (lower) bounds, as displayed in the left panel of Figure 5.
For the particular case of PSR J0737-3039A, with M 0737 = 1.3381M , the bounds on R 0737 from a measurement of I 0737 become tighter (right panel of Figure 5). Ref. [37] The binding energy is another structural quantity that either can be measured or is important in analyzing astrophysical observations. Indeed, the first multi-messenger event was SN 1987A, not GW170817, in which about 20 neutrinos in a 10 s burst were observed in the Kamiokande and IMB detectors from a gravitational-collapse supernova in the Large Magellanic Cloud on 23 February 1987. 5 The next galactic supernova will likely result in the detection of tens of thousands of neutrinos lasting for many tens of seconds. The most straightforward observable is the total radiated neutrino energy, all of which originates from the gravitational binding energy BE of the newly-formed neutron star, about GM 2 /R 3 × 10 53 erg. Just as for the moment of inertia, relatively EOS-independent bounds for BE can be found from piecewise polytropes [37] BE/M = −0.00921 + 0.6603β + 0.2651β 2 + (−1.60 ± 4.00)β 4 , valid for β ≥ 0.12 as shown in Figure 6.    [38], and the solid line shows the minimum permitted by the unitary gas constraint [38].
It is useful to generate bounds on binding energy formulated entirely in terms of M, even if its uncertainties are larger than that of Equation (18) The gravitational wave signal from binary neutron star mergers differs from that of black holes because of finite-size effects, revealed primarily through tidal deformabilities. These are described by the tidal Love number k 2 [39][40][41], the proportionality constant between an external tidal field and the quadrupole deformation of a star. It is found by solving an additional first-order differential equation [42] simultaneously integrated with Equation (2), and has values ranging between 0.05 and 0.15 for neutron stars. For black holes, k 2 = 0. It is more useful to define a dimensionless deformability by where β is the compactness parameter. Λ is highly correlated with both M and R (right panel of Figure 6), showing, as for moments of inertia and binding energies, a quasi-universal behavior that is insensitive to the EOS. The piecewise polytrope results show that not only does Λ for a configuration decrease with a high power of the mass, but, for a given mass, it also increases with the same high power of the radius. The results of Ref. [42] suggested for M > 1M , or β ≥ 0.1 that k 2 is roughly proportional to 1/β, so, from the definition of Λ, Equation (20), it is to be expected that Λ ∝ β −6 . Although this relation breaks down for both small and large masses, M < 1M or M > 1.8M , for the expected mass ranges expected in compact mergers, Ref. [38] found with the bounds a = 0.0086 ± 0.0011.

Applications to and Constraints from GW170817
The marvelous multi-messenger observation of the binary neutron star merger GW170817 provided unprecedented information on neutron star properties that can tightly constrain viable EOSs. First, the gravitational wave data alone set limits on the total mass of the inspiralling stars and their radii. Second, this information, coupled with the detections of a gamma-ray burst 1.7 s after the merger, an optical/infrared afterglow thought to be due to significant mass ejection of very neutron-rich matter resulting in the synthesis of high-opacity heavy elements (the r-process), and radio jets at later times, allows an upper limit to M max to be estimated [43].

Inferences from Gravitational Waves
First, we consider what is learned from the gravitational wave data. Tidal deformation results in excess dissipation of orbital energy, speeding up the final stages of the inspiral. Both neutron stars are deformed, and the change in the waveform phase due to tides, to lowest order, is [40] δΦ t (t) = − 117 256 where q = M 2 /M 1 ≤ 1 is the binary mass ratio, M = (M 1 M 2 ) 3/5 (M 1 + M 2 ) −1/5 is the so-called chirp mass, and f (t) is the time-dependent gravitational wave frequency, which is twice the orbital frequency.Λ is a particular combination of the two stars' deformabilities: where Λ 1 and Λ 2 refer to the two stars. GW170817 provided the accurate measurement [44] M = 1.188 ± 0.001M . With the additional GW170817 constraint 0.7 < q < 1, this implies that 2.73M < M 1 + M 2 < 3.05M . The behavior ofΛ as a function of M is very similar to that of Λ as a function of M and is shown in Figure 7. This figure shows thatΛ decreases with a large power of M, approximately 5 to 6, and that EOSs with larger radii, as parameterized by Thus, a precision measurement of M coupled with even a poor estimate ofΛ has the potential of constraining the neutron star radius rather well.
One can understand the general behavior ofΛ using Equation (21), in which case Equation (23) becomes R is a reference radius. A notable result from piecewise polytropes is that stellar radii do not vary appreciably, for a given EOS, in the mass range 1.1M − 1.6M (the average variation is about 1% [38]). With the approximation which is remarkably insensitive to q. For example, ∂Λ/∂q = 0 when q = 1 or q = 0.434, andΛ(q = 0.8)/Λ(q = 1) = 1.013. Although Equation (26) is based on a number of assumptions, it is substantially the same as Equation (24). It predicts, using q = 1 andR = R 1.4 , the constant to be 2 6/5 a = 0.00374 ± 0.000479, which has the same mean value although nearly double the uncertainty. The piecewise polytrope result R 1 R 2 is important because it now follows that where α ∼ 1, which is useful in the analysis of gravitational wave data in which very broad priors of Λ 1 and Λ 2 are normally assumed. Zhao and Lattimer [38] used piecewise polytropes, without the assumption that R 1 = R 2 , to find, in the case that M = 1.188M and q = 0.7 that α = 1 +0.65 (of course, α ≡ 1 for q = 1). In comparison, if no correlation is assumed, α is virtually unconstrained.
The use of Equation (27) effectively reduces the dimensionality of a Markov Chain Monte Carlo gravitational-wave analysis by 1, with a substantial increase in computational speed. The gravitational waves from the recently observed merger of two neutron stars, GW170817, were initially analyzed by the LVC collaboration [44]. The signal was fitted to the Taylor F2 post-Newtonian aligned-spin model [45][46][47][48][49][50] which has 13 parameters. Seven of those parameters are extrinsic, including the sky location, the source's distance, polarization angle and inclination, and the coalescence phase and time. The remaining six parameters are intrinsic, including the masses M 1 and M 2 , dimensionless tidal deformabilities Λ 1 and Λ 2 , and the component's aligned spins χ 1 and χ 2 .
The most accurate result is the measurement of M. If spin and deformation effects are small, to lowest order in a post-Newtonian expansion, the orbital behavior only depends on the frequency dependence of the inspiral waveform such that Here,ḟ is the time derivative of the frequency f (t). The GW170817 analysis by the LVC consortium [44] found a relatively precise measurement M = 1.188 ± 0.001M . 6 However, only a vague restriction of the mass ratio, q > ∼ 0.65, is obtained since this parameter only weakly affects waveform models. Incidentally, the additional measurement of the wave amplitude h(t) allows an estimate of the source distance LVC also determined, to 90% confidence, thatΛ < 800, but established no lower bound. 7 The initial LVC analysis employed conservative assumptions and did not utilize information obtained from electromagnetic observations (which constrain the source location and distance). Furthermore, correlations between Λ 1 and Λ 2 were not considered. This is tantamount to the assumption that the two stars did not necessarily have the same EOS. However, only under somewhat contrived circumstances is it likely that the EOSs describing the two stars would differ appreciably. For example, one star could be a black hole, which, however, the measured chirp mass makes exceedingly unlikely as there is no known formation mechanism to produce a stellar black hole of less mass than M max , or one star could be a hybrid star and the other a conventional star [51]. This would generally require the transition mass to be in between the masses of the two stars; otherwise, both would be hybrid or both would be conventional, and both would have correlated Λs according to Λ 1 q 6 Λ 2 . For a close q, this would be especially unlikely. However, even if only one star is a hybrid star, the EOS of the hybrid star below the hybrid transition density would still be identical to that of the hadronic star, and its EOS at higher densities has to be constrained to satisfy the 2M maximum mass constraint. The masses and deformabilities of the two stars would still retain significant correlations. Another possibility is that the stars occupied separate branches in the M − R diagram because one branch (the hadronic branch) is metastable above the phase transition density. Then, the hybrid branch does not have to satisfy the 2M constraint. However, this would require that, for some reason, the lower mass star would be the one to have transitioned to the quark branch while the more massive star did not, opposite to conventional wisdom. More discussion of twin stars, and earlier references, can be found in Refs. [52,53]. 6 This value takes into account cosmological redshift corrections. 7 This result actually corresponds to the upper 80% confidence bound, as it was defined to be the level above which 10% of the probability remains.
By utilizing the additional information described above, De et al. [54] were able to improve the GW170817 analysis, resulting in tighter constraints on M andΛ (including the discovery of a lower bound toΛ). Three parameters were eliminated by using the electromagnetic source's distance and two-dimensional sky position. Additionally, the prior ranges chosen for the spin parameters χ 1 and χ 2 were taken to be smaller than in the LVC analysis, as motivated by observed double neutron star systems. Finally, Ref. [54] correlated the tidal deformabilites using Equation (27).
As seen in Figure 8, the 80% uncertainty region forΛ, 105 <Λ < 485 is considerably improved compared toΛ < 800 found by Ref. [44]. 8 Later, LVC reanalyzed GW170817 with the EOS correlations as suggested by Ref. [54], together with improved waveform models, and concluded the 90% confidence interval was reduced to 190 <Λ < 630. This compares to 85 <Λ < 640 for Ref. [54]. Figure 8 shows the important result that about 36.2% of theΛ posteriors violate the unitary gas conjecture. It is noteworthy that the EOS correlations introduced by Ref. [55] were determined from EOSs that implicitly satisfy the unitary gas conjecture. If the analysis of Ref. [54] were to also include the unitary gas conjecture, the new 90% confidence interval would change to 240 <Λ < 720, putting the two results into clearer agreement. Figure 8. Neutron star parameters inferred for GW170817 [54], assuming that Λ 1 = Λ 2 q 6 . Contours refer to 68%, 80%, 90% and 95% confidence levels. Left panel: The binary chirp mass M and mass ratio q. Right panel: M and the binary deformabilityΛ. The blue dashed line shows the predicted minimum value from the unitary gas constraint [37]. Attached side panels show the integrated probability distributions over M orΛ; dashed lines show where integrated probabilities reach 2.5%, 5%, 10%, 16%, 50%, 84%, 90%, 95% and 97.5%, respectively.
The uncertainties inΛ and q are correlated. Asymmetric mass ratios are predicted to advance the waveform phase δΦ, while a finite deformability retards the phase (Equation (22) and Figure 8). In addition, positive spins are predicted to advance the waveform phase, so fits to the waveform suffer from a triple degeneracy. Fortunately, it is apparent from Equation (23) thatΛ itself is relatively insensitive to q, so thatΛ is a good proxy for the radii. Assuming q ∼ 1 and inverting Equation (24), one finds For GW170817, one finds R 1.4 (12.5 ± 0.2) Λ /500 1/6 km. For the initial LVC result,Λ < 800, one finds R 1.4 < 13.9 km. Because the EOS uncertainties in Equation (30) are represented as absolute 8 Ref.
[44] used the notation "90% upper bound", above which 10% of the posteriors lay, to refer to the equivalent upper 80% bound shown in Figure 8 since 10% of the posteriors lie above it and 10% lie below the lower 80% bound.
bounds, we can use the upper limit without introducing additional uncertainty. The correlated-EOS LVC 90% upper bound from Ref. [55] is R 1.4 < 13.3 km, while the De et al. [54] result is R 1.4 < 13.6 km. The implications for the masses and radii for the two neutron stars in GW170817 are summarized in Figure 9. The left (right) panel of this figure shows the inferred deformabilites (radii) of the stars as a function their masses. The regions with M > (<)1.365M refer to the more (less) massive star.

Inferences from Multi-Messenger Observations
The remnant formed in the immediate aftermath of the GW170817 merger is believed to have been differentially rotating and containing sufficient angular momentum to be near its mass-shedding limit. Such an object can be hydrodynamically unstable, metastable, or absolutely stable against gravitational collapse depending on its mass and the EOS. The maximum mass thought to be supportable by differential rotation is about M DR > ∼ 1.5M max , while uniform rotation can support less, a maximum of M UR = ξ M max , where M max is the maximum mass of a non-rotating neutron star and ξ ∼ 1.2 [56]. Because of binding energies, which depend on both M and the EOS (Equations (18) or (19)), we consider the baryon mass M B as opposed to the gravitational mass M for this discussion. We then define M B,UR /M B,max = ξ b . It has been found that the factors ξ and ξ b depend relatively weakly on the EOS [57]. The aftermath of the merger depends mostly on the relative values of the total inspiralling baryon mass M BT = M B1 + M B2 , compared to M B,DR and M B,UR , modulo the ejected mass.
If M BT > M B,DR , a prompt collapse to a black hole results, with almost no likelihood for the formation of a gamma-ray burst, mass ejection or radio jets. If M BT < M B,UR , on the other hand, the remnant will be indefinitely stable, as the loss of angular momentum from a uniformly-rotating star is thought to occur on much long timescales. Although a gamma-ray burst, radio jets and ejected mass are all likely in this case, the ejected mass may be "poisoned" by neutrinos emitted from the remnant over long timescales which convert neutrons to protons, rendering synthesis of very heavy elements problematical. The afterglow in this case would tend to be due to "blue" radiation flowing from low-atomic weight, low-opacity material, as opposed to "red" radiation flowing from high-opacity matter consisting of very heavy elements, which is what was observed. In addition, a long-lived remnant will emit vast amounts of dipole radiation from spin-down and magnetic fields, which would enhance the lumninosity radiated by the ejecta. This was apparently not seen in the afterglow from GW170817. We thus infer M B,UR < M BT < M B,DR [43], again, modulo the ejected mass. It is convenient to use Equation (18) for the binding energy since it depends only on M. This situation now