# Resonant Shattering Flares as Multimessenger Probes of the Nuclear Symmetry
Energy

Duncan Neill,1 William G. Newton,2 David Tsang1
1Department of Physics, University of Bath, Claverton Down, Bath, BA1 1AL
2Department of Physics and Astronomy, Texas A&M University-Commerce, Commerce,
TX, 75429-3011 E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The behaviour of the nuclear symmetry energy near saturation density is
important for our understanding of dense nuclear matter. This density
dependence can be parameterised by the nuclear symmetry energy and its
derivatives evaluated at nuclear saturation density. In this work we show that
the core-crust interface mode of a neutron star is sensitive to these
parameters, through the (density-weighted) shear-speed within the crust, which
is in turn dependent on the symmetry energy profile of dense matter. We
calculate the frequency at which the neutron star quadrupole ($\ell=2$) crust-
core interface mode must be driven by the tidal field of its binary partner to
trigger a Resonant Shattering Flare (RSF). We demonstrate that coincident
multimessenger timing of an RSF and gravitational wave chirp from a neutron
star merger would enable us to place constraints on the symmetry energy
parameters that are competitive with those from current nuclear experiments.

###### keywords:

dense matter – stars: neutron – stars: oscillations – neutron star mergers –
gravitational waves – equation of state

††pubyear: 2015††pagerange: Resonant Shattering Flares as Multimessenger
Probes of the Nuclear Symmetry Energy–Resonant Shattering Flares as
Multimessenger Probes of the Nuclear Symmetry Energy

##  1 Introduction

Neutron stars contain the most extreme matter in the universe. They act as
natural laboratories to investigate nuclear physics, allowing us to study the
physics of matter at nuclear densities. To investigate the internal structure
of these compact stars we must probe them using observational phenomena.

Short gamma-ray bursts (SGRBs) (Kouveliotou et al., 1993; D’Avanzo, 2015)
likely originate from the merging of binary neutron stars (Eichler et al.,
1989; Fong et al., 2010). These bursts are characterised by a large peak in
the gamma-ray count that lasts for $\sim 2$ second or less. Around $\sim
3-10$% of SGRBs are preceded by a ‘precursor’ flare (Zhong et al., 2019; Troja
et al., 2010). Precursor flares can be identified as a lower, separate peak in
the gamma-ray count a short time ($\sim 0.1-5.0$ s) before the main peak
(Zhong et al., 2019). Resonant Shattering Flares (Tsang et al., 2012b) are
relatively isotropic, short ($\sim 0.1$s duration) gamma-ray flares that are
triggered by a tidal resonance of the binary, and can appear as either
precursor flares, or orphan flares if the main SGRB is beamed away from the
observer.

Recent observations of neutron star merger GW170817 (Abbott et al., 2017;
Goldstein et al., 2017) have begun a new era of multi-messenger astronomy
involving gravitational waves and counterparts across the electromagnetic
spectrum, allowing an unprecedented probe into the physics of binary neutron
star mergers (see e.g. Raithel, 2019, and references therein). Here, we will
show how multi-messenger coincident timing of a gravitational-wave chirp and
the prompt-gamma ray emission from a Resonant Shattering Flare can be used to
determine the frequency of a particular neutron star oscillation mode and
hence constrain fundamental parameters in nuclear physics.

###  1.1 Resonant Shattering Flares

During the gravitational-wave induced inspiral of neutron star binaries, the
normal modes of a neutron star can become excited by resonant tidal
interactions with its binary partner. As the binary orbit shrinks, the
frequency of the orbit (and hence the gravitational wave frequency) increases.
When the orbital frequency sweeps through the appropriate resonance windows,
the normal modes are excited, causing their oscillations to rapidly grow in
magnitude (Tsang et al., 2012b; Tsang, 2013; Lai, 1994).

If a resonant mode is sufficiently large to cause a deformation that exceeds
the breaking strain in the neutron star crust, the crust will fracture,
depositing seismic energy into the crust. One such mode is the quadrupole
crust-core interface ($i$) mode. This mode is caused by the discontinuity in
bulk material properties between the crust and core, and has a sufficient
overlap with the tidal field (Tsang et al., 2012b) for the resonance to
quickly deposit energy into the mode. The fractures and seismic waves continue
to be driven until the crust reaches its elastic limit, causing it to shatter.
This scatters the seismic waves to high frequencies that are able to couple to
the star’s magnetic field, depositing energy into the magnetosphere and
sparking a pair-photon fireball. Multiple colliding shells may be emitted over
the course of the resonance window ($t_{\rm res}\sim 0.1$ seconds), leading to
a single non-thermal burst with a duration of the resonance time, and maximum
luminosity determined by the surface magnetic field strength (Tsang, 2013).

For sufficiently strong surface magnetic field, these Resonant Shattering
Flares (RSFs) can be detectable well beyond the Advanced LIGO horizon.
Coincident timing of the RSF prompt emission and the gravitational-wave chirp
provides a precise measurement of the resonant mode frequency. Therefore, by
using a model of the crust and core with consistent nuclear physics to
calculate the $i$-mode frequency, the parameters of that model could be
restricted to ranges which results in mode frequencies that closely match the
observed gravitational wave frequency during the RSF.

In this paper, we will explore the implications of a coincident multi-
messenger detection of an RSF on constraining fundamental nuclear physics
parameters. In particular we will explore the relationship between the
$i$-mode frequency, neutron star structure, and the nuclear symmetry energy
parameters that determines the behaviour of matter near nuclear saturation
density.

###  1.2 Nuclear Symmetry Energy and Neutron Star Structure

Normal mode frequencies are dependent on the neutron star equation of state
(EOS) (the relationship between the energy density and pressure within the
star) and the composition of the crust as a function of density. At low
densities the EOS can be accurately calculated by using the properties of
experimentally measured nuclei (Baym et al., 1971). However, at the extreme
densities found in neutron stars we must rely on nuclear models to calculate
the EOS, with different models giving significantly different EOSs due to the
relative lack of experimental data that probes neutron-rich matter. This
uncertainty is conveniently captured by the nuclear symmetry energy $E_{\rm
sym}$, which encodes the average energy per nucleon required to convert
symmetric nuclear matter (number of protons $Z$ = number of neutrons $N$) into
pure neutron matter ($Z$=0). At nuclear saturation density (the density of
nucleons in a nucleus) the symmetry energy can be parameterised by the
coefficients of its density expansion, the first three of which are: the
magnitude of the symmetry energy at saturation density

$\displaystyle J=E_{\rm sym}(n_{0}),$ (1)

the slope of the symmetry energy

$\displaystyle L=3n_{0}\frac{\partial E_{\rm sym}(n_{\rm b})}{\partial n_{\rm
b}}\biggr{\rvert}_{n_{\rm b}=n_{0}},$ (2)

and the curvature

$\displaystyle K_{\rm sym}=9n_{0}^{2}\frac{\partial^{2}E_{\rm sym}(n_{\rm
b})}{\partial n_{\rm b}^{2}}\biggr{\rvert}_{n_{\rm b}=n_{0}},$ (3)

where $n_{0}\approx 0.16$ fm-3 is nuclear saturation density and $n_{\rm b}$
is the baryon density.

The symmetry energy parameters have some simple physical implications. The
magnitude of the symmetry energy $J$ controls the proton fraction at
saturation density, and the slope $L$ correlates with the pressure at that
density. $K_{\rm sym}$ controls the derivative of pressure, thus determining
how the pressure changes as one moves away from saturation density, and plays
an important role determining the stability of matter and its compressibility.
Through these effects, neutron star structure and composition are sensitive to
the symmetry energy and related nuclear observables (Brown, 2000; Horowitz &
Piekarewicz, 2001; Steiner et al., 2008; Fattoyev et al., 2018), and therefore
astrophysical observations provide information on $E_{\rm sym}$ that
complement nuclear experiment. Constraining the symmetry energy has been a
priority in the field of nuclear physics over the past two decades (Tsang et
al., 2009; Li et al., 2014; Horowitz et al., 2014), and the burgeoning field
of multimessenger astronomy provides exciting opportunities to synthesise
astrophysical observation and nuclear experiment to learn more about dense
matter (Tsang et al., 2019). Crust properties and their observational
manifestations are particularly sensitive to the symmetry energy (Newton et
al., 2014), including shear waves in the dense regions of the inner crust
(Steiner & Watts, 2009; Sotani et al., 2012, 2013; Gearheart et al., 2011).
RSFs are therefore a strong candidate to probe the symmetry energy (Tsang et
al., 2012b; Newton et al., 2014).

##  2 Parameterised Neutron Star Equation of State and Composition

###  2.1 The Nuclear model

A consistent nuclear physics description of a neutron star requires an
underlying nuclear model. We use an extended Skyrme mean field model for the
uniform matter EOS (Holt & Lim, 2018). Three parameters of the model affect
only the pure neutron matter EOS, and can be written in terms of $J$, $L$ and
$K_{\rm sym}$ (Newton & Crocombe, 2020; Balliet et al., 2020). This way, a
given set of values {$J$, $L$, $K_{\rm sym}$} can be converted into a set of
Skyrme models (each of which labeled by the 3 symmetry energy parameters) and
corresponding crust models and equations of state.

We explore the distributions of $i$-mode frequencies predicted for three
different sets of symmetry energy parameters.

Figure 1: The ranges used for the first two symmetry energy parameters $J$ and
$L$. The ranges covered by the Uniform, MSL0 and PNM parameter sets are shown
in grey, red and blue respectively. The Uniform and MSL0 sets are truncated in
the high $L$, low $J$ region because stable neutron star crust models do not
exist in that region of parameter space. Also pictured is the region in the
intersection of a number of nuclear experimental (Kortelainen et al., 2010;
Chen et al., 2010; Piekarewicz et al., 2012; Tsang et al., 2009), theoretical
(Hebeler et al., 2010; Gandolfi et al., 2012) and astrophysical (Steiner et
al., 2010, 2013) constraints, as synthesized by Lattimer & Lim (2013).

(i) Our most conservative range of EOSs uses symmetry energy parameters
uniformly distributed over the ranges $25\leq J\leq 35$ MeV, $20\leq L\leq 80$
MeV, and $-200\leq K_{\rm sym}\leq 40$ MeV. These ranges were chosen to cover
the ranges inferred from a variety of experimental probes (Liu et al., 2010;
Tsang et al., 2012a; Lattimer & Lim, 2013). This will be referred to as our
“Uniform” distribution, and is shown in the shaded blue region of Figure 1 in
$J$-$L$ space. Note that this set, as well as the MSL0 set below, are
truncated in the high$L$-low$J$ region of parameter space as stable crust
models do not exist for those combinations of parameters, as explained in
section 2.2.

