Measurement of the nuclear symmetry energy parameters from gravitational wave events

TL;DR

This paper uses gravitational wave data from GW170817 to constrain nuclear symmetry energy parameters, especially the slope L_0, revealing lower values than previously reported and providing new insights into nuclear matter properties.

Contribution

It introduces novel constraints on the symmetry energy parameters from gravitational wave observations, focusing on the slope L_0 and deriving bounds on higher-order nuclear terms.

Findings

Gravitational wave data are sensitive to the slope L_0 of the symmetry energy.

The analysis favors lower L_0 values around 20 MeV.

New analytic constraints on higher-order nuclear symmetry energy terms.

Abstract

The nuclear symmetry energy plays a role in determining both the nuclear properties of terrestrial matter as well as the astrophysical properties of neutron stars. The first measurement of the neutron star tidal deformability, from gravitational wave event GW170817, provides a new way of probing the symmetry energy. In this work, we report on new constraints on the symmetry energy from GW170817. We focus in particular on the low-order coefficients: namely, the value of the symmetry energy at the nuclear saturation density, S_0, and the slope of the symmetry energy, L_0. We find that the gravitational wave data are relatively insensitive to S_0, but that they depend strongly on L_0 and point to lower values of L_0 than have previously been reported, with a peak likelihood near L_0 ~ 20 MeV. Finally, we use the inferred posteriors on L_0 to derive new analytic constraints on higher-order nuclear terms.

Equations

$$E_b(n,Y_p)=E_0(n)+E_{\mathrm{sym}}(n)(1-2Y_p{)}^2,$$

$$E_0(n)=B_0+\frac{K_0}{18}{u}^2+\frac{Q_0}{162}{u}^3+\mathcal{O}(u^4),$$

$$E_{\mathrm{sym}}(n)=S_0+\frac{L_0}{3}u+\frac{K_{\mathrm{sym}}}{18}{u}^2+\frac{Q_{\mathrm{sym}}}{162}{u}^3+\mathcal{O}(u^4),$$

$$P(n,Y_{p})=n^{2}\left\{\frac{\partial[E_{b}(n,Y_{p})+E_{e}(n,Y_{p})]}{\partial n}\right\}\biggr{\rvert}_{Y_{p},S},$$

$$\begin{gathered} P(n,Y_p)=\left(\frac{n^2}{3n_{\mathrm{sat}}}\right)\times \\ \left[\frac{K_0}{3}u+\frac{Q_0}{18}{u}^2+\left(L_0+\frac{K_{\mathrm{sym}}}{3}u+\frac{Q_{\mathrm{sym}}}{18}{u}^2\right)(1-2Y_p{)}^2\right]. \end{gathered}$$

$$\frac{Y_p}{(1-2Y_p{)}^3}=\frac{64E_{\mathrm{sym}}(n{)}^3}{3{\pi }^{2}n(\hslash c{)}^3}$$

$$(1-2Y_p{)}^2\approx a+bu+c{u}^2+\mathcal{O}(u^3),$$

$$\begin{gathered} P(n,Y_p)=\left(\frac{n^2}{3n_{\mathrm{sat}}}\right)\{a{L}_0+\left(b{L}_0+\frac{K_0+a{K}_{\mathrm{sym}}}{3}\right)u \\ +\left(c{L}_0+\frac{b{K}_{\mathrm{sym}}}{3}+\frac{Q_0+a{Q}_{\mathrm{sym}}}{18}\right)u^2\}. \end{gathered}$$

$$P(n)=K_{\mathrm{poly}}{n}^{\Gamma },$$

$$P(n)=\frac{a{L}_0{n}_{\mathrm{sat}}}{3}{\left(\frac{n}{n_{\mathrm{sat}}}\right)}^{\Gamma }.$$

$${\Lambda }_i=\frac{2}{3}{k}_2{\left(\frac{G{m}_i}{R_i{c}^2}\right)}^{-5},$$

$$\tilde{\Lambda }=\frac{16}{13}\frac{(m_1+12m_2)m_1^4{\Lambda }_1+(m_2+12m_1)m_2^4{\Lambda }_2}{(m_1+m_2{)}^5},$$

$${\mathcal{M}}_c=\frac{(m_1{m}_2{)}^{3/5}}{(m_1+m_2{)}^{1/5}}=m_1\frac{q^{3/5}}{(1+q{)}^{1/5}},$$

$$\mathcal{P}(L_0)=\mathcal{P}(\tilde{\Lambda })\left(\frac{\partial \tilde{\Lambda }}{\partial L_0}\right),$$