(ii) The properties of pure neutron matter are the most important ingredient
in the crust EOS, so relevant constraints for neutron stars comes from pure
neutron matter theory. Great strides have been made in modeling pure neutron
matter from first principles, particularly using chiral effective field
theory. A useful way to parameterise these models is through a Taylor
expansion of the Fermi liquid parameters that characterize the two-neutron
interaction energy Holt & Kaiser (2017); Holt & Lim (2018), three of which are
shown to be sufficient to paramterise range of predictions for the PNM EOS
from $ab-initio$ calculations. We translate this range into $J$, $L$ and
$K_{\rm sym}$ space as detailed in Newton & Crocombe (2020) using a $9\times
9\times 9$ grid uniform in the Fermi liquid parameter space, and use this to
form our second set of neutron star models. In $J$-$L$ space, this range is
shown as the blue shaded region in Figure 1. This will be referred to as our
“PNM” distribution.

(iii) In many Skyrme models, such as Sly4 (Chabanat et al., 1998) and SkI6
(Reinhard & Flocard, 1995; Nazarewicz et al., 1996), there are only two
parameters that control the PNM EOS alone, and therefore $K_{\rm sym}$ is not
a free parameter. Likewise, only $J$ and $L$ are independent in many
extractions of symmetry energy constraints. It will be useful to examine how
important it is to take into account the $K_{\rm sym}$ degree of freedom, so
for our third parameter set we emulate Skyrme models that have only $J$ and
$L$ as free parameters. We take the dependence of $K_{\rm sym}$ from the MSL0
parameterisation of the Skyrme model (Chen et al., 2009), which has previously
been used to extract symmetry energy constraints from nuclear experiment. In
that model, $K_{\rm sym}$ is related to $J$ and $L$ by (Newton & Crocombe,
2020)

$\displaystyle K_{\rm sym}=3.71L-11.13J+11.93\;\rm{MeV},$ (4)

which restricts us to a single plane in the $J$,$L$,$K_{\rm sym}$ parameter
space. We will refer to this as our “MSL0” distribution. Note that there is
nothing physically special about this particular choice of relation between
the symmetry energy parameters.

We show the regions of $J$-$L$ space we sample for each of the three
distributions in Figure 1. The regions we cover encompass a number of nuclear
experimental constraints, the intersection of which is also shown as the white
hatched region Lattimer & Lim (2013).

###  2.2 The crust model

To calculate the crust composition and EOS, we use our sets of extended Skyrme
models in a compressible liquid drop model (CLDM) (Newton et al., 2013;
Balliet et al., 2020). The model assumes a lattice consisting a single species
of nucleus immersed in a neutron gas in a repeating unit cell (the Wigner-
Seitz approximation). By minimising the energy of the unit cell with respect
to the physical parameters of the cell - the neutron gas density $n_{\rm n}$,
the cell radius $r_{\rm c}$, and the mass and charge number of the nuclear
cluster $A,Z$ \- we obtain the ground state composition and EOS of the crust
at a given density. We can then calculate quantities required to model the
normal modes in the crust, for example the shear modulus and frozen-
composition adiabatic index as a function of baryon density $n_{\rm b}$
(Strohmayer et al., 1991; Chugunov & Horowitz, 2010).

The shear modulus in the crust is given by

$\displaystyle\mu=\frac{0.1106}{1+17810\left(\frac{ak_{b}T}{\left(Ze\right)^{2}}\right)^{2}}\frac{n_{i}\left(Ze\right)^{2}}{a},$
(5)

where $T$ is the temperature, $n_{i}$ is the ion density, $Z$ is the proton
number of the nuclei, and

$\displaystyle a=\left(\frac{3}{4\pi n_{i}}\right)^{\frac{1}{3}}.$ (6)

We conduct our calculations of the neutron star modes in the zero temperature
limit. The ion number density can be written in terms of the fraction of
nucleons in the neutron gas $X_{\rm n}$ through the relation (Newton et al.,
2013)

$X_{\rm N}=1-X_{\rm n}=\frac{n_{\rm i}}{n_{b}}A,$ (7)

where $X_{\rm N}$ is the fraction of nucleons in the nucleus. We can therefore
re-write the shear modulus as a function of $X_{\rm n}$ (at zero temperature)
as (Steiner et al., 2008; Gearheart et al., 2011)

$\mu=0.1106\left(\frac{4\pi}{3}\right)^{1/3}A^{-4/3}n_{\rm b}^{4/3}(1-X_{\rm
n})^{4/3}(Ze)^{2}.$ (8)

The adiabatic index at constant composition is given by

$\Gamma_{\rm 1}=\frac{n_{\rm b}}{P}\frac{dP}{dn_{\rm b}}\bigg{|}_{\rm
constant\;composition}=\frac{n_{\rm b}}{P}\bigg{[}\frac{dP_{n}}{dn_{\rm
n}}+x\frac{dP_{e}}{dn_{\rm e}}\bigg{]},$ (9)

where $x$ is the average proton fraction and $P$ the total pressure. $P_{\rm
n,e}$ are the pressures and $n_{\rm n,e}$ the number densities of dripped
neutrons and electrons respectively.

We calculate crust models over the full range of each of our three symmetry
energy ranges. However, EOSs with low values of $J$ and high values of $L$ do
not result in stable crust models. Such EOSs have a symmetry energy that falls
rapidly with decreasing density. If the magnitude of symmetry energy at
saturation density is already small, then at sub-saturation densities the
slope of the symmetry energy - and hence the pressure of pure neutron matter -
must become very small, or even become negative. Since neutron matter pressure
supports the inner crust in hydrostatic equilibrium, such crust models will be
inherently unstable. This is the reason our ranges of $J$ and $L$ in Figure 1
are truncated at high $L$, low $J$.

###  2.3 The core model

The Skyrme model is designed to describe nuclear interactions around nuclear
saturation density. As one moves into the neutron star core, the increasing
importance of relativistic effects, the possible appearance of hyperons at
supersaturation density, and the likely transition from nucleonic to quark
degrees of freedom in the inner core mean that the model is not suited to
describe matter beyond about twice saturation density, and the symmetry energy
loses its physical meaning. In order to explore the symmetry energy effects on
RSFs, we do need to control the core EOS. We use the piecewise polytrope
method (Read et al., 2009a, b; Steiner et al., 2010, 2013; Özel & Psaltis,
2009; Özel et al., 2010; Özel et al., 2016): we fit two polytropes at
supersaturation densities, one at a density of $n_{1}$=1.5$n_{0}$ and one at
$n_{2}$=2.7$n_{0}$, as detailed in Newton et al. (2018). We then have three
regions of the star: the crust and outer core, in which the pressure and
energy density are given by the Skyrme EOS, and the two polytropic regions in
which the pressures are given by

$\displaystyle P=P_{\rm Skyrme}$ $\displaystyle\;\;\;\;\;n<n_{1}$
$\displaystyle P_{1}=K_{1}n^{\gamma_{1}}$
$\displaystyle\;\;\;\;\;n_{1}<n<n_{2}$ $\displaystyle
P_{2}=K_{2}n^{\gamma_{2}}$ $\displaystyle\;\;\;\;\;n_{2}<n$ (10)

where continuity of pressure determines the constants $K_{1}$ and $K_{2}$. The
energy density in the three density regions is obtained by integrating the
first law of thermodynamics:

$\epsilon_{i}=(1+a_{i})n+\frac{K_{i}n^{\gamma_{i}}}{\gamma_{1}-1};\;\;\;\;\;\;\;a_{i}=\frac{\epsilon_{i-1}(n_{i})}{n_{i}}-\frac{K_{i}n^{\gamma_{i}-1}}{\gamma_{i}-1}-1$
(11)

where $a_{i}$ are constants of integration, $i$={1,2} and the subscript 0
labels the Skyrme EOS.

The speed of sound is

$\frac{c_{\rm
s,i}(n)}{c}=\bigg{(}\frac{\gamma_{i}P}{P+\epsilon}\bigg{)}^{1/2}$ (12)

In the eventuality that the EOS becomes acausal ($c_{\rm s}>c$) at a given
density $n_{\rm acausal}$, we transition to a causal EOS:

$P_{\rm causal}=\epsilon=bn^{1/3}\;\;\;\;\;n_{\rm acausal}<n$ (13)

where $\epsilon$ is the energy density, $b$ is a constant given by

$b=\frac{1+a}{n_{\rm acausal}}+K\frac{n_{\rm acausal}^{\gamma-2}}{\gamma-1}$
(14)

and $a$ is either $a_{1}$ or $a_{2}$ depending on which region the EOS becomes
acausal in.

Each equation of state we generate is characterised by 5 parameters: the three
symmetry energy coefficients $J$, $L$ and $K_{\rm sym}$ for the Skyrme-EOS,
and the polytropic parameters $\gamma_{1}$ and $\gamma_{2}$. $\gamma_{2}$,
which controls the high density part of the EOS, can be tuned to give a
desired maximum mass. $\gamma_{1}$, which controls the EOS at intermediate
densities in the core, can be tuned to give a particular moment of inertia of
a 1.4$M_{\odot}$ star ($I_{1.4}$) while keeping the other parameters fixed. We
can thus parameterise each EOS by $J$, $L$, $K_{\rm sym}$, $I_{1.4}$ and
$M_{\rm max}$.

In this work we want to concentrate on how the crust models and their
uncertainty affect the $i$-mode frequency, and so we fix the high density
degrees of freedom $I_{1.4}$ and $M_{\rm max}$. As we shall see, the $i$-mode
frequency is relatively insensitive to the stellar radius and therefore the
high density EOS. We choose to fix $M_{\rm max}$ at 2.2$M_{\odot}$,
comfortably above the maximum accurately measured pulsar mass (Cromartie et
al., 2020) and consistent with maximum masses inferred from modeling of the
binary neutron star merger resulting in GW170817. Given a value of $J$, $L$,
$K_{\rm sym}$, and $M_{\rm max}$ the moment of inertia $I_{1.4}$ can be
systematically varied between the minimum and maximum values allowed by
causality. As we demonstrate in Figure 18, this variation has only a small
effect on the $i$-mode frequencies calculated, and so for each value of $J$,
$L$, $K_{\rm sym}$, and $M_{\rm max}$ we choose the EOS model whose moment of
inertia $I_{1.4}$ is the average of the maximum and minimum possible values of
$I_{1.4}$ as the representative EOS.