$$\begin{gathered} a{L}_0(u+1{)}^{\Gamma -2}=a{L}_0+\left(b{L}_0+\frac{K_0+a{K}_{\mathrm{sym}}}{3}\right)u \\ +\left(c{L}_0+\frac{b{K}_{\mathrm{sym}}}{3}+\frac{Q_0+a{Q}_{\mathrm{sym}}}{18}\right)u^2. \end{gathered}$$

$$(a-b)L_0-\frac{K_0+a{K}_{\mathrm{sym}}}{3}=0$$

$$c{L}_0+\frac{b{K}_{\mathrm{sym}}}{3}+\frac{Q_0+a{Q}_{\mathrm{sym}}}{18}=0.$$

$$(2a-b)L_0-\frac{K_0+a{K}_{\mathrm{sym}}}{3}=0$$

$$(c-a)L_0+\frac{b{K}_{\mathrm{sym}}}{3}+\frac{Q_0+a{Q}_{\mathrm{sym}}}{18}=0.$$

$$\left(\frac{3a}{2}-b\right)L_0-\frac{K_0+a{K}_{\mathrm{sym}}}{3}=0$$

$$\left(c-\frac{3a}{8}\right)L_0+\frac{b{K}_{\mathrm{sym}}}{3}+\frac{Q_0+a{Q}_{\mathrm{sym}}}{18}=0,$$

$$\Gamma =3.0:$$

$$\Gamma =3.5:$$

$$\Gamma =4.0:$$

$$\Gamma =3.0:$$

$$\Gamma =3.5:$$

$$\Gamma =4.0:$$

Full text

1. Introduction

Determining the nuclear symmetry energy is one of the main goals of modern nuclear physics. The symmetry energy, which characterizes the difference in energy between pure neutron matter and matter with equal numbers of protons and neutrons, is typically represented as a series expansion in density, with coefficients that represent the value of the symmetry energy at the nuclear saturation density, $S_{0}$, the slope, $L_{0}$, the curvature, $K_{\rm sym}$, and the skewness $Q_{\rm sym}$, as well as higher-order terms. The symmetry energy is one of the two main components in nuclear formulations of the dense-matter equation of state (EOS); the other being the energy of symmetric matter, which can similarly be broken down into nuclear expansion terms.

Of these expansion terms, only the low-order parameters can be experimentally constrained, as a result of the limited densities and energies that can be reached in laboratory-based experiments. For example, experimental constraints on $S_{0}$ and $L_{0}$ have been inferred by fitting nuclear masses, by measuring the neutron skin thickness, the giant dipole resonance, and electric dipole polarizability of 208Pb, and by observing isospin diffusion or multifragmentation in heavy ion collisions (Tsang et al., 2012; Lattimer & Lim, 2013; Oertel et al., 2017). However, there exist only limited experimental constraints on $K_{\rm sym}$ and no direct constraints on $Q_{\rm sym}$ (Lattimer & Lim, 2013).

The symmetry energy also plays a key role in a number of astrophysical phenomena, from determining the neutron star radius (Lattimer & Prakash, 2001), to affecting the gravitational wave emission during neutron star mergers (e.g., Fattoyev et al. 2013), r-process nucleosynthesis in merger ejecta (Nikolov et al., 2011), and the outcomes of core-collapse supernovae (e.g., Fischer et al. 2014). In the new gravitational wave era, measurements of the tidal deformability of neutron stars offer a promising way to observationally constrain the symmetry energy. The first detection of gravitational waves from a neutron star-neutron star merger, GW170817, constrained the effective tidal deformability of the binary system to $\widetilde{\Lambda}\lesssim 900$ (Abbott et al., 2017). Subsequent work refined these constraints to $\widetilde{\Lambda}=300\begin{subarray}{c}{}^{+420}_{-230}\end{subarray}$ (90% highest posterior density) for a system with chirp mass $\mathcal{M}_{c}=1.186\begin{subarray}{c}{}^{+0.001}_{-0.001}\end{subarray}~{}M_{\odot}$, for low-spin priors (Abbott et al., 2019).111Throughout this paper, we will exclusively use the low-spin prior results for GW170817, as is most relevant for binary neutron stars in our Galaxy.

Already, several analyses have set initial constraints on nuclear parameters using the tidal deformability of GW170817. Malik et al. (2018) found evidence of correlations between linear combinations of nuclear parameters and the neutron star radius, tidal love number, and tidal deformability, for a wide range of EOS. They used the inferred bounds on $\Lambda_{1.4}$ from GW170817, i.e., the tidal deformability of a 1.4 $M_{\odot}$ neutron star, to constrain the symmetric nuclear parameters as well as $K_{\rm sym}$, for given choices of $L_{0}$. Carson et al. (2019) expanded this work to include a broader set of EOS and used updated posteriors on $\widetilde{\Lambda}$ to calculate the posteriors for various nuclear parameters. In quoting final constraints on the high-order nuclear parameters, the authors of both studies either limit $L_{0}$ to a predetermined range or marginalize over priors on $L_{0}$. However, one might expect $\widetilde{\Lambda}$ to be particularly sensitive to $L_{0}$, as a result of the direct mapping between $\widetilde{\Lambda}$ and the neutron star radius (Raithel et al., 2018; De et al., 2018; Raithel, 2019) and the tight correlation between the radius and $L_{0}$ (Lattimer & Prakash, 2001).

In a more general analysis that allowed for variable $L_{0}$, Zhang & Li (2019) showed that a precision measurement of $\Lambda_{1.4}$ maps to a plane of constraints on $L_{0}$, $K_{\rm sym}$, and $Q_{\rm sym}$. They found that the tidal deformability is sensitive to the higher-order symmetry terms ($K_{\rm sym}$ and $Q_{\rm sym}$), and conclude that there is no unique mapping between $\Lambda_{1.4}$ and $L_{0}$. Krastev & Li (2018) extended this work and showed that the mapping gets even more complicated when the isotriplet/isosinglet interaction is allowed to vary. They found that assuming different density-dependences of the symmetry energy can result in identical values of the $\Lambda_{1.4}$, implying that there can be no one-to-one mapping between the tidal deformability and individual nuclear parameters.

While there may be no unique mapping of $\Lambda$ to individual nuclear parameters when the parameters are allowed to vary fully independently, we find that a more restricted parameter space is often sufficient to reproduce a wide range of EOS. With a well-motivated parameter reduction, we will show that it becomes possible to directly map from the tidal deformability to nuclear parameters.

In this paper, we introduce a framework to reduce the allowed space of nuclear parameters and we show that constraints can, indeed, be placed directly on the slope of the symmetry energy. Our method does not rely on priors from nuclear experiments. We only assume that the density-dependence of the EOS can be represented with a single-polytrope around the nuclear saturation density, which decreases the parameter space significantly. This dimension-reduction allows us to map from observed constraints on $\widetilde{\Lambda}$ directly to $L_{0}$, independently of any nuclear priors. We find that the tidal deformability is relatively insensitive to $S_{0}$, as has been assumed in the above analyses. However, we find that $\widetilde{\Lambda}$ is quite sensitive to $L_{0}$ and that GW170817 implies a relatively small value of $9.0\lesssim L_{0}\lesssim 65.4$ MeV, with a most likely value of $L_{0}\approx 22.5$ MeV. These constraints are approximate, but they point to values of $L_{0}$ that are significantly lower than those inferred from nuclear physics experiments and theory. Finally, we use the inferred posterior on $L_{0}$ to analytically constrain combinations of the higher-order nuclear terms. We find that combinations of $K_{\rm sym}$ and $K_{0}$ can be constrained by GW170817, and that combinations of $K_{\rm sym}$, $Q_{0}$, and $Q_{\rm sym}$ can also be constrained.

We start in $\S$2 with an overview of the nuclear EOS formalism that we will use in this paper. We introduce our polytropic approximation of these EOS in $\S$3. In $\S$4, we map the measurement of $\widetilde{\Lambda}$ from GW170817 to posteriors over $L_{0}$. Finally, in $\S$5, we use the posterior on $L_{0}$ to constrain linear combinations of the higher-order nuclear parameters.

2. Nuclear expansion of the equation of state

We start by introducing the EOS formalism that we will use to connect the tidal deformability from a gravitational wave event to nuclear parameters. As discussed in the introduction, we use this standard formalism to decompose the EOS into a symmetric matter part and the symmetry energy, which we can generically write as

[TABLE]

where $E_{b}(n,Y_{p})$ is the energy per baryon for a given density $n$ and proton fraction $Y_{p}$, $E_{0}(n)$ is the energy of symmetric matter, and $E_{\rm sym}(n)$ is the symmetry energy.

We represent the symmetric energy term with a series expansion and keep terms to third order, i.e.,

[TABLE]

where the expansion is performed around the nuclear saturation density, $n_{\rm sat}$, and $u\equiv(n/n_{\rm sat})-1$. Here, $B_{0}$ is the bulk binding energy of symmetric matter at $n_{\rm sat}$, and $K_{0}$ and $Q_{0}$ represent the incompressibility and skewness of symmetric matter.

Similarly, we expand the symmetry energy around $u$ and write

[TABLE]