Figure 2: The equations of state used in this work. The different colours are
for the three different sets of symmetry energy parameter ranges we have used,
with the black being for our Uniform ranges, blue for our PNM ranges, and red
for our MSL0 ranges. Outside the density range shown here, the equations of
state are indentical.

The resultant EOSs are shown in Figure 2. These are used with the Tolman-
Oppenheimer-Volkoff (TOV) equations (Oppenheimer & Volkoff, 1939; Tolman,
1939) to determine the grid of stellar models for which the $i$-mode
frequencies are calculated.

##  3 Calculation of the Normal Modes

We calculate the frequencies and radial/transverse displacements of the modes
of a neutron star by linearly perturbing the equations defining the
equilibrium state of neutron star matter. We assume that the binary lifetime
is much larger than the spin-down and cooling times of the individual neutron
stars, and so we may ignore rotation and high-temperature effects. We will
also assume that the frequencies of the modes that we are interested in ($\sim
10-100$ Hz) result in oscillations that are significantly faster than the
beta-equilibrium timescale. Therefore, weak-interactions do not have time to
change the composition of a displaced mass element to more closely match the
local composition.

###  3.1 Basic Equations

To construct a 2-component neutron star model (consisting of a solid crust and
fluid core) we follow McDermott et al. (1988), beginning with the mass
continuity equation, momentum conservation equation, and Poisson’s equation:

$\displaystyle\frac{\partial\rho}{\partial t}+\nabla\cdot(\rho v)=0,$ (15)
$\displaystyle\frac{\partial v}{\partial
t}=(v\cdot\nabla)v=\frac{1}{\rho}\nabla\cdot\sigma-\nabla\Phi,$ (16)
$\displaystyle\nabla^{2}\Phi=4\pi G\rho,$ (17)

where $\rho$ is the energy density, $v$ is the velocity of the matter,
$\sigma$ is the stress tensor, and $\Phi$ is the gravitational potential. By
combining the linear perturbations of these equations, and taking the Cowling
approximation by ignoring perturbations of the gravitational potential
(Cowling, 1941), we obtain the wave equation

$\displaystyle\omega^{2}u=-\nabla\left(\frac{\Gamma_{1}p}{\rho}\nabla\cdot
u\right)-\nabla\left(\frac{1}{\rho}u\cdot\nabla
p\right)-\hat{r}A\frac{\Gamma_{1}p}{\rho}\nabla\cdot u$
$\displaystyle+\frac{1}{\rho}\biggr{(}\nabla\left(\frac{2}{3}\mu\nabla\cdot
u\right)-\left(\nabla\mu\cdot\nabla\right)u-\nabla\left(u\cdot\nabla\mu\right)$
$\displaystyle+\left(u\cdot\nabla\right)\nabla\mu-\mu\left(\nabla^{2}u+\nabla\left(\nabla\cdot
u\right)\right)\biggr{)},$ (18)

where $u(x,t)$ is the Lagrangian displacement,

$\displaystyle
A=\frac{1}{\rho}\frac{d\rho}{dr}-\frac{1}{\Gamma_{1}P}\frac{dP}{dr}$ (19)

is the Schwarzschild discriminant, and $\Gamma_{1}$ is the adiabatic index
defined in equation (9). The non-diagonal terms of the stress tensor are given
by $\sigma_{ij}=\mu\nabla_{i}u_{j}$, assuming the isotropic shear modulus
$\mu$ from equation (8).

Taking the perturbations to have time dependence of the form $e^{i\omega t}$,
where $\omega$ is the mode frequency, we have

$\displaystyle u(x,t)=\xi(x)e^{i\omega t}.$ (20)

In spherical coordinates this can be further separated into radial and
transverse components:

$\displaystyle\xi_{r}=U(r)Y_{lm},\;\;\;\xi_{\theta}=V(r)\frac{\partial
Y_{lm}}{\partial\theta},\;\;\;\xi_{\phi}=\frac{V(r)}{\sin(\theta)}\frac{\partial
Y_{lm}}{\partial\phi},$ (21)

where $U(r)$ is the radial displacement, $V(r)$ is the transverse
displacement, and $Y_{lm}$ are the spherical harmonics.

By using the separation of variables given in equations (20) and (21), the
wave equation can be rewritten in terms of $U$ and $V$:111The version of
equation (22) in McDermott et al. (1988) has a typographical error.

$\displaystyle\rho\omega^{2}U=\rho\frac{d\hat{\chi}}{dr}-A\Gamma_{1}p\hat{\alpha}-\frac{d}{dr}\left(\frac{1}{3}\mu\hat{\alpha}\right)+\frac{d\mu}{dr}\left(\hat{\alpha}-2\frac{dU}{dr}\right)$
$\displaystyle-\mu\left(\frac{1}{r^{2}}\frac{d}{dr}\left(r^{2}\frac{dU}{dr}\right)-\frac{\ell(\ell+1)}{r^{2}}U+\frac{2\ell(\ell+1)}{r^{2}}V-\frac{2}{r^{2}}U\right),$
(22)
$\displaystyle\rho\omega^{2}V=\rho\frac{\hat{\chi}}{r}-\frac{1}{3}\frac{\mu\hat{\alpha}}{r}-\frac{d\mu}{dr}\left(\frac{dV}{dr}-\frac{V}{r}+\frac{U}{r}\right)$
$\displaystyle-\mu\left(\frac{1}{r^{2}}\frac{d}{dr}\left(r^{2}\frac{dV}{dr}\right)-\frac{\ell(\ell+1)}{r^{2}}V+\frac{2}{r^{2}}U\right),$
(23)

where:222The version of equation (25) in McDermott et al. (1988) has a
typographical error.

$\displaystyle\hat{\alpha}=\frac{1}{r^{2}}\frac{d}{dr}(r^{2}U)-\frac{\ell(\ell+1)}{r}V,$
(24)
$\displaystyle\hat{\chi}=-\frac{\Gamma_{1}p}{\rho}\hat{\alpha}-\frac{1}{\rho}\frac{\partial
p}{\partial r}U.$ (25)

In the crust these equations can be solved as a set of four first-order
differential equations, whereas in the core they are simplified by the
requirement that $\mu=0$, resulting in two first-order differential equations.

In this work we have applied Newtonian perturbations to a relativistic
equilibrium stellar structure, resulting in a hybrid model. In order for this
model to be usable, the modes must be orthogonal such that any perturbation
can be expressed as a unique linear combination of the modes. The
eigenfunction and eigenvalue of any mode can be defined in relation to the
oscillation operator, $\mathcal{H}$, as shown by Reisenegger (1994):

$\displaystyle\mathcal{H}\xi=-\omega^{2}\xi.$ (26)

Newtonian perturbations will only result in orthogonal modes if $\mathcal{H}$
is Hermitian with respect to the inner product of two vector fields (any two
displacements of matter within the star), ie:

$\displaystyle\int_{*}\rho_{0}(r)\zeta^{*}(x)\cdot\mathcal{H}\Psi(x)d^{3}x=\int_{*}\rho_{0}(r)\mathcal{H}\zeta^{*}(x)\cdot\Psi(x)d^{3}x.$
(27)

For a Newtonian stellar model the oscillation operator is Hermitian, and
therefore applying Newtonian perturbations results in orthogonal modes.

For a relativistic stellar model, the oscillation operator is not Hermitian.
However, this does not pose a problem for a relativistic perturbation approach
(see e.g. Yoshida & Lee, 2002) because their eigenfrequencies are complex
numbers, with the imaginary component arising from the damping of the mode due
to the emission of gravitational waves. This imaginary component cancels out
the deviation from orthogonality that arises from $\mathcal{H}$ not being
Hermitian, and so relativistic perturbations can be applied to a relativistic
stellar model to obtain orthogonal modes.

The hybrid model we have adopted can cause problems, because the stellar model
does not give a Hermitian oscillation operator and the eigenfrequencies of
Newtonian oscillations do not have the imaginary component required to cancel
out the modes’ deviation from orthogonality. To fix this, we follow
Reisenegger (1994) and define the local acceleration due to gravity within the
star as

$\displaystyle g=-\frac{1}{\rho}\frac{dP}{dr}.$ (28)

This form for the gravity makes the oscillation operator Hermitian within the
relativistic model, and therefore we can apply Newtonian perturbations to this
modified relativistic star to obtain orthogonal modes.

At the boundary between the solid crust and the fluid core, there are three
jump conditions that must be satisfied:333The version of the left-hand side of
equation (30) in McDermott et al. (1988) has a typographical error.

$\displaystyle U\rvert_{r=R_{\rm cc}^{+}}=U\rvert_{r=R_{\rm cc}^{-}},$ (29)
$\displaystyle\frac{1}{p}\left(\lambda\hat{\alpha}+2\mu\frac{dU}{dr}\right)\biggr{\rvert}_{r=R_{\rm
cc}^{+}}=\tilde{V}\left(\frac{U}{r}-\frac{\omega^{2}V}{g}\right)\biggr{\rvert}_{r=R_{\rm
cc}^{-}},$ (30)
$\displaystyle\frac{\mu}{p}\left(\frac{dV}{dr}-\frac{V}{r}+\frac{U}{r}\right)\biggr{\rvert}_{r=R_{\rm
cc}^{+}}=0,$ (31)

where

$\displaystyle\tilde{V}=-\frac{d\ln\left(P\right)}{d\ln\left(r\right)}=\frac{\rho
gr}{p},$ (32) $\displaystyle\lambda=\Gamma_{1}P-\frac{2}{3}\mu$ (33)

is the Lamé coefficient,

$\displaystyle M_{r}=4\pi\int^{r}_{0}r^{\prime 2}\rho(r^{\prime})dr^{\prime}$
(34)

is the mass contained within radius $r$, and $r=R_{\rm cc}^{+}$ ($r=R_{\rm
cc}^{-}$) indicates that the value is evaluated at the boundary when
approached from the crust (core) of the star. The three conditions require
different properties to be continuous across the crust-core boundary. The
first is for the radial displacement, the second is for the pressure, and the
third is for the transverse traction (which must be zero in the fluid core).
The first two jump conditions can be combined to cancel out the arbitrary
magnitude of $U$ and $V$ (which is different in the crust and the core),
giving us

$\displaystyle\frac{r}{p}\frac{\left(\lambda\hat{\alpha}+2\mu\frac{dU}{dr}\right)}{U}\biggr{\rvert}_{r=R_{\rm
cc}^{+}}=\tilde{V}\left(1-\frac{\sigma^{2}r}{g}\frac{V}{U}\right)\biggr{\rvert}_{r=R_{\rm
cc}^{-}}.$ (35)