where $S_{0}$ represents the symmetry energy at $n_{\rm sat}$ and $L_{0}$, $K_{\rm sym}$ and $Q_{\rm sym}$ give the slope, curvature, and skewness of the symmetry energy, respectively.

Such expansions are commonly used in representing neutron star matter because the coefficients can be linked to nuclear physics parameters near the saturation density. While experimental constraints on certain of these parameters exist, in order to be as general as possible, we will only assume knowledge of the bulk binding energy term and fix it to $B=-15.8$ MeV (Margueron et al., 2018a). We will leave the remaining six parameters ($K_{0},Q_{0},S_{0},L_{0},K_{\rm sym},Q_{\rm sym}$) free.

We can convert from the energy per particle to the pressure using the standard thermodynamic relation,

[TABLE]

where $S$ is the entropy, and we have formally included the electron contribution to the total energy, $E_{e}(n,Y_{p})$. However, for the current analysis, we neglect the contribution of electrons and assume that the total energy is dominated by the baryons.

The pressure for our nuclear expansion is then

[TABLE]

For a cold star in $\beta$-equilibrium, the proton fraction is uniquely determined by the density and the symmetry energy, according to

[TABLE]

where $\hbar$ is the Planck constant and $c$ is the speed of light. (For a derivation of this relationship and an analytic solution for $Y_{p}$, see Appendix A of Raithel et al. 2019).

In order to simplify the subsequent calculations, we perform an additional series expansion on the neutron excess parameter and define

[TABLE]

keeping terms up to second-order, as in eq. (5). The coefficients of this expansion depend on the symmetry energy parameters of up to the same order, i.e., $a=a(S_{0}),b=b(S_{0},L_{0})$, and $c=c(S_{0},L_{0},K_{\rm sym})$.222We provide a Mathematica notebook to calculate these coefficients, along with corresponding C routines, at https://github.com/craithel/Symmetry-Energy. Thus, keeping terms to second order, we can write the nuclear expansion of the pressure as

[TABLE]

3. Polytropic approximation

While the expansion derived in $\S$2 is useful for its direct connection to nuclear parameters, it is also complicated. The pressure of eq. (8) depends on 6 nuclear parameters: $S_{0}$, $K_{0}$, $Q_{0}$, $L_{0}$, $K_{\rm sym}$, and $Q_{\rm sym}$. However, many studies have shown that a wide range of EOS can be approximated with piecewise polytropic parametrizations (e.g., Read et al. 2009; Özel & Psaltis 2009; Steiner et al. 2010; Raithel et al. 2016). The pressure of a single-polytrope is given by

[TABLE]

where the polytropic constant $K_{\rm poly}$ and index $\Gamma$ are free parameters. The possibility of modeling the EOS with a few number of polytropes motivated us to explore whether the pressure in eq. (8) truly depends on all six nuclear parameters independently, or whether, as we will show, the parameter space can be further restricted. In this section, we will show that modeling the full pressure of eq. (8) with a single polytrope near the nuclear saturation density reasonably captures the density dependence. We will then use this simplified model to derive constraints on nuclear parameters using data from GW170817.

Our goal is to approximate the nuclear expansion pressure of eq. (8) with the polytropic pressure of eq. (9). We require that these two expressions match at $n_{\rm sat}$ and then extrapolate to higher densities using the polytropic index. This requirement uniquely determines the polytropic constant, so that our simplified nuclear pressure can be written as

[TABLE]

At low densities of $n\leq 0.5~{}n_{\rm sat}$, we fix the EOS to the nuclear EOS SLy (Douchin & Haensel, 2001). For $0.5~{}n_{\rm sat}\leq n<n_{\rm sat}$, we perform a power-law interpolation, to ensure matching between SLy and the polytropic approximation.

In order to test whether this simplified model of the pressure reasonably captures the density-dependence of the full nuclear expansion, we generate a sample of 1,000 test EOS using eq. (8). The EOS are created by drawing independent values of each of the six nuclear parameters ($K_{0},Q_{0},S_{0},L_{0},K_{\rm sym},$ and $Q_{\rm sym}$) from the experimentally-constrained distributions reported in Table I of Margueron et al. (2018b) (a similar approach was taken in Carson et al. 2019). We exclude any EOS that become hydrostatically unstable or that have superluminal sound speeds across a density range of $n\sim 0.01-10~{}n_{\rm sat}$. Additionally, we require that the analytic expression for $Y_{p}$ derived from eq. (6) be positive and less than 0.5 (i.e., neutron rich) across the same density range.