This leaves us with two jump conditions, equations (31) and (35), and two
eigenvalues, $\Omega$ and $V(R_{*})$.

Every mode must satisfy the boundary conditions at the centre and surface of
the star. The conditions at the surface are based on the requirements that the
Lagrangian pressure perturbation and transverse traction go to zero at the
surface:

$\displaystyle\frac{1}{p}\left(\lambda\hat{\alpha}+2\mu\frac{dU}{dr}\right)\biggr{\rvert}_{r=R_{*}}=\tilde{V}\left(\frac{\tilde{V}-c_{1}\Omega^{2}-4+\tilde{U}}{\frac{\ell(\ell+1)}{c_{1}\Omega^{2}}-\tilde{V}}+1\right)\frac{U}{r}\biggr{\rvert}_{r=R_{*}},$
(36)
$\displaystyle\frac{\mu}{p}\left(\frac{dV}{dr}-\frac{V}{r}+\frac{U}{r}\right)\biggr{\rvert}_{r=R_{*}}=0,$
(37)

where every quantity is evaluated at the surface of the star ($r=R_{*}$), and

$\displaystyle\tilde{U}=\frac{d\ln\left(M_{r}\right)}{d\ln\left(r\right)}=4\pi
r^{2}\rho\frac{r}{M_{r}},\;\;\;c_{1}=\left(\frac{r}{R_{*}}\right)^{3}\frac{M_{*}}{M_{r}},\;\;\;\Omega^{2}=\frac{\omega^{2}R_{*}^{3}}{GM_{*}}$

are equilibrium properties of the star. The condition at the centre of the
star follows from the requirement that $U$ and $V$ be regular there:

$\displaystyle\left(\frac{c_{1}\Omega^{2}}{l}U-\frac{\sigma^{2}r}{g}V\right)\biggr{\rvert}_{r=0}=0,$
(38)

where every quantity is evaluated at the centre of the star ($r=0$).

###  3.2 The Crust-Core Interface Mode

We numerically solve for the eigenvalues by adjusting trial eigenvalues and
solving equations (22) and (23) until the jump conditions are satisfied,
indicating that a mode has been found. For $J=30$ MeV, $L=50$ MeV, and $K_{\rm
sym}=-80$ MeV the interface mode is found when the eigenvalues are $f=134.3$
Hz and $V(R_{*})=-7.72$, resulting in the radial and transverse displacements
shown in Figure 3. This mode has a distinctive peak in radial displacement at
the crust-core boundary, which is expected since the $i$-mode is caused by the
discontinuity between the crust and core. The transverse displacement in the
core is relatively small, with the discontinuity separating it from the larger
displacement in the crust. Thus, a larger fraction of the mode energy goes
into deforming the crust, helping it to reach the breaking strain faster. This
makes the crust-core $i$-mode a good candidate to power an RSF.

Figure 3: The quadrupole crust-core interface mode for the EOS parameterised
by $J=30$ MeV, $L=50$ MeV, and $K_{\rm sym}=-80$ MeV. $V(r)$ has been reduced
by an order of magnitude so that it can be plotted alongside $U(r)$.

##  4 Results

###  4.1 The Impact of Symmetry Energy Parameters on Stellar Structure

Figure 4 shows how the relationships between mass and stellar radius and
between mass and crust-core transition radius change with the symmetry energy
parameters $J$, $L$ and $K_{\rm sym}$. The maximum neutron star mass ($M_{\rm
max}$) and moment of inertia at $1.4M_{\odot}$ ($I_{1.4}$), which primarily
control the core EOS, have been fixed as described in Section 2.3. For many
different types of mode, the frequency is dependant on both the stellar radius
and the crust-core transition radius. Therefore, if this is the only thing we
consider, Figure 4 tells us that we would expect mode frequencies of $1.4$ M⊙
neutron stars to vary the most with $K_{\rm sym}$ and the least with $J$.
However, we must also consider the impact of the symmetry energy parameters on
the restoring forces that cause the modes to oscillate. For the $i$-mode, the
restoring forces are dominated by shear forces, and therefore we would expect
the impact of the symmetry energy parameters on the $i$-mode frequency to be
closely related to their impact on the shear speed,
$c_{t}=\sqrt{\frac{\mu}{\rho}}$. Figure 5 shows how the relationship between
the stellar mass and the density-weighted average of the shear speed in the
crust changes with the symmetry energy parameters. This figure shows that $J$
and $L$ have larger impacts on the average shear speed than $K_{\rm sym}$, and
that the shear speed is strongly dependent on all three of the symmetry energy
parameters.

Parameter | Change in | Change in | Change in
---|---|---|---
varied |  $R_{*}$ (Km) |  $R_{\rm cc}$ (Km) |  $\bar{c}_{t}$ (cm/s)
|  |  |
$J$ ($25\rightarrow 35$ MeV) | $0.21$ | $-0.23$ | $3.8\times 10^{6}$
|  |  |
$L$ ($20\rightarrow 70$ MeV) | $0.47$ | $0.48$ | $-4.9\times 10^{6}$
|  |  |
$K_{\rm sym}$ ($-200\rightarrow 40$ MeV) | $0.87$ | $0.79$ | $1.8\times 10^{6}$
|  |  |
$I_{1.4}$(Imin $\rightarrow$ Imax) | $1.40$ | $1.12$ | $4.1\times 10^{3}$
Table 1: Typical changes in stellar radius, crust-core transition radius, and
density-weighted shear speed of a $1.4$ M⊙ NS caused by varying each of the
symmetry energy parameters. For each parameter, we vary it over the specified
range while holding the others constant in the middle of their Uniform ranges.
For comparison, the EOS in the middle of all of the Uniform ranges ($J=30$,
$L=50$ and $K_{\rm sym}=-80$ MeV) results in $R_{*}=11.92$ Km, $R_{cc}=10.95$
Km, and $\bar{c}_{t}=4.8\times 10^{6}$ cm/s. Also shown are the changes due to
varying the moment of inertia of a $1.4$ M⊙ NS (which determines the core EOS)
across the full range of values allowed by causality while holding the three
symmetry energy parameters constant and keeping the maximum mass fixed to
M${}_{\rm max}=2.2$ M⊙.

In Table 1 we quantify the typical impact that varying the symmetry energy
parameters has on the properties of a $1.4$ M⊙ neutron star. From this table,
and the trends of Figures 4 and 5, we see that varying the symmetry energy
parameters causes a fractional change in the average shear speed that is
significantly greater than the fractional change in the stellar radius or the
transition radius. Therefore, we expect that the symmetry energy parameters’
relationships to the $i$-mode frequency will be dominated by their relative
contributions to the average shear speed. We will investigate this further in
Section 4.3, after we have calculated the dependence of the frequency on the
symmetry energy parameters.

So far we have ignored the uncertainty in the parameters which control the
core EOS. To address this, in Table 1 we also give the changes in stellar
properties caused by varying the $1.4$ M⊙ moment of inertia over an extremely
conservative range. We find that this causes $\sim 10$% changes in the stellar
radius and crust-core transition radius. These changes, while significant, are
much smaller than the order unity changes in shear speed caused by varying the
symmetry energy parameters. We therefore expect that, when compared to the
symmetry energy parameters, the moment of inertia (and thus core EOS) has
little impact on the $i$-mode frequency, and so we shall keep it fixed. The
validity of this choice will be discussed further in Section 5. The maximum
mass is more sensitive to the EOS at higher densities than the moment of
inertia is, and so our results will be less sensitive to the choice of the
maximum mass.

Figure 4: The relationships between the neutron star mass and the stellar
radius (solid lines) and the crust-core transition radius (dashed lines) for
different EOSs. Each plot varies a different symmetry energy parameter over a
wide range of values, with the lines being labelled with the varied parameter.
The red (middle) lines of each plot are the same.

Figure 5: The relationship between the neutron star mass and the density-
weighted average of the shear speed for different EOSs. Each plot varies a
different symmetry energy parameter over a wide range of values, with the
lines being labelled with the varied parameter. The red (middle) lines of each
plot are the same.

###  4.2 Interface Mode dependence on Nuclear Parameters

For the remainder of this paper, unless otherwise stated, we focus our results
on $1.4M_{\odot}$ neutron stars. We explored the three different ranges of
symmetry energy parameters described in Section 2.1. In Section 4.2.1, we use
our uniform (weakly-constrained) $J$, $L$ and $K_{\rm sym}$ ranges, in order
to avoid tying our results to those of previous works. In Section 4.2.2 we use
our PNM parameter ranges, where the parameters are consistent with the results
of pure neutron matter calculations as this is the most relevant constraint
for neutron star matter, which is extremely neutron rich. Finally, in Section
4.2.3 we use our MSL0-like parameter ranges, where $K_{\rm sym}$ is defined as
a particular function of $J$ and $L$. This lets us more directly compare with
previous works that have only allowed the first two symmetry energy parameters
to vary, such as Chen et al. (2010), Steiner & Gandolfi (2012) and Tsang et
al. (2009).

####  4.2.1 Uniform (weakly-constrained) $J$, $L$ and $K_{\rm sym}$ ranges

We constructed a set of EOSs which had symmetry energy parameters evenly
spaced in the three-dimensional parameter space defined for our uniform
distribution in Section 2.1 (we used a $J$ spacing of $1$ MeV, $L$ spacing of
$10$ MeV, and $K_{\rm sym}$ spacing of $40$ MeV). After using the TOV
equations to obtain a stellar model for each EOS, we calculated their $\ell=2$
$i$-mode frequencies. We then interpolated between these frequencies to find
surfaces of constant frequency in the $J$,$L$,$K_{\rm sym}$ parameter space,
shown in Figure 6.

Figure 6: Surfaces of constant $i$-mode frequency in the weakly-constrained
symmetry energy parameter space of our uniform distribution. The surfaces are
$40$ Hz apart. The data to manipulate this plot can be found at
https://github.com/davtsang/RSFSymmetry/.

In order to better understand Figure 6 we plot its two-dimensional projection
on the $J-L$ plane, shown in Figure 7. This figure shows the values of $J$ and
$L$ that can result in the $i$-mode having the chosen frequency, with the
spread in $L$ at any given $J$ being due to the range of possible $K_{\rm
sym}$ values. We could also plot the projections on the $J-K_{\rm sym}$ and
$L-K_{\rm sym}$ planes. However, we find that the strong dependence of the
frequency on $L$ and $J$ means that these plots are uninformative, since the
variation in the projected parameter can cause the $i$-mode frequency regions
to cover almost the entire parameter space.