We fit each EOS in our sample with the simplified pressure model of eq. (10), fixing $L_{0}$ and $S_{0}$ to their drawn values. We perform the fit across the density range $n=1-3n_{\rm sat}$, in order to most strongly weight the density regime which is responsible for determining the neutron star radius, and hence the effective tidal deformability (Lattimer & Prakash, 2001; Raithel et al., 2018). We show the resulting distribution of $\Gamma$ values in Fig. 1. We find that nearly all of the EOS constructed with the nuclear expansion formalism can be represented with a polytropic index of $\Gamma=3-4$, with a most common fit value of $\Gamma\sim 3.5$. Moreover, we find that the residuals between the full nuclear expansion EOS and our polytropic approximation are small. Figure 2 shows the cumulative distribution of the residuals from the EOS fits at a range of densities. The residuals are smallest at low densities, where the tidal deformability is expected to be determined and where the nuclear expansion formalism still applies. At $2~{}n_{\rm sat}$, the error introduced by our polytropic approximation is $\lesssim$15% for 90% of the EOS sample. For completeness, we also show in Fig. 2 the residuals at core densities of $6-8~{}n_{\rm sat}$ (thin, dashed lines), even though the symmetry energy expansion is expected to break down at the these high densities. At $8~{}n_{\rm sat}$, the residuals of our approximation are still $\lesssim$ 50% for 90% of the sample.

We, therefore, find that the single-polytrope approximation reasonably recreates $P(n)$ for most combinations of the nuclear parameters. While this approximation is not exact, it is a useful technique that will allow us to explore the parameter-dependence of $\widetilde{\Lambda}$ in a new way. Our simplified model depends only on $S_{0}$ and $L_{0}$, thereby reducing a six-dimensional parameter space to two dimensions. This will allow us to directly map from $\widetilde{\Lambda}$ to $S_{0}$ and $L_{0}$, without requiring us to fix or marginalize over the higher-order terms.

4. Relating the tidal deformability to the symmetry energy

Using the framework for pressure introduced in $\S$3, we can now connect the observed constraints on $\widetilde{\Lambda}$ from a gravitational wave event to nuclear parameters. We start with the general expression for the tidal deformability of a single star,

[TABLE]

where $m_{i}$ is the mass of the star, $R_{i}$ is the stellar radius, and, following the convention of Flanagan & Hinderer (2008), we call $k_{2}$ the tidal apsidal constant. The tidal apsidal constant depends both on the compactness of the star, as well as the overall density gradient of the particular EOS (Hinderer, 2008; Hinderer et al., 2010; Postnikov et al., 2010).

We follow the method outlined in Hinderer et al. (2010) for constructing a set of augmented Oppenheimer-Volkoff equations. We integrate these stellar structure equations to calculate the stellar mass, radius, and tidal apsidal constant for a given central density. We then compute the effective tidal deformability of the binary system as

[TABLE]

where the subscripts indicate the component stars in the binary system.

We show the effect of $L_{0}$ on each of these stellar properties in Fig. 3. For demonstrative purposes, in this figure we have fixed $S_{0}$=32 MeV and $\Gamma=3.5$, as well as the component masses for calculating $\widetilde{\Lambda}$ (see below for the effect of varying each of these assumptions). We show the mass-radius relations for a variety of $L_{0}$ in the left panel of Fig. 3. The middle panel panel shows the tidal apsidal constant as a function of the stellar compactness, for the same set of $L_{0}$ values. We find that smaller values $L_{0}$ lead to both smaller radii and to smaller tidal apsidal constants. Both of these trends act to reduce the tidal deformability of the star (see eq. 11), as shown in the right panel of Fig. 3. We find that the dependence on $L_{0}$ persists both for $\Lambda_{1.4}$ and for the binary tidal deformability, $\widetilde{\Lambda}$. Thus, we expect that the measurement of $\widetilde{\Lambda}$ from a gravitational wave event should have significant constraining power on $L_{0}$.

From eq. (12) and Fig. 3, it is clear that $\widetilde{\Lambda}$ depends on the stellar masses and radii, with an additional dependence on the EOS through $k_{2}$. By introducing the chirp mass,

[TABLE]

and the mass ratio $q\equiv m_{2}/m_{1}$, we can explicitly write the dependences of $\widetilde{\Lambda}$ as $\widetilde{\Lambda}(\mathcal{M}_{c},q,R_{1},R_{2},\rm{EOS})$. This is a particularly convenient choice because, for gravitational wave events, we expect the chirp mass to be precisely measured and the mass ratio to also be constrained. This was indeed the case for GW170817, for which the the chirp mass was determined to be $\mathcal{M}_{c}=1.186\begin{subarray}{c}+0.001\\ -0.001\end{subarray}~{}M_{\odot}$ and the mass ratio was constrained to $q\in(0.73,1.00)$ at the 90% confidence level (Abbott et al., 2019).