Figure 7: The two-dimensional projection of Figure 6 on the $J-L$ plane,
showing the ranges of the symmetry energy parameters in which five example
$i$-mode frequencies can be obtained. The widths of the regions are caused by
$K_{\rm sym}$ having any value between $-200$ MeV and $40$ MeV.

One of the uncertainties affecting the constraints we could put on the
symmetry energy parameters is the timescale over which resonant excitation of
modes can occur. This can be calculated as (Tsang et al., 2012b)

$\displaystyle t_{\rm res}\sim 8\times
10^{-2}\text{s}\left(\frac{\mathcal{M}}{1.2\text{M}_{\odot}}\right)^{\frac{-5}{6}}\left(\frac{f_{\rm
mode}}{100\text{ Hz}}\right)^{\frac{-11}{6}},$ (39)

where

$\displaystyle\mathcal{M}=\frac{M_{1}^{\frac{3}{5}}M_{2}^{\frac{3}{5}}}{(M_{1}+M_{2})^{\frac{1}{5}}}$
(40)

is the chirp mass (Cutler & Flanagan, 1994). This timescale can be combined
with the rate of change of the gravitational wave frequency

$\displaystyle\dot{f}_{\rm gw}=\frac{f_{\rm gw}}{4.7\times
10^{-3}\text{s}}\left(\frac{\mathcal{M}}{1.2\text{M}_{\odot}}\right)^{\frac{5}{3}}\left(\frac{f_{\rm
gw}}{1000\text{ Hz}}\right)^{\frac{8}{3}}$ (41)

(where $f_{\rm gw}\approx f_{\rm mode}$) to obtain a simple estimate of the
range of frequencies over which resonance can occur:

$\displaystyle\delta f\sim t_{\rm res}\dot{f}_{\rm gw}\sim
3.7\text{Hz}\left(\frac{\mathcal{M}}{1.2\text{M}_{\odot}}\right)^{\frac{5}{6}}\left(\frac{f_{\rm
mode}}{100\text{ Hz}}\right)^{\frac{11}{6}}.$ (42)

From this we see that the width of the resonance window increases with the
frequency, with $\delta f$ scaling as $f_{\rm mode}^{\frac{11}{6}}$. For a
chirp mass of $1.2M_{\odot}$ and a resonance at $100$ Hz, we get a frequency
range of $\delta f\sim 3.7$ Hz, and for a resonance at $160$ Hz we get a range
of $\delta f\sim 8.8$ Hz. This means that the spread of the frequency regions
in the $J-L$ plane (seen in Figure 8) is quite small ($\delta L\lesssim 5$
MeV), and therefore the impact of the resonance window is significantly less
than that of the $K_{\rm sym}$ range, which causes a spread of $\sim 20$ MeV
in $L$. It should be noted that this is a very conservative estimate, and that
by more accurately calculating both the rate at which energy is transferred
into the modes and the breaking strain of the crust, this uncertainty could be
significantly reduced by calculating the time it takes for the crust to
shatter.

Figure 8: The spread in the constraints on $J$ and $L$ caused by the resonance
window (ie: the $i$-mode frequency can have any value between $f_{GW}-\delta
f$ and $f_{GW}+\delta f$). We have fixed $K_{\rm sym}=-80$ MeV to isolate the
effect of the resonance window.

In Figure 9 we plot regions in the $J-L$ plane that result in $i$-modes that
can be resonantly excited by certain chosen GW frequencies when considering
both the $K_{\rm sym}$ range (shown in Figure 7) and the resonance window
(given by equation (42)). These regions are compared to the combined
experimental nuclear constraints given in Lattimer & Lim (2013) (which
includes constrains from: fits to nuclear masses (Kortelainen et al., 2010),
neutron skin thickness (Chen et al., 2010), dipole polarisability (Piekarewicz
et al., 2012), giant dipole resonances (Trippa et al., 2008), and isotope
diffusion in heavy ion collisions (Tsang et al., 2009)). From these results,
in order to be consistent with the combined experimental nuclear constraints,
we could expect to observe precursor flares in the range $120\lesssim
f_{GW}\lesssim 280$ Hz. This range is very wide as it is based on our most
conservative constraints on $J$ and $L$, with the upper bounds of Figure 9
using $f_{i}=f_{GW}-\delta f$ and $K_{\rm sym}=40$ MeV (or the maximum $K_{\rm
sym}$ with a stable crust), and the lower bounds using $f_{i}=f_{GW}+\delta f$
and $K_{\rm sym}=-200$ MeV (or the minimum $K_{\rm sym}$ with a stable crust).
In order to find more useful constraints on the symmetry energy parameters, we
can reduce the $J$,$L$,$K{\rm sym}$ parameter space used to generate the EOSs.

Figure 9: Examples of constraints that could be applied to $J$ and $L$, with
each constraint being labelled with the GW frequency at the time of the RSF.
The width of the constraints comes from allowing $K_{\rm sym}$ to have any
value in the range $-200$ to $40$ MeV at all $J$ and $L$ values, and from
approximating the uncertainty in the $i$-mode frequency due to the resonance
window with equation (42). The hatched area in the centre of the plot
represents the combined experimental nuclear constraints on $J$ and $L$
(Lattimer & Lim, 2013).

####  4.2.2 $J$, $L$ and $K_{\rm sym}$ constrained using pure neutron matter
theory (PNM)

There are many different constraints on the nuclear symmetry energy parameters
that we could use to reduce the $J$,$L$,$K_{\rm sym}$ parameter space. To keep
our results conservative, we consider only the most relevant constraints so as
to not make our results overly dependent on other works. For neutron star
matter, which is extremely neutron rich, one such constraint comes from
calculations of the properties of pure neutron matter (see Section 2.1 and the
PNM ranges discussed in Newton & Crocombe (2020)). With this additional
constraint on the symmetry energy parameter ranges, we repeat our method from
Section 4.2.1 for obtaining surfaces of constant frequency in the
$J$,$L$,$K_{\rm sym}$ parameter space, resulting in Figure 10.

Figure 10: Surfaces of constant $i$-mode frequency in the $J$,$L$,$K_{\rm
sym}$ parameter space, where the $J$, $L$ and $K_{\rm sym}$ ranges are
constrained by pure neutron matter theory. The surfaces are 20 Hz apart. The
viewing angle is the same as for Figure 6, which causes the surfaces to
collapse into lines. The data to manipulate this plot can be found at
https://github.com/davtsang/RSFSymmetry/.

Figure 11 shows two-dimensional projections of Figure 10. Similar to Figure 7,
the first plot of this figure is the projection on the $J-L$ plane. However,
by constraining the parameter space with the results of pure neutron matter
theory we have reduced the range of $K_{\rm sym}$ values, and therefore the
widths of the frequency regions are much smaller. Similarly, the ranges of $J$
and $L$ have been reduced, and therefore the projections of Figure 11 on the
$J-K_{\rm sym}$ and $L-K_{\rm sym}$ planes are now informative. These are
shown in the second and third plots of Figure 11, where the widths of the
frequency regions are determined by the ranges of $L$ and $J$ (respectively).
In these three plots, the widths of the frequency regions show the impact of
the projected symmetry energy parameter; the wider the regions, the more
significant the uncertainty in the projected parameter is to the mode
frequency.

Figure 11: The two-dimensional projections of Figure 10, showing the ranges of
the symmetry energy parameters in which five example $i$-mode frequencies can
be obtained while being consistent with pure neutron matter theory.

To calculate the constraints that we could place on the symmetry energy
parameters, we combine the frequency regions shown in Figure 11 with the
uncertainty in the $i$-mode frequency due to the resonance window, given by
equation (42). This results in Figure 12, which shows example constraints that
could be applied to the symmetry energy parameters in the event of RSF
detections at certain GW frequencies, alongside the combined experimental
nuclear physics constraints on $J$ and $L$ from Lattimer & Lim (2013). By
comparing Figure 9 and the first plot of Figure 12 we can see that restricting
the symmetry energy parameters to the ranges predicted by pure neutron matter
theory has significantly tightened our constraints, making them competitive
with the experimental constraints. By inverting our method, we find that
$120\lesssim f_{GW}\lesssim 180$ Hz results in constraints on $J$ and $L$ that
are consistent with the combined experimental nuclear constraints.

Figure 12: The constraints that we could apply to the symmetry energy
parameters in the event of an RSF detection at different GW frequencies, where
the $J$, $L$ and $K_{\rm sym}$ ranges are constrained by pure neutron matter
theory. The first (second) [third] plot shows the constraints on $J$ and $L$
($J$ and $K_{\rm sym}$) [$L$ and $K_{\rm sym}$], with the strength of each
constraint being determined by the $K_{\rm sym}$ ($L$) [$J$] range and the
width of the resonance window. The hatched area in the centre of the first
plot indicates the combined experimental nuclear constraints on $J$ and $L$.

####  4.2.3 Ksym as a function of J and L (MSL0)

To compare our work to more restricted two-parameter Skyrme models, we
reproduce the MSL0 model’s $K_{\rm sym}$ dependence on $J$ and $L$. Note that
this dependence does not have any special physical significance, and there are
other relationships between the symmetry energy parameters that are equally
plausible. Using a similar grid of $J$ and $L$ values as in Section 4.2.1, we
calculated the $i$-mode frequencies for a set of stellar models to obtain
Figure 13. This figure also shows the approximate relationship between
frequency and time before coalescence, given by (Tsang et al., 2012b;
Blanchet, 2006)

$\displaystyle t_{c}-t=\frac{3}{8}t_{GW}=1.76\times
10^{-3}\rm{s}\left(\frac{\mathcal{M}}{1.2M_{\odot}}\right)^{-\frac{5}{3}}\left(\frac{f_{GW}}{1000\rm{Hz}}\right)^{-\frac{8}{3}}.$
(43)

Figure 13: The relationship between $J$, $L$ and the $i$-mode frequency for
EOSs with the MSL0-like $K_{\rm sym}$ dependence. As there is no $K_{\rm sym}$
range, this plot shows every EOS used in Section 4.2.3. The right axis of the
plot shows the approximate times before coalescence at which the frequencies
occur.

We interpolated between the grid of $J$ and $L$ values to obtain frequency
contours in the $J-L$ plane. These contours are spread by the resonance window
calculated with equation (42), resulting in the constraints shown in Figure 14
(where we have also plotted the nuclear physics constraints on $J$ and $L$).
These results represent a best case scenario for the $K_{\rm sym}$ range, as
there is no uncertainty in its value at all $J$ and $L$ values. From this
figure we can see that a precursor flare detected when $130\lesssim
f_{GW}\lesssim 170$ Hz would provide constraints on $J$ and $L$ consistent
with those from nuclear physics.

Figure 14: Frequency contours in the $J-L$ plane that have been spread by the resonance window (with the MSL0-like $K_{\rm sym}$ dependence). The hatched area in the centre of the plot indicates the combined experimental nuclear constraints on $J$ and $L$. $J$, $L$, $K_{\rm sym}$ | GW frequency | Time before
---|---|---
ranges | of the RSF (Hz) | coalescence (s)
|  |
Uniform |  $120$-$280$ |  $0.50$-$0.05$
|  |
PNM |  $120$-$180$ |  $0.50$-$0.17$
|  |
MSL0 |  $130$-$170$ |  $0.41$-$0.20$
Table 2: Summary of, for each of our data sets, the approximate range of
gravitational wave frequencies (and corresponding times before coalescence,
assuming a chirp mass of $1.2$ $M_{\odot}$) in which a Resonant Shattering
Flare needs to occur in order to be consistent with the combined experimental
nuclear constraints on $J$ and $L$ given in Lattimer & Lim (2013). These
results are for the Skyrme model described in Newton & Crocombe (2020), using
Newtonian perturbations and a $M_{*}=1.4$ M⊙ neutron star. They do not include
the impact of the core parameterisation, which is fixed as described in
Section 2.3.

Table 2 inverts our method for constraining the symmetry energy parameters by
showing the approximate range of gravitational wave frequency in which we
would expect to observe an RSF in order for our constraints on $J$ and $L$ to
be consistent with the combined experimental constraints (Lattimer & Lim,
2013) (ie: they have a non-zero overlap). From this we can see that, for the
model used in this work and a $1.4$ M⊙ neutron star, we expect to observe RSFs
at gravitational wave frequencies of around $150$ Hz, or approximately $0.3$ s
before coalescence. This is similar to the time before the main SGRB that many
precursors are observed ($~{}0.1-5.0$ s), providing evidence that these
precursors are RSFs.

###  4.3 Shear speed

In order to determine the cause of the change in $i$-mode frequency due to
variations in $J$, $L$ and $K_{\rm sym}$, we investigated how the properties
of the star we discussed in Section 4.1 relate to the frequency. As we
predicted, the frequency was strongly dependent on the density-weighted
average of the shear speed in the crust. This is shown in Figure 15, which
relates the frequency and the average shear speed for our three sets of EOSs.

Figure 15: The relationship between the $i$-mode frequency and the density-
weighted average shear speed for each of our EOSs. The black markers are for
EOSs in our Uniform set, the blue for EOSs in our PNM set, and the red for
EOSs in our MSL0 set.

In a similar way to how the first plot of Figure 11 shows the $J$ and $L$
values that can result in the $i$-mode having chosen frequencies, Figure 16
shows the $J$ and $L$ values that can result in stars with chosen average
shear speeds. The similarities between the regions shown in this figure and in
Figure 11 indicate that the $i$-mode frequency and average shear speed are
closely linked. At higher $J$ and $L$ values these figures become less
similar, suggesting that the significance of other stellar properties
increases with $J$ and $L$. In this figure we have only shown the results for
our PNM set of EOSs, since all three sets of EOSs give the same qualitative
results.

Figure 16: Similar to the first plot of Figure 11, but for the density-
weighted average shear speed. The spread is due to the $K_{\rm sym}$ range,
which is consistent with pure neutron matter theory, and the regions are
labelled in $10^{6}\rm{cm/s}$.

##  5 Discussion

For all three sets of EOSs, our constraints on $J$ and $L$ (shown in Figures
9, 12 and 14) are angled in the same direction as the combined constraints
from other works (Lattimer & Lim, 2013). Therefore, a detection of an RSF at a
frequency in the middle of the range that is consistent with these constraints
would provide a small improvement to our knowledge of $J$ and $L$. However, if
an RSF were to be detected at a higher or lower frequency our constraints
could be more interesting due to their overlap with the combined experimental
nuclear constraints being smaller.

The shear speed increasing with $J$ and $K_{\rm sym}$ and decreasing with $L$
is correlated with the impact of these changes on the symmetry energy in the
crust (where $n_{\rm b}<n_{0}$); increasing $J$ and $K_{\rm sym}$ and
decreasing $L$ causes the symmetry energy at crustal densities to increase, as
can be seen in Figure 17. Here the markers indicate the crust-core transition,
although note that the crust-core transition density does not strongly
correlate with the crust thickness - see e.g. (Ducoin et al., 2011).
Increasing the symmetry energy increases the energy cost of creating more
neutrons, and therefore decreases the fraction of dripped neutrons in the
crust. Equation (7) shows that as the dripped neutron fraction decreases, the
ion number density increases (for a fixed mass number of nucleus in the
crust). This leads to an increase in the shear modulus, as can be seen in
equation (6) or (8). By calculating the density-averaged neutron fraction
$\bar{X}_{\rm n}$, mass number $\bar{A}$ and charge number $\bar{Z}$ we have
confirmed that changes in $\bar{X}_{\rm n}$ are the dominant outcome of
varying any of the three symmetry energy parameters.

Together, Figures 17 and 15 provide a qualitative physical understanding of
the $i$-mode frequency dependence shown in Figures 9, 12, and 14. The symmetry
energy profile of the star determines the composition and shear modulus in the
inner crust, on which the $i$-mode frequency is strongly dependent.

Figure 17: The symmetry energy below nuclear saturation density for different
combinations of symmetry energy parameters. Each plot varies a different
parameter over a wide range of values, with the lines being labelled with the
varied parameter. The markers indicate the crust-core transition density,
showing how it varies with the symmetry energy parameters. The red (middle)
lines of each plot are the same.

It is important to note that equation (6) is a fit to calculations of the
shear modulus for an ionic lattice with no dripped neutrons and ionic
separations typical of the outer crust. The shear modulus of the deep inner
crust [including the nuclear pasta layers (Pethick & Potekhin, 1998)] remains
an important outstanding problem, the result of which might significantly
affect our results. However, the fact that the shear modulus depends on the
ion separation, which in turn depends on the fraction of dripped neutrons,
means that the relationship between the $i$-mode frequency and the symmetry
energy is likely to persist.

To test that our results are not significantly affected by the choice of core
EOS, we can investigate the impact of allowing the moment of inertia parameter
(I1.4) to vary between the minimum and maximum values allowed by causality.
Having a range of I1.4 values does not noticeably affect the shear speed, but
it does increase the range of stellar radii obtained with our sets of EOSs. In
Table 1 we show the typical changes in relevant stellar properties caused by
increasing I1.4 from its minimum to maximum value while holding the symmetry
energy parameters constant. If the symmetry energy parameters are all allowed
to vary in their Uniform ranges, and I1.4 is held constant at the average of
the minimum and maximum values allowed by causality, the radius of a $1.4$ M⊙
NS ranges from $10.6$ to $12.7$ Km. If we also allow I1.4 to vary between its
minimum and maximum values, the radius ranges from $9.90$ to $13.6$ Km.Both of
these radius ranges are obtained while assuming that the NS maximum mass is
$2.2$ M⊙. From Figure 18 we can see that the impact of the core EOS on the
$i$-mode frequency is negligible at low $L$ and $J$, and is still small at
higher values. This is because the change in radius caused by the core EOS is
far less significant than the change in shear speed caused by the symmetry
energy parameters. This illustrates that resonant shattering flares mainly
probe the EOS and composition of the neutron star crust, in contrast to tidal
deformability measurements that give information about the core EOS.

We can quantify the impact of the core EOS by calculating the change in
$i$-mode frequency caused by varying I1.4. For all $J$, $L$ and $K_{\rm sym}$
values in our ‘uniform’ ranges, compared to the average value of I1.4 we find
that the maximum and minimum I1.4 cause approximately $-5\%$ and $+6\%$
changes in $i$-mode frequency (respectively). As the I1.4 range used here is
only constrained by causality, it is extremely conservative and therefore the
choice of core EOS does not significantly affect our results. We also note
that the same event that results in a coincident detection of an RSF can be
used to extract the tidal deformability. This parameter constrains the core
EOS in a complimentary manner to the constraints explored in this work, with
the added restriction that neutron star masses and EOSs are the same in the
description of each phenomenon.

Figure 18: The effect of changing I1.4 (which is used to parameterise the core
EOS) on the symmetry energy parameters that give a frequency of $150$ Hz,
where we have fixed $K_{\rm sym}=-80$ MeV to isolate the impact of the core
EOS. Imin is the minimum moment of inertia allowed by causality, Imax is the
maximum, and Imid is their average. For the set of EOSs plotted here, Imin
ranges from 1.03$\times 10^{45}$ g cm2 to 1.2$\times 10^{45}$ g cm2 and Imax
ranges from 1.44$\times 10^{45}$ g cm2 to 1.57$\times 10^{45}$ g cm2. The
moment of inertia increases with $J$ and $L$.

We have assumed a neutron star mass of $1.4$ M⊙, but from Figure 19 we can see
that a realistic degree of uncertainty in the stellar mass (Abbott et al.,
2017) has a noticeable impact on the symmetry energy parameters that give a
chosen $i$-mode frequency. The change in the chosen frequency contour is
similar to the impact of the core EOS shown in Figure 18. However, the moment
of inertia range used in Figure 18 is very conservative, and so the
uncertainty in the neutron star mass measurement is likely to have a more
significant impact on the symmetry energy parameter constraints than the
uncertainty in the core EOS. Uncertainty in the mass of an RSF’s source should
be considered when calculating constraints on the symmetry energy parameters,
as its impact is similar to that of the resonance window (as can be seen by
comparing Figures 19 and 8).

Figure 19: The effect of changing the total NS mass on the symmetry energy
parameters that give a frequency of $150$ Hz, where we have fixed $K_{\rm
sym}=-80$ MeV to isolate the impact of the stellar mass. These mass values are
the 90% confidence ranges calculated by Abbott et al. (2017) for the stars in
GW170817 (using low-spin priors).

In this work we have used the hybrid approach of non-relativistic
perturbations of a relativistic star to obtain the wave-equation (ignoring
dynamic perturbations of the gravitational potential). Using relativistic
perturbations (while still ignoring metric perturbations) as in Yoshida & Lee
(2002) can result in $\sim 10$% changes in the mode frequencies. This is
significant when compared to the width of our constraints on $J$ and $L$, and
we will explore this effect on the constraints in a future work.