Additionally, we note that, given the masses of each star, the EOS can be used to uniquely determine the corresponding radii.333While there do exist some EOS for which the mass-to-radius mapping is not unique (notably, the so-called “twin-stars,” which can have identical masses and different radii; see, e.g., Glendenning & Kettner 2000), these EOS have complex structure that cannot be represented with single polytropes. We, therefore, neglect these special cases for the present study. Within the polytropic approximation of eq. (10), the EOS depends only on $S_{0}$, $L_{0}$, and $\Gamma$, where $\Gamma$ is narrowly constrained to be $\sim 3-4$ for a wide range of realistic EOS. We can, therefore, summarize the dependences of the tidal deformability as $\widetilde{\Lambda}=\widetilde{\Lambda}(\mathcal{M}_{c},q,S_{0},L_{0},\Gamma)$.

We have already shown that $\widetilde{\Lambda}$ depends sensitively on $L_{0}$ for fixed $S_{0}$ and $\Gamma$. In order to explore the full, more general dependences of $\widetilde{\Lambda}$, we perform a grid search across the range $S_{0}\in(26,38)$ MeV and $L_{0}\in(10,120)$ MeV. For each set of values, we construct an EOS according to eq. (10), fixing the mass ratio to $q=0.7,0.87,$ or 1.0 and fixing $\Gamma$ to 3, 3.5, or 4. In all cases, we fix the chirp mass to the central value from GW170817 of $1.186~{}M_{\odot}$. For each combination of parameters, we compute the mass, radius, and tidal apsidal constant by numerically integrating the augmented TOV equations and then compute $\widetilde{\Lambda}$ using eqs. (11-13). We show the resulting contours of $\widetilde{\Lambda}$ as a function of $S_{0}$ and $L_{0}$ in Fig. 4. The three panels correspond to three different choices of $\Gamma$, with fixed $q=0.87$. We find that the particular choice of $q$ does not significantly affect these or our later results, so we fix $q$ to the central value of 0.87 from GW170817 for the remainder of this analysis.

We find that $\widetilde{\Lambda}$ is only weakly dependent on $S_{0}$, especially for smaller values of $\widetilde{\Lambda}$, as are preferred by the current gravitational wave data. In contrast, $\widetilde{\Lambda}$ depends quite sensitively on $L_{0}$. We, therefore, focus on $L_{0}$ in the following analysis and fix $S_{0}$ to a characteristic value of 32 MeV (Li & Han, 2013; Oertel et al., 2017).

This final simplification renders $\widetilde{\Lambda}$ as a function only of $L_{0}$, for fixed $\Gamma$. We can, therefore, transform the measured posterior on $\widetilde{\Lambda}$ to a posterior on $L_{0}$, according to

[TABLE]

where we calculate the Jacobian term numerically.

We show the resulting one-dimensional posteriors on $L_{0}$ in Fig. 5. Figure 5 also shows two current sets of constraints on $L_{0}$, in dark and light green from Lattimer & Lim (2013) and Oertel et al. (2017), respectively, which are based on a combination of astrophysical observations of neutron stars, nuclear experiments, and theory. Earlier constraints on $43<L_{0}<52$ MeV (68% confidence) were calculated using neutron star radii alone (Steiner & Gandolfi, 2012). On the other hand, theoretical calculations of the neutron matter EOS using quantum Monte Carlo methods (Gandolfi et al., 2012) or chiral effective field theory (Hebeler et al., 2013) produce comparable constraints, of $L_{0}=31.3-63.6$ MeV and $L_{0}=32.4-57.0$ MeV, respectively. We use the summary results from Lattimer & Lim (2013) and Oertel et al. (2017) to encompass these theoretical, observational, and experimental constraints.

We find that the gravitational wave data imply smaller values of $L_{0}$ than these previous studies have found. In particular, for $\Gamma=3.5$, we find a 90% highest-posterior density interval of $9.0<L_{0}<65.4$ MeV, with a peak likelihood at $L_{0}\approx 22.5$ MeV. There is a small correlation between the choice of $\Gamma$ and the inferred constraints on $L_{0}$, with choices of larger values of $\Gamma$ leading to lower values of $L_{0}$. Nevertheless, for $\Gamma=3.5$ or 4, the peak likelihoods in $L_{0}$ lie outside the allowed constraints from both Lattimer & Lim (2013) and Oertel et al. (2017). For $\Gamma=3$, the peak likelihood falls at the lower limit of the constraint from Lattimer & Lim (2013).

Several recent studies connecting GW170817 to the nuclear EOS have either restricted $L_{0}$ to similar priors or marginalized over them (Malik et al., 2018; Carson et al., 2019). While the posteriors on $L_{0}$ presented here are not exact, they do suggest that GW170817 points toward small values of $L_{0}$ that may be in tension with such priors. We, therefore, conclude that it is important to explore the dependence of gravitational wave data on $L_{0}$ directly, in order to gain new information on $L_{0}$ itself as well as to avoid biasing the interpretation of higher-order parameters.

5. Constraints on higher-order nuclear parameters

In $\S$4, we showed that GW170817 directly maps to constraints on $L_{0}$ using our polytropic approximation of the nuclear expansion. In this section, we turn to the higher-order nuclear terms. In particular, we will show that by taking the polytropic approximation in eq. (10), we can place constraints on the allowed combinations of the remaining four nuclear parameters.

We start by equating the polytropic approximation and the full nuclear expansion, and match terms of equivalent order, i.e.,

[TABLE]

For this expression to be true at all densities, the terms of equivalent order must all sum to zero. For example, setting $\Gamma=3$ implies the constraints

[TABLE]

Likewise, setting $\Gamma=4$ implies

[TABLE]

For $\Gamma=3.5$, we introduce one final series expansion on the left-hand side of eq. (15) to simplify $(u+1)^{1.5}\approx 1+(3/2)u+(3/8)u^{2}+\mathcal{O}(u^{3})$. Using this approximation, we can again require terms of the same order to sum to zero, and we find

[TABLE]

where we recall that $a=a(S_{0}),b=b(S_{0},L_{0}),$ and $c=c(S_{0},L_{0},K_{\rm sym})$.

Thus, for a given choice of $\Gamma$, we have two independent sets of constraints: the first connects the parameter set $\{S_{0}$, $L_{0}$, $K_{0}$, and $K_{\rm sym}\}$, while the second connects $\{S_{0}$, $L_{0}$, $K_{\rm sym}$, $Q_{0}$, $Q_{\rm sym}\}$. In the following, we will use the posterior on $L_{0}$ from $\S$4 to constrain the remaining combinations of higher-order terms using these relationships.

We start with the first set of constraints, on $\{S_{0}$, $L_{0}$, $K_{0}$, and $K_{\rm sym}\}$. These constraints correspond to eqs. (16a), (17a), and (18a). As in $\S$ 4, we fix $S_{0}$=32 MeV and find that this choice does not strongly affect the results. We then use the 90%-credible interval on $L_{0}$ from Fig. 5 to bound the allowed range of $K_{\rm sym}-K_{0}$ values. We show the resulting constraints in Fig. 6 for each $\Gamma$. In this figure, the light shaded regions represent the bounds on the $K_{\rm sym}-K_{0}$ relationship allowed by the 90%-credible interval on $L_{0}$. The dark solid lines indicate the $K_{\rm sym}-K_{0}$ relationship corresponding to the most likely value of $L_{0}$, for a given $\Gamma$.

We find that GW170817 places tight constraints on linear combinations of $K_{0}$ and $K_{\rm sym}$. For the three values of $\Gamma$, fixing $L_{0}$ to its maximum likelihood value from Fig. 5 yields

[TABLE]

These equations correspond to the dark, solid lines in Fig. 6. The coefficient in front of $K_{0}$ is the same in all three cases because it is given simply by $1/a$ and hence depends only on $S_{0}$. The constant term depends on $S_{0}$, $L_{0}$, and the coefficients of $L_{0}$ in eqs. (16a), (17a), and (18a), and thus varies slightly with the choice of $\Gamma$.

Previous studies have constrained $K_{0}$ by fitting nuclear models to measurements of the isoscalar giant monopole resonance. Depending on the analysis methods, the results range from quite narrow, $K_{0}=248\pm 8$ MeV (Piekarewicz, 2004) and $K_{0}=240\pm 20$ MeV (Shlomo et al., 2006), to broader bounds of $K_{0}=250-315$ MeV (Stone et al., 2014). If we take the broadest range of these allowed values and assume $220<K_{0}<315$ and combine this with our results in Fig. 6, we find that $-375\lesssim K_{\rm sym}\lesssim 45$, at 90% confidence. This constraint on $K_{\rm sym}$ is broader than, but consistent with previous results, including the constraint of $K_{\rm sym}=-111.8\pm 71.3$ MeV that was derived from universal relations between $K_{\rm sym}$ and lower-order expansion terms (Mondal et al., 2017).