Accurately calculating the Schwarzschild discriminant is not simple, and as it
has little impact on the $i$-mode we have set it to zero (i.e. we have assumed
the star is barotropic). We have also assumed that the binary lifetime is much
longer than the neutron stars’ spin-down and cooling timescales, and so we
have ignored rotational and high temperature effects. Finally, we have not
considered the impact of superfluidity in the core of a neutron star.
Superfluidity allows protons and neutrons to move somewhat independently of
each-other, introducing a new set of counter-moving normal modes (Andersson &
Comer, 2001), as well as modifying the frequency of modes that mainly
oscillate within the core of the star. In the inner crust, partial entrainment
of the superfluid may change the shear speed by reducing the effective density
that accelerates due to shear forces.

##  6 Conclusion

We have calculated the relationship between the neutron star interface mode
frequency and the first three parameters that characterise the density
dependence of the nuclear symmetry energy at saturation density ($J$, $L$ and
$K_{\rm sym}$). This was done by using an extended Skyrme mean-field model for
the crust and outer core of the star, supplemented by two polytropes that
controlled the high density EOS. We have used this to present potential
constraints on the symmetry energy parameters that could be obtained by
coincident multimessenger detection of a Resonant Shattering Flare and
gravitational wave chirp during a binary neutron star inspiral. These
constraints have been shown to be competitive with current nuclear
experimental constraints.

Previous works have shown (Abbott et al., 2017; Bauswein et al., 2017; Abbott
et al., 2018; De et al., 2018; Abbott et al., 2019) that the gravitational
wave chirp from a binary neutron star inspiral (with sufficient signal-to-
noise ratio) can constrain the tidal deformabilities, masses, and radii of the
stars. These, in turn, place constraints on the neutron star equation of state
(Read et al., 2009b; Read et al., 2013; Lackey & Wade, 2015; Margalit &
Metzger, 2017; Annala et al., 2018; Lim & Holt, 2018; Most et al., 2018;
Fattoyev et al., 2018; Carson et al., 2019; Zhang & Li, 2019; Landry & Essick,
2019; Essick et al., 2020), primarily in the core. In this work, we have
examined the nuclear physics constraints (in particular on the nuclear
symmetry energy parameters $J$, $L$, and $K_{\rm sym}$) that could be obtained
by a future detection of a Resonant Shattering Flare along with a
gravitational-wave chirp. Timing of the RSF relative to the GW chirp can
provide a direct measurement of the resonant frequency of the $\ell=2$ core-
crust interface mode (Tsang et al., 2012b). This frequency is dependent on
properties of the neutron star near the core-crust boundary, and is thus
sensitive to the nuclear symmetry energy parameters which determine (in a
model dependent way), the properties of the neutron star near nuclear
saturation. The measurement of an $i$-mode frequency through coincident timing
of an RSF would provide astrophyiscal constraints orthogonal to those
sensitive mainly to the core EOS.

Following Newton & Crocombe (2020), we constructed three sets of EOSs
parameterised by $J$, $L$ and $K_{\rm sym}$, with each set allowing these
parameters to have different ranges. The high density EOS parameters were
fixed by choosing a reasonable value for $M_{\rm max}$ and a representative
value of the moment of inertia of a 1.4Modot star, $I_{1.4}$. Solving for the
$i$-mode frequencies, we were able to determine the region in the parameter
space to which $J$, $L$ and $K_{\rm sym}$ could be constrained given
measurements of different frequency values.

Multimessenger coincident timing of an RSF would give the $i$-mode frequency
to a precision roughly determined by the duration of the flare. Additionally,
taking the conservative assumption that the nuclear symmetry energy parameters
are consistent with the results of pure neutron matter theory provides
constraints on $J$, $L$, and $K_{\rm sym}$ that are competitive with
Kortelainen et al. (2010), Chen et al. (2010) and Tsang et al. (2009).
Conversely, we can use the constraints found by other works to obtain the
range of frequencies in which we would expect to observe an RSF for a $1.4$ M⊙
neutron star. For the models used in this work, the range predicted by pure
neutron matter theory is $\sim 120-180$ Hz.

We have shown that it is important to take into account the variation of the
third symmetry energy parameter, $K_{\rm sym}$, independent of $J$ and $L$.
For example, if we allow all three to vary, our predicted range of $i$-mode
frequencies is $120$-$280$ Hz, while if $K_{\rm sym}$ is restricted by a
choice of model, an artificially smaller range is predicted ($130$-$170$ MeV
in the case of the MSL0 model considered here). Conversely, experimental
measurements of $K_{\rm sym}$ will constrain the predicted range of
frequencies.

In Figure 15 we showed that the $i$-mode frequency (a global property of the
NS) is strongly dependent on the average (density-weighted) shear speed within
the crust (a local material property of the crust). Therefore, the dependence
of the frequency on the symmetry energy parameters is dominated by their
effects on the shear modulus within the crust, and in particular near the
crust-core boundary. Figures 5 and 16 related the shear speed to the symmetry
energy parameters (similar to Figure 1 of Steiner & Watts, 2009), connecting
changes in these nuclear physics parameters to their impact on the average
shear speed in the crust. While other global properties of the stellar
structure (e.g. neutron star radius, radius of the core-crust transition)
which vary with the model parameters ($J$, $L$, $K_{\rm sym}$, $M_{\rm max}$
and $I_{1.4}$) also play a role, we found that the $i$-mode frequency depends
most strongly on the average shear speed, as can be seen from the similarities
between Figure 16 and Figure 12.

The quantitative results presented in this work are model dependant. Our focus
is on constraining the symmetry energy parameters that characterize the crust
and out core EOS, so that is where we span the widest range of the available
parameter space by independently varying the first three parameters in the
Taylor expansion around saturation density. However, we do restrict the
parameter space of $J$ and $L$ to that spanned by nuclear experimental
constraints; notably, this means values of $L$ are, for the most part, below
90 MeV. This excludes some of the stiffest EOSs, and therefore the neutron
star models with the largest possible radii. A recent measurement of the
neutron skin of 208Pb suggests that the slope of the symmetry energy $L$ may
be significantly above 100 MeV (Reed et al., 2021), which, although at odds
with most other experimental results, is a reminder we should not rule out
stiffer EOSs. The parameter space of the high density EOS, consisting of two
polytropes, is restricted to a maximum neutron star mass of $2.2$ M⊙, and a
moment of inertia of a 1.4$M_{\odot}$ star in the middle of the range allowed
by causality. Using different EOS models may result in significantly different
frequencies, with Tsang et al. (2012b) showing $i$-mode frequencies as low as
$30$ Hz. However, we have investigated the impact of variation in the moment
of inertia parameter and found that it had little impact on the $i$-mode, as
it did not significantly affect the shear speed. Therefore, when choosing
models for use in the analysis of RSFs, the description of the neutron star
crust is the most important input. A exploration of the wider parameter space
including high-$L$ EOSs will be the subject of future work.

A number of upcoming nuclear experiments promise to constrain the symmetry
energy further. We highlight the ongoing efforts to extract the neutron skin
of neutron rich nuclei from measurements of the parity-violating asymmetry in
the electron scattering cross-section caused by the weak interaction
(Abrahamyan et al., 2012) at Jefferson Lab and Mainz Superconducting
Accelerator (Horowitz et al., 2014; Becker et al., 2018; Thiel et al., 2019).
The latter is responsible for the recent measurement of the neutron skin
mentioned above. As illustrated in figure 1 of Steiner & Watts (2009), neutron
skins provide a constraint on the symmetry energy that is orthogonal to those
provided by the constraints on the shear speed and hence the $i$-mode
frequency. Powerful constraints may be obtained in the future by combining
these weak, EM and gravitational-wave observations to probe the strong force
in multi-messenger nuclear astrophysics.

Using upcoming LIGO/Virgo/KAGRA observing runs (Abbott et al., 2020), and
existing Gamma-ray burst monitors such as Swift/BAT (Barthelmy et al., 2005)
and Fermi/GBM (Meegan et al., 2009) to provide coincident timing, the
detection of a Resonant Shattering Flare during a binary neutron star inspiral
can provide a new complementary astrophysical constraint on nuclear physics
parameters by probing the bulk properties of neutron star matter near the
crust/core transition. The rates of RSFs are currently uncertain, with
precursor flares estimated to occur for $\sim 3-10$% of SGRBs. However, the
recent coincident detection of an (off-axis) SGRB and the chirp from GW170817
suggests a rate of NS mergers such that we may soon be able to obtain these
powerful constraints.

## Acknowledgements

DN is supported by a University Research Studentship Allowance from the
University of Bath, and WGN acknowledges support from NASA grant
80NSSC18K1019. We would like to thank the anonymous referee whose insightful
suggestions have greatly improved the clarity and presentation of this work.

## Data Availability

The code to calculate the stellar models and $i$-mode frequencies, along with
the tabulated EOSs and compositions used for the grids of symmetry energy
parameters are provided via https://github.com/davtsang/RSFSymmetry/.