Using the tidal deformability from GW170817, Malik et al. (2018) found constraints of $-112<K_{\rm sym}<-52$ MeV or $-140<K_{\rm sym}<16$ MeV, depending on their choice of prior for $L_{0}$. In a similar analysis, Carson et al. (2019) derived constraints of $-259\leq K_{\rm sym}\leq 32$ MeV, after marginalizing over $L_{0}$. Our results, which are derived with the polytropic approximation and no priors on $L_{0}$, are consistent with both of these analyses. However, if we take the $K_{\rm sym}(K_{0})$ relationship that corresponds to the maximum likelihood in $L_{0}$ (i.e., the dark solid lines in Fig. 6), we find that the data point to smaller values of $K_{\rm sym}$, below the lower bound from the Malik et al. (2018) study and on the lower end of the Carson et al. (2019) constraints. This is likely a consequence of the fact that these studies both used priors that either forbade or disfavored low values of $L_{0}$, such as those we find in this paper.

Finally, we turn to the second set of constraints on the parameters $\{S_{0}$, $L_{0}$, $K_{\rm sym}$, $Q_{0}$, $Q_{\rm sym}\}$. As we did above, we will fix $S_{0}=32$ MeV and use the maximum likelihood in $L_{0}$. We can then use eqs. (16b), (17b), and (18b) to calculate the relationship between the remaining four parameters. For the most likely value of $L_{0}$, we find

[TABLE]

To our knowledge, no nuclear experiments have constrained $Q_{0}$ or $Q_{\rm sym}$ and only broad theoretical bounds have been calculated. For example, Zhang et al. (2017) found $-800<Q_{0}<400$ MeV based on analyses of energy density functionals. Nevertheless, future experiments or astrophysical observations may one day provide stricter bounds on $Q_{0}$ or $Q_{\rm sym}$. Within our polytropic framework, any such measurements can then be used to constrain the correlated parameters using these analytic relationships.

6. Conclusions

In this paper, we have introduced a new approximation of the nuclear EOS which allows for a direct mapping from measured $\widetilde{\Lambda}$ constraints to the symmetry energy parameters. We have shown that a wide sample of nuclear EOS can be reasonably represented with a single-polytrope approximation in the density range of interest, which simplifies the EOS to depend only on $S_{0}$, $L_{0}$, and $\Gamma$, rather than the full six nuclear parameters. Moreover, we find that $\widetilde{\Lambda}$ is relatively insensitive to $S_{0}$. With many future gravitational wave detections expected in the coming years, this framework will make it possible to map the gravitational wave event directly to $L_{0}$ or to combinations of higher-order nuclear parameters.

With this parameter-space reduction and focusing on the existing measurement of $\widetilde{\Lambda}$ from GW170817, we were able to map from the full posterior on $\widetilde{\Lambda}=300\begin{subarray}{c}{}^{+420}_{-230}\end{subarray}$ (Abbott et al., 2019) to posteriors on $L_{0}$. We find that GW170817 points to significantly smaller values of $L_{0}$ than have been previously been reported, with a peak likelihood of $L_{0}\sim 20$ MeV. We additionally use these posteriors on $L_{0}$ to constrain combinations of higher-order nuclear parameters, finding tight constraints on the allowed combinations of $K_{\rm sym}$ and $K_{0}$, as well as constraints on $K_{\rm sym},Q_{0},$ and $Q_{\rm sym}$.

We note that the final constraints on $L_{0}$ depend slightly on the choice of $\Gamma$ and, of course, will depend on the robustness of our polytropic approximation. If the true combination of $K_{0},Q_{0},S_{0},L_{0},K_{\rm sym},$ and $Q_{\rm sym}$ produce an EOS with significant sub-structure, then our single-polytrope approximation is not the optimal approach. Moreover, if the dense-matter EOS contains a phase transition to quark matter, then the polytropic approximation will be inadequate and, depending on the particular formulation, the relationship between the tidal deformability and $L_{0}$ may be significantly weaker as well (e.g., Zhu et al. 2018).

Nevertheless, the results in this paper indicate that gravitational wave data can significantly constrain the slope of the symmetry energy for nuclear EOS. This is an important point. Previous studies connecting GW170817 to the nuclear EOS have either fixed the allowed range of $L_{0}$ or marginalized over $L_{0}$ (Malik et al., 2018; Carson et al., 2019), using priors from nuclear physics that we find to be in modest conflict with the values inferred from GW170817. As the LIGO/Virgo team continue to observe new gravitational wave events and further pin down the tidal deformability of neutron stars, it will become increasingly important to develop robust approaches to constrain the nuclear parameters in model-independent ways.

Acknowledgments. We thank Dimitrios Psaltis, Kent Yagi, Andrew Steiner, and Zack Carson for useful conversations related to this work. This work is supported by NSF Graduate Research Fellowship Program Grant DGE-1746060 and support from NASA grant NNX16AC56G.