## References

  * Abbott et al. (2017) Abbott B. P., et al., 2017, Phys. Rev. Lett., 119, 161101
  * Abbott et al. (2018) Abbott B. P., et al., 2018, Phys. Rev. Lett., 121, 161101
  * Abbott et al. (2019) Abbott B. P., et al., 2019, Physical Review X, 9, 011001
  * Abbott et al. (2020) Abbott B. P., et al., 2020, Living Reviews in Relativity, 23, 3
  * Abrahamyan et al. (2012) Abrahamyan S., et al., 2012, Phys. Rev. Lett., 108, 112502
  * Andersson & Comer (2001) Andersson N., Comer G. L., 2001, MNRAS, 328, 1129
  * Annala et al. (2018) Annala E., Gorda T., Kurkela A., Vuorinen A., 2018, Phys. Rev. Lett., 120, 172703
  * Balliet et al. (2020) Balliet L., Newton W., Cantu S., Budimir S., 2020, arXiv e-prints, p. arXiv:2009.07696
  * Barthelmy et al. (2005) Barthelmy S. D., et al., 2005, Space Sci. Rev., 120, 143
  * Bauswein et al. (2017) Bauswein A., Just O., Janka H.-T., Stergioulas N., 2017, ApJ, 850, L34
  * Baym et al. (1971) Baym G., Pethick C., Sutherland P., 1971, ApJ, 170, 299
  * Becker et al. (2018) Becker D., et al., 2018, European Physical Journal A, 54, 208
  * Blanchet (2006) Blanchet L., 2006, Living Reviews in Relativity, 9, 4
  * Brown (2000) Brown B. A., 2000, Phys. Rev. Lett., 85, 5296
  * Carson et al. (2019) Carson Z., Steiner A. W., Yagi K., 2019, Phys. Rev. D, 99, 043010
  * Chabanat et al. (1998) Chabanat E., Bonche P., Haensel P., Meyer J., Schaeffer R., 1998, Nuclear Phys. A, 635, 231
  * Chen et al. (2009) Chen L.-W., Cai B.-J., Ko C. M., Li B.-A., Shen C., Xu J., 2009, Phys. Rev. C, 80, 014322
  * Chen et al. (2010) Chen L.-W., Ko C. M., Li B.-A., Xu J., 2010, Phys. Rev. C, 82, 024321
  * Chugunov & Horowitz (2010) Chugunov A. I., Horowitz C. J., 2010, Monthly Notices of the Royal Astronomical Society: Letters, 407, L54
  * Cowling (1941) Cowling T. G., 1941, MNRAS, 101, 367
  * Cromartie et al. (2020) Cromartie H. T., et al., 2020, Nature Astronomy, 4, 72
  * Cutler & Flanagan (1994) Cutler C., Flanagan É. E., 1994, Phys. Rev. D, 49, 2658
  * D’Avanzo (2015) D’Avanzo P., 2015, Journal of High Energy Astrophysics, 7, 73
  * De et al. (2018) De S., Finstad D., Lattimer J. M., Brown D. A., Berger E., Biwer C. M., 2018, Phys. Rev. Lett., 121, 091102
  * Ducoin et al. (2011) Ducoin C., Margueron J., Providência C. m. c., Vidaña I., 2011, Phys. Rev. C, 83, 045810
  * Eichler et al. (1989) Eichler D., Livio M., Piran T., Schramm D. N., 1989, Nature, 340, 126
  * Essick et al. (2020) Essick R., Landry P., Holz D. E., 2020, Phys. Rev. D, 101, 063007
  * Fattoyev et al. (2018) Fattoyev F. J., Piekarewicz J., Horowitz C. J., 2018, Phys. Rev. Lett., 120, 172702
  * Fong et al. (2010) Fong W., Berger E., Fox D. B., 2010, ApJ, 708, 9
  * Gandolfi et al. (2012) Gandolfi S., Carlson J., Reddy S., 2012, Phys. Rev. C, 85, 032801
  * Gearheart et al. (2011) Gearheart M., Newton W. G., Hooker J., Li B.-A., 2011, MNRAS, 418, 2343
  * Goldstein et al. (2017) Goldstein A., et al., 2017, ApJ, 848, L14
  * Hebeler et al. (2010) Hebeler K., Lattimer J. M., Pethick C. J., Schwenk A., 2010, Physical Review Letters, 105, 161102
  * Holt & Kaiser (2017) Holt J. W., Kaiser N., 2017, Phys. Rev. C, 95, 034326
  * Holt & Lim (2018) Holt J. W., Lim Y., 2018, Physics Letters B, 784, 77
  * Horowitz & Piekarewicz (2001) Horowitz C. J., Piekarewicz J., 2001, Phys. Rev. Lett., 86, 5647
  * Horowitz et al. (2014) Horowitz C. J., et al., 2014, Journal of Physics G Nuclear Physics, 41, 093001
  * Kortelainen et al. (2010) Kortelainen M., Lesinski T., Moré J., Nazarewicz W., Sarich J., Schunck N., Stoitsov M. V., Wild S., 2010, Phys. Rev. C, 82, 024313
  * Kouveliotou et al. (1993) Kouveliotou C., Meegan C. A., Fishman G. J., Bhat N. P., Briggs M. S., Koshut T. M., Paciesas W. S., Pendleton G. N., 1993, ApJ, 413, L101
  * Lackey & Wade (2015) Lackey B. D., Wade L., 2015, Phys. Rev. D, 91, 043002
  * Lai (1994) Lai D., 1994, MNRAS, 270, 611
  * Landry & Essick (2019) Landry P., Essick R., 2019, Phys. Rev. D, 99, 084049
  * Lattimer & Lim (2013) Lattimer J. M., Lim Y., 2013, ApJ, 771, 51
  * Li et al. (2014) Li B.-A., Ramos À., Verde G., Vidaña I., 2014, European Physical Journal A, 50, 9
  * Lim & Holt (2018) Lim Y., Holt J. W., 2018, Phys. Rev. Lett., 121, 062701
  * Liu et al. (2010) Liu M., Wang N., Li Z.-X., Zhang F.-S., 2010, Phys. Rev. C, 82, 064306
  * Margalit & Metzger (2017) Margalit B., Metzger B. D., 2017, ApJ, 850, L19
  * McDermott et al. (1988) McDermott P. N., van Horn H. M., Hansen C. J., 1988, ApJ, 325, 725
  * Meegan et al. (2009) Meegan C., et al., 2009, ApJ, 702, 791
  * Most et al. (2018) Most E. R., Weih L. R., Rezzolla L., Schaffner-Bielich J., 2018, Phys. Rev. Lett., 120, 261103
  * Nazarewicz et al. (1996) Nazarewicz W., et al., 1996, Phys. Rev. C, 53, 740
  * Newton & Crocombe (2020) Newton W. G., Crocombe G., 2020, arXiv e-prints, p. arXiv:2008.00042
  * Newton et al. (2013) Newton W. G., Gearheart M., Li B.-A., 2013, ApJS, 204, 9
  * Newton et al. (2014) Newton W. G., Hooker J., Gearheart M., Murphy K., Wen D.-H., Fattoyev F. J., Li B.-A., 2014, European Physical Journal A, 50, 41
  * Newton et al. (2018) Newton W. G., Steiner A. W., Yagi K., 2018, ApJ, 856, 19
  * Oppenheimer & Volkoff (1939) Oppenheimer J. R., Volkoff G. M., 1939, Physical Review, 55, 374
  * Özel & Psaltis (2009) Özel F., Psaltis D., 2009, Phys. Rev. D, 80, 103003
  * Özel et al. (2010) Özel F., Baym G., Güver T., 2010, Phys. Rev. D, 82, 101301
  * Özel et al. (2016) Özel F., Psaltis D., Güver T., Baym G., Heinke C., Guillot S., 2016, ApJ, 820, 28
  * Pethick & Potekhin (1998) Pethick C., Potekhin A., 1998, Physics Letters B, 427, 7
  * Piekarewicz et al. (2012) Piekarewicz J., Agrawal B. K., Colò G., Nazarewicz W., Paar N., Reinhard P. G., Roca-Maza X., Vretenar D., 2012, Phys. Rev. C, 85, 041302
  * Raithel (2019) Raithel C. A., 2019, European Physical Journal A, 55, 80
  * Read et al. (2009a) Read J. S., Lackey B. D., Owen B. J., Friedman J. L., 2009a, Phys. Rev. D, 79, 124032
  * Read et al. (2009b) Read J. S., Markakis C., Shibata M., Uryū K., Creighton J. D. E., Friedman J. L., 2009b, Phys. Rev. D, 79, 124033
  * Read et al. (2013) Read J. S., et al., 2013, Phys. Rev. D, 88, 044042
  * Reed et al. (2021) Reed B. T., Fattoyev F. J., Horowitz C. J., Piekarewicz J., 2021, arXiv e-prints, p. arXiv:2101.03193
  * Reinhard & Flocard (1995) Reinhard P. G., Flocard H., 1995, Nuclear Phys. A, 584, 467
  * Reisenegger (1994) Reisenegger A., 1994, ApJ, 432, 296
  * Sotani et al. (2012) Sotani H., Nakazato K., Iida K., Oyamatsu K., 2012, Physical Review Letters, 108, 201101
  * Sotani et al. (2013) Sotani H., Nakazato K., Iida K., Oyamatsu K., 2013, MNRAS, 434, 2060
  * Steiner & Gandolfi (2012) Steiner A. W., Gandolfi S., 2012, Phys. Rev. Lett., 108, 081102
  * Steiner & Watts (2009) Steiner A. W., Watts A. L., 2009, Phys. Rev. Lett., 103, 181101
  * Steiner et al. (2008) Steiner A. W., Li B. A., Prakash M., 2008, in Exotic States of Nuclear Matter. pp 47–54 (arXiv:0711.4652), doi:10.1142/9789812797049_0008
  * Steiner et al. (2010) Steiner A. W., Lattimer J. M., Brown E. F., 2010, ApJ, 722, 33
  * Steiner et al. (2013) Steiner A. W., Lattimer J. M., Brown E. F., 2013, ApJ, 765, L5
  * Strohmayer et al. (1991) Strohmayer T., Ogata S., Iyetomi H., Ichimaru S., van Horn H. M., 1991, ApJ, 375, 679
  * Thiel et al. (2019) Thiel M., Sfienti C., Piekarewicz J., Horowitz C. J., Vanderhaeghen M., 2019, Journal of Physics G Nuclear Physics, 46, 093003
  * Tolman (1939) Tolman R. C., 1939, Physical Review, 55, 364
  * Trippa et al. (2008) Trippa L., Colò G., Vigezzi E., 2008, Phys. Rev. C, 77, 061304
  * Troja et al. (2010) Troja E., Rosswog S., Gehrels N., 2010, ApJ, 723, 1711
  * Tsang (2013) Tsang D., 2013, ApJ, 777, 103
  * Tsang et al. (2009) Tsang M. B., Zhang Y., Danielewicz P., Famiano M., Li Z., Lynch W. G., Steiner A. W., 2009, Phys. Rev. Lett., 102, 122701
  * Tsang et al. (2012a) Tsang M. B., et al., 2012a, Phys. Rev. C, 86, 015803
  * Tsang et al. (2012b) Tsang D., Read J. S., Hinderer T., Piro A. L., Bondarescu R., 2012b, Phys. Rev. Lett., 108, 011102
  * Tsang et al. (2019) Tsang M. B., Lynch W. G., Danielewicz P., Tsang C. Y., 2019, Physics Letters B, 795, 533
  * Yoshida & Lee (2002) Yoshida S., Lee U., 2002, A&A, 395, 201
  * Zhang & Li (2019) Zhang N.-B., Li B.-A., 2019, European Physical Journal A, 55, 39
  * Zhong et al. (2019) Zhong S.-Q., Dai Z.-G., Cheng J.-G., Lan L., Zhang H.-M., 2019, ApJ, 884, 25