# g-mode oscillations in hybrid stars:
A tale of two sounds

Prashanth Jaikumar<EMAIL_ADDRESS>Department of Physics and
Astronomy, California State University Long Beach, Long Beach, California
90840, USA  Alexandra Semposki<EMAIL_ADDRESS>Madappa Prakash
<EMAIL_ADDRESS>Department of Physics and Astronomy, Ohio University, Athens,
Ohio 45701, USA  Constantinos Constantinou<EMAIL_ADDRESS>INFN-
TIFPA, Trento Institute of Fundamental Physics and Applications, Povo, 38123
Trento, Italy  European Centre for Theoretical Studies in Nuclear Physics and
Related Areas, Villazzano, 38123 Trento, Italy

###### Abstract

We study the principal core $g$-mode oscillation in hybrid stars containing
quark matter and find that they have an unusually large frequency range
($\approx$ 200 - 600 Hz) compared to ordinary neutron stars or self-bound
quark stars of the same mass. Theoretical arguments and numerical calculations
that trace this effect to the difference in the behaviour of the equilibrium
and adiabatic sound speeds in the mixed phase of quarks and nucleons are
provided. We propose that the sensitivity of core $g$-mode oscillations to
non-nucleonic matter in neutron stars could be due to the presence of a mixed
quark-nucleon phase. Based on our analysis, we conclude that for binary
mergers where one or both components may be a hybrid star, the fraction of
tidal energy pumped into resonant $g$-modes in hybrid stars can exceed that of
a normal neutron star by a factor of 2-3, although resonance occurs during the
last stages of inspiral. A self-bound star, on the other hand, has a much
weaker tidal overlap with the $g$-mode. The cumulative tidal phase error in
hybrid stars, $\Delta\phi\cong$ 0.5 rad, is comparable to that from tides in
ordinary neutron stars, presenting a challenge in distinguishing between the
two cases. However, should the principal $g$-mode be excited to sufficient
amplitude for detection in a postmerger remnant with quark matter in its
interior, its frequency would be a possible indication for the existence of
non-nucleonic matter in neutron stars.

††preprint: APS/123-QED

##  I Introduction

The core of a neutron star (NS) can, in principle, support phases of dense
deconfined quark matter (QM) [1]. Confirmation of the presence of quarks in
NSs, however, has not been possible either through observations or from
lattice-gauge calculations of finite baryon density matter starting from the
Lagrangian of Quantum Chromodynamics (QCD). Although perturbative calculations
of QM have been performed [2, 3, 4], their applicability is limited to baryon
densities $n_{B}\gtrsim 40n_{s}$ [5], where $n_{s}\simeq 0.16~{}\rm{fm^{-3}}$
is the nuclear matter saturation density. Such densities, however, lie well
beyond the central densities $n_{c}$ in the range 3-8$n_{s}$ of observed NSs.
In view of this conundrum, theoretical studies of QM in NSs have been
exploratory in nature by positing either a sharp 1st-order or a smooth
crossover hadron-to-quark phase transition. Depending on the treatment of the
phase transition and the equations of state (EOSs) of hadronic and quark
matter, either a phase of pure QM or a phase in which hadrons are admixed with
quarks can be realized (for a detailed account, see Ref. [6] and an extensive
list of references therein). In either case, stars with quarks are difficult
to distinguish from normal NSs based on the knowledge of masses and radii
alone as similar results can be obtained with both. While the long-term
cooling of a NS can be affected by the presence of quarks, cooling data are
relatively sparse and gathered over decades [7, 8]. Gravitational wave
observations from compact binary mergers can be another probe of the EOS, but
currently, constraints on tidal polarizability [9, 10, 11] from gravitational
wave data [12] are consistent with both normal and quark-containing stars,
depending on the theoretical assumptions made [6, 13, 14].

In this paper, we are particularly interested in how NS oscillations can shed
light on the presence of QM in stars that contain an admixture of nucleons and
quarks (termed hybrid stars). Andersson and Kokkotas [15] have proposed that
NS oscillations (in particular, the $f,p$ modes) could be a “fingerprint” for
the supra-nuclear EOS in gravitational wave data. A review of potential
signatures of QM in NSs in the multi-messenger era, including the role of
their oscillations, can be found in [16]. Along these lines, we offer in this
work a new diagnostic of deconfined QM in NSs based on asteroseismology. We
show that a steep rise in the frequency of the principal $g$-mode (gravity
mode) occurs as soon as QM appears in a mixed phase in the core, exceeding the
typical core $g$-mode frequency of a nucleonic star by a factor of two or
more. This rise is essentially driven by a drop in the local equilibrium speed
of sound at the onset of the phase transition, while the adiabatic sound speed
changes only slightly. If this $g$-mode becomes resonant with the tidal force
during the late stages of a binary inspiral, the resulting energy transfer
from the orbital motion to the star via tidal coupling can affect the phase of
the gravitational waveform, and potentially signal a hybrid star.

NS oscillations are categorized by the nature of the dominant restoring force
for the perturbation in question. Several types of modes can be supported by a
star and it is desirable to investigate as many of them as possible in detail.
These modes are typically excited and sustained in different regions of the
star and their amplitudes and damping rates are subject to considerable
uncertainty. Here, for reasons that will become apparent, we focus our
attention on the $g$-mode and its coupling to gravitational waves.

A $g$-mode is a specific kind of non-radial fluid oscillation initiated when a
parcel of fluid is displaced against the background of a stratified
environment [17, 18]. While pressure equilibrium is rapidly restored via sound
waves, chemical equilibrium can take longer causing buoyancy forces to oppose
the displacement. Since cold NSs are not convective, the opposing force sets
up stable oscillations throughout the core, with a typical frequency, called
the (local) Brunt-Väisälä frequency [19], which depends on the difference
between the equilibrium and adiabatic sound speeds as well as the local metric
coefficients. Convectively stable $g$-modes exist for a wide range of models
of the EOS [20]. Though the $g$-mode in NSs has been studied before [21, 22,
23], with recent works incorporating additional features like hyperonic matter
and superfluidity [24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34], the novel
aspect of our work is the investigation of the $g$-mode frequency in a phase
transition from nuclear matter to a mixed phase of quarks and nucleons. We
point out that a similar result was obtained in [30] for superfluid hyperonic
stars. However, the calculations presented there were for a fixed NS mass of
about 1.64$M_{\odot}$ with a radius of 13.6 km. While their chosen hyperonic
EOS does support a maximum mass of $2.015M_{\odot}$ [35], the nuclear and
quark EOSs chosen in our work satisfy additional observational and
experimental constraints, and presents the effect for a wide range of masses
up to the observed maximum.

This paper also extends the results of [36] by incorporating aspects of
General Relativity in the fluid oscillation equations (while remaining within
the Cowling approximation), updating the nuclear EOS to include consistency
with radius constraints from tidal polarizability [37] and NICER data [38, 39]
on the radius of $\simeq 1.4M_{\odot}$ NS. We also provide new analytical
results for the two sound speeds in a mixed phase of quarks and nucleons. Our
study can be of practical interest to investigations of the sound speed in
NSs, which is attracting renewed attention [40, 41, 42, 43]. The matter of
detecting $g$-modes from the pre- or post-coalescence phase of binary NS
mergers is not addressed in detail here, but we present an estimate of its
impact on the accumulated tidal phase up to the merger.

It is pertinent to note that oscillation modes other than $g$-modes can also
potentially be affected by the presence of QM in NSs. Radial oscillation
modes, which however do not couple to gravitational waves, were studied in
[44]. Among non-radial modes, the $i$-mode (interface mode) has been recently
investigated for the special case of crystalline QM surrounded by a hadronic
envelope in [45] and its frequency can range from (300-1500) Hz, which can be
probed by current or third generation interferometric detectors. The $r$-mode
(Rossby mode) frequency and its damping rate for NSs containing QM also
differs from a purely nucleonic one [46, 47, 48]. The $s$-mode (shear mode)
can be excited in a mixed phase of quarks and nucleons and is sensitive to the
shear modulus of structures in the mixed phase [49], probing the surface
tension of QM. The $g$-mode oscillation in stars containing QM has been
studied in the case of a sharp interface [50, 51] between hadronic and quark
matter, yielding the spectrum of so-called discontinuity $g$-modes, but these
works assume a strong 1st-order phase transition and a large value of the
surface tension for QM, while we study the case of a mixed phase of
significant extent that would be favored if the same surface tension were
small enough111Here, we do not explicitly study surface and curvature effects
or the impact of a non-trivially structured mixed phase on the oscillation
spectrum [52, 49, 53], but a more complete treatment of $g$-modes in hybrid
stars should address these issues.. The $g$-mode for baryonic stars with
superfluid effects was studied in [24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34]
highlighting the subtle role of temperature and composition gradients in
driving the $g$-mode. In this work, we investigate the composition gradients
induced by the mixed phase of quarks and nucleons, which supports unusually
high-frequency $g$-modes through its effect on the adiabatic and equilibrium
sound speeds.

The organization of this paper is as follows. In Sec. II, we introduce the
governing equations of the $g$-mode and outline its relation to the two sound
speeds. In Sec. III, we present the EOS models in the nucleonic and quark
sectors chosen for our study of the $g$-mode spectrum. The rationale for our
parameter choices and the basic features of these models are highlighted here
for orientation. We stress that our choices are representative, but not
exhaustive. In Sec. IV, we derive expressions for the two sound speeds in the
mixed phase of quarks and nucleons. Results for the sound speeds, the Brunt-
Väisälä frequency and the $g$-mode frequency in nucleonic, quark and hybrid
stars are gathered and interpreted in Sec. V. The tidal forcing and phase
error due to $g$-mode excitation are estimated in Sec. VI. Our summary,
conclusions and outlook are contained in Sec. VII. The appendices contain
details about the determination of parameters in the nuclear EOS model, the
resulting NS structural properties and the two sound speeds.

##  II $g$-mode Oscillations

In this section, we outline the equations for fluid oscillations and non-
rotating NS structure that were used to determine the eigenfrequencies of the
$g$-mode. In general, the oscillatory displacement of a fluid element in a
spherically symmetric star is represented by a vector field
${\vec{\xi}}^{nlm}(\vec{r}){\rm e}^{-i\omega t}$ with $n,l$ and $m$ denoting
the radial, azimuthal and magnetic mode indices. To be precise, the
frequencies also carry subscripts $nlm$ implicitly understood, with
degeneracies that are broken in more realistic cases such as with rotation or
magnetic fields. For even-parity or spheroidal modes, separation into radial
and tangential components yields $\xi_{r}^{nlm}(\vec{r})$ =
$\eta_{r}^{nl}(r)Y_{lm}(\theta,\phi)$ and $\xi_{\perp}^{nlm}(\vec{r})$ =
$r\eta_{\perp}^{nl}(r)\nabla_{\perp}Y_{lm}(\theta,\phi)$, respectively, where
$Y_{lm}(\theta,\phi)$ are the spherical harmonics. From the perturbed
continuity equation for the fluid, the tangential function $\eta_{\perp}$ can
be traded for fluid variables as $\delta p/\epsilon$ =
$\omega^{2}r\eta_{\perp}(r)Y_{lm}(\theta,\phi){\rm e}^{-i\omega t}$, where
$\delta p$ is the corresponding local (Eulerian) pressure perturbation and
$\epsilon$ the local energy density. Within the relativistic Cowling
approximation222The Cowling approximation neglects the back reaction of the
gravitational potential and reduces the number of equations we have to solve.
While this approximation is not strictly consistent with our fully general
relativistic (GR) treatment of the equilibrium structure of the star, it does
not change our conclusions qualitatively or even quantitatively, since this
approximation is accurate for $g$-mode frequencies at the few % level [54].,
the equations of motion to be solved to determine the frequency of a
particular mode are [55, 17, 29]

$\displaystyle-\frac{1}{\mathrm{e}^{\lambda/2}r^{2}}\frac{\partial}{\partial
r}\left[\mathrm{e}^{\lambda/2}r^{2}\xi_{r}\right]+\frac{l(l+1)\mathrm{e}^{\nu}}{r^{2}\omega^{2}}\frac{\delta
p}{p+\epsilon}-\frac{\Delta p}{\gamma p}=0$ $\displaystyle\frac{\partial\delta
p}{\partial r}+g\left(1+\frac{1}{c_{\mathrm{s}}^{2}}\right)\delta
p+\mathrm{e}^{\lambda-\nu}h\left(N^{2}-\omega^{2}\right)\xi_{r}=0\,,$ (1)

where $h=p+\epsilon$, and we have suppressed the indices on $\omega$ and
$\xi$. Equation (1) involves thermodynamic quantities that follow from the
specific EOS. Specifically, $p$ denotes pressure, $\epsilon$ energy density,
and $\gamma$ the adiabatic index of the fluid. The Lagrangian variation of the
pressure enters as $\Delta p$, and is related to the Eulerian variation
$\delta p$ through the operator relation $\Delta\equiv\delta+\xi\cdot\nabla$.
The symbol $c_{s}$ denotes the adiabatic sound speed, which is related to the
adiabatic index as $c_{s}^{2}=\gamma p/(\mu_{n}n_{B})$ where $\mu_{n}$ is the
neutron chemical potential333In beta-equilibrated charge neutral matter, the
neutron chemical potential is sufficient to determine all other chemical
potentials. and $n_{B}$ the local baryon density. The equilibrium sound speed
enters through the Brunt-Väisälä frequency ($N$) which is given by

$\displaystyle N^{2}\equiv
g^{2}\Big{(}\frac{1}{c_{e}^{2}}-\frac{1}{c_{s}^{2}}\Big{)}{\rm
e}^{\nu-\lambda}\,,$ (2)

where $g=-\nabla\phi=-\nabla p/h$ with $h=\epsilon+p$ the enthalpy of the
fluid. Finally, $\nu(r)$ and $\lambda(r)$ are metric functions of the
unperturbed star which features in the Schwarzschild interior metric ($r<R$):

$\displaystyle-\mathrm{d}s^{2}\equiv\,g_{\alpha\beta}\mathrm{d}x^{\alpha}\mathrm{d}x^{\beta}=$
$\displaystyle-\mathrm{e}^{\nu(r)}\mathrm{d}t^{2}+\mathrm{e}^{\lambda(r)}\mathrm{d}r^{2}$
(3)
$\displaystyle+r^{2}\left(\mathrm{d}\theta^{2}+\sin^{2}\theta\mathrm{d}\varphi^{2}\right).$

Explicitly,

$\displaystyle e^{\lambda(r)}=\frac{1}{1-\left(\frac{2Gm(r)}{c^{2}r}\right)}$
(4)

and

$\displaystyle
e^{\nu(r)}=\exp\bigg{[}-\frac{2G}{c^{2}}\int_{0}^{r}\left(\frac{\left(m(r^{\prime})+\frac{4\pi
p(r^{\prime})r^{\prime
3}}{c^{2}}\right)}{r^{\prime}\left(r^{\prime}-\frac{2m(r^{\prime})G}{c^{2}}\right)}\right)dr^{\prime}\bigg{]}{\rm
e}^{\nu_{0}},$ (5)

where $m(r^{\prime})$ is the enclosed mass of the star at $r^{\prime}$. These
metric functions must match to their exterior values at the surface $r=R$,
hence the constant factor ${\rm e}^{\nu_{0}}$ [56].

In this work, we study the fundamental $g$-mode with $n$ = 1 and fix the
mode’s multipolarity at $l$ = 2. For the non-rotating stars we consider here,
solutions are degenerate in $m$. Note that our definition of the “fundamental”
mode refers to the lowest radial order of the $g$-mode which also has the
highest frequency. This should not be confused with the qualitatively
different $f$-mode which is also referred to sometimes as the fundamental
mode. Furthermore, overtones with lower frequency exist, but we do not perform
any computations with them here, since the fundamental $g$-mode has the
highest frequency and will be excited during the final stage of the pre-merger
phase when tidal forces are strongest. The system of equations in Eq. (1)
cannot be solved analytically even with a simple model of a neutron star. Our
aim will be to solve this numerically as an eigenvalue system for the $g$-mode
frequency $\omega$. Physically, the solution to this system of equations,
under the boundary conditions $\Delta p=0$ at the surface and
$\xi_{r},\,\delta p/\epsilon$ regular at the center, only exists for discrete
values of the mode frequency $\omega$. These values represent the $g$-mode
spectrum for a chosen stellar model. Because we have employed the Cowling
approximation and ignored the perturbations of the metric that must accompany
fluid perturbations, we cannot compute the imaginary part of the
eigenfrequency (damping time) of the $g$-mode444The damping time of $g$-modes
due to viscosity and gravitational wave emission, crudely estimated in [57,
36], suggests that the $g$-mode can become secularly unstable for temperatures
$10^{8}~{}{\rm K}<T<10^{9}~{}{\rm K}$ for rotational speeds exceeding twice
the $g$-mode frequency of a static star.. We turn now to discuss the EOS
models for nucleonic and quark matter employed in this work.

##  III Models for the Equation of State

The EOS models chosen in this work were predicated on the requirement that the
squared sound speeds $c_{e}^{2}$ (see Eq.(25)) and $c_{s}^{2}$ could be
calculated straightforwardly. In the nucleonic sector, we employ the model of
Zhao and Lattimer (ZL) [43] which is consistent with nuclear systematics at
and below the nuclear saturation density $n_{s}$. With suitable choices of the
slope of the nuclear symmetry energy at $n_{s}$ (see below), this EOS is also
consistent with the recent chiral effective theory calculations of Drischler
et al. [58] in which uncertainty estimates of the EOS up to $2n_{s}$ were
provided (see Fig. 2 of this reference). In addition, the ZL EOS is able to
support $\simeq 2M_{\odot}$ stars required by mass measurements of heavy NSs
[59], and is consistent with the recent radius measurements of $\sim
1.4M_{\odot}$ stars [38, 39] and the tidal deformability estimates from the
binary neutron star merger GW170817 [37].

Among the many models and treatments available in the quark sector [6], we
utilize the vMIT model of Gomes et al. [60] as a caricature of strongly
interacting quarks at the densities attained within NSs. Such interactions
between quarks are required to satisfy astrophysical data, particularly those
of heavy mass NSs. For the treatment of the nucleon-to-quark transition at
supra-nuclear densities, we employ the Gibbs construction [61] which renders
the transition to be smooth. Alternative models and treatments that feature
strong first- or second-order phase transitions will be undertaken in
subsequent work.

###  III.1 The ZL EOS for Nucleonic Matter

For completeness and to set the stage for the calculation of the two sound
speeds in the next section, relevant details of the ZL model are provided
below. The total energy density of interacting nucleons in neutron star matter
(NSM) is

$\displaystyle\epsilon_{B}=$
$\displaystyle\sum_{i=n,p}\frac{1}{\pi^{2}}\int_{0}^{k_{Fi}}k^{2}\sqrt{M_{B}^{2}+k^{2}}\,dk$
$\displaystyle+n_{B}V(n_{n},n_{p})\,,$ (6)

where the Fermi momenta $k_{Fi}=(3\pi^{2}n_{i})^{1/3}$ with $i=n,p$ and
$n_{B}=n_{n}+n_{p}$, and $M_{B}$ is the baryon mass in vacuum. In the ZL
model, interactions between nucleons are written as

$\displaystyle V(n_{n},n_{p})\equiv V(u,x)=$
$\displaystyle\,4x(1-x)(a_{0}u+b_{0}u^{\gamma})$
$\displaystyle+(1-2x)^{2}(a_{1}u+b_{1}u^{\gamma_{1}}),$ (7)

where $u=n_{B}/n_{s}$ and the proton fraction $x=n_{p}/n_{B}$. Adding and
subtracting $a_{0}u+b_{0}u^{\gamma}$, the above equation can be rewritten as

$\displaystyle V(u,x)$ $\displaystyle=V_{0}+S_{2i}(u)(1-2x)^{2}\quad$ with
$\displaystyle V_{0}$ $\displaystyle=a_{0}u+b_{0}u^{\gamma},\quad$
$\displaystyle\quad S_{2i}(u)$
$\displaystyle=(a_{1}-a_{0})u+b_{1}u^{\gamma_{1}}-b_{0}u^{\gamma}\,,$ (8)

where the subscript “$2i$” in $S_{2i}$ refers to the interacting part of the
total symmetry energy $S_{2}=S_{2k}+S_{2i}$, with $S_{2k}$ representing the
kinetic part. Expanding the kinetic part in Eq. (6) to order $(1-2x)^{2}$, we
obtain the result555For the derivation of the kinetic part of the symmetry
energy and its derivatives, see the Appendix.

$\displaystyle
S_{2k}=\frac{1}{8}\left[\frac{1}{n}\frac{\partial^{2}\epsilon_{Bk}}{\partial
x^{2}}\right]_{x=\frac{1}{2}}=\frac{k_{F}^{2}}{6E_{F}}\,,$ (9)

where $k_{F}=(3\pi^{2}n_{B}/2)^{1/3}$ is the Fermi wave number of symmetric
nuclear matter (SNM) and $E_{F}=\sqrt{k_{F}^{2}+M_{B}^{2}}$. Collecting the
results, the energy per baryon relative to $M_{B}$ is given by

$\displaystyle\frac{\epsilon_{B}}{n_{B}}-M_{B}$ $\displaystyle=$
$\displaystyle E(u,x)=E_{{\rm SNM}}+S_{2}(u)(1-2x)^{2}$ where $\displaystyle
E_{{\rm SNM}}$ $\displaystyle=$ $\displaystyle
T_{1/2}+V_{1/2}=T_{1/2}+(a_{0}u+b_{0}u^{\gamma}),$ $\displaystyle S_{2}(u)$
$\displaystyle=$
$\displaystyle\frac{k_{F}^{2}}{6E_{F}}+(a_{1}-a_{0})u+b_{1}u^{\gamma_{1}}-b_{0}u^{\gamma}\,.$
(10)

The kinetic energy per baryon $T_{1/2}$ in SNM ($x=1/2$) is given by the
expression

$\displaystyle T_{1/2}$ $\displaystyle=$
$\displaystyle\frac{\epsilon_{1/2}^{\rm kin}}{n_{B}}-M_{B}\quad{\rm with}$
$\displaystyle\epsilon_{1/2}^{\rm kin}$ $\displaystyle=$
$\displaystyle\frac{2}{4\pi^{2}}\bigg{[}k_{F}E_{F}\left(k_{F}^{2}+\frac{M_{B}^{2}}{2}\right)$
(11) $\displaystyle-$
$\displaystyle\frac{1}{2}M_{B}^{4}\ln\left(\frac{k_{F}+E_{F}}{M_{B}}\right)\bigg{]}\,,$

where $n_{B}$, $k_{F}$ and $E_{F}$ refer to those in SNM.

The baryon pressure $p_{B}$ is

$\displaystyle p_{B}=n_{s}u^{2}\frac{dE_{B}}{du}=p_{\rm
SNM}+n_{s}u^{2}(1-2x)^{2}\frac{dS_{2}(u)}{du}\,,$ (12)

where

$\displaystyle p_{SNM}$ $\displaystyle=$ $\displaystyle p_{1/2}^{\rm
kin}+n_{s}~{}(a_{0}u^{2}+\gamma b_{0}u^{\gamma+1}),\quad{\rm with}$
$\displaystyle p_{1/2}^{\rm kin}$ $\displaystyle=$
$\displaystyle\frac{2}{12\pi^{2}}\bigg{[}k_{F}E_{F}\left(k_{F}^{2}-\frac{3}{2}M_{B}^{2}\right)$
$\displaystyle+$
$\displaystyle\frac{3}{2}M_{B}^{4}\ln\left(\frac{k_{F}+E_{F}}{M_{B}}\right)\bigg{]},\quad{\rm
and}$ $\displaystyle u\frac{dS_{2}(u)}{du}$ $\displaystyle=$
$\displaystyle\frac{2}{3}S_{2k}\left[1-18\left(\frac{S_{2k}}{k_{F}}\right)^{2}\right]+(a_{1}-a_{0})u$
(13) $\displaystyle+$ $\displaystyle b_{1}\gamma_{1}u^{\gamma_{1}}-b_{0}\gamma
u^{\gamma}.$

The incompressibility $K_{B}$ in SNM is obtained from

$\displaystyle K_{B}$ $\displaystyle=$ $\displaystyle
9\frac{dp_{B}}{dn_{B}}=9\left[2u\frac{dE_{B}}{du}+u^{2}\frac{d^{2}E_{B}}{du^{2}}\right]$
(14) $\displaystyle=$ $\displaystyle
9\left\\{\frac{k_{F}^{2}}{3E_{F}}+\left[2a_{0}u+\gamma(\gamma+1)b_{0}u^{\gamma}\right]\right\\}$

The energy per baryon in pure neutron matter (PNM in which $x=0$) relative to
the baryon mass is

$\displaystyle E_{{\rm PNM}}$ $\displaystyle=$ $\displaystyle
T_{0}+V_{0}=T_{0}+(a_{1}u+b_{1}u^{\gamma}_{1})$ $\displaystyle T_{0}$
$\displaystyle=$ $\displaystyle\frac{\epsilon_{0}^{\rm
kin}}{n_{B}}-M_{B}\quad{\rm with}$ $\displaystyle\epsilon_{0}^{\rm kin}$
$\displaystyle=$
$\displaystyle\frac{1}{4\pi^{2}}\bigg{[}k_{Fn}E_{Fn}\left(k_{Fn}^{2}+\frac{M_{B}^{2}}{2}\right)$
(15) $\displaystyle-$
$\displaystyle\frac{1}{2}M_{B}^{4}\ln\left(\frac{k_{Fn}+E_{Fn}}{M_{B}}\right)\bigg{]}\,,$

where now $n_{B}=n_{n}$, $k_{Fn}=(3\pi^{2}n_{n})^{1/3}=2^{1/3}k_{F}$, and
$E_{Fn}=\sqrt{k_{Fn}^{2}+M_{B}^{2}}$.

The determination of the EOS constants in SNM and PNM, and relevant NS
structural properties are summarized in Appendix A.

###  III.2 The vMIT Equation of State for Quark Matter

In recent years, variations of the original bag model [62] have been adopted
[60, 63] to calculate the structure of NSs with quarks in their cores to
account for $\geq 2\mathrm{M}_{\odot}$ maximum-mass stars. Termed as vMIT or
vBag models, the QCD perturbative results are dropped and replaced by
repulsive vector interactions between quarks in such works. We will provide
some numerical examples of the vMIT model for contrast with other models as
those of the vBag model turn out to be qualitatively similar.

The Lagrangian density of the vMIT bag model is

$\mathcal{L}=\sum_{i}\left[\bar{\psi}_{i}\left(i\not{\partial}-m_{i}-B\right)\psi_{i}+\mathcal{L}_{\mathrm{int}}\right]\Theta,$
(16)

where
$\mathcal{L}_{\mathrm{int}}=\mathcal{L}_{\mathrm{pert}}+\mathcal{L}_{\mathrm{vec}}$
describes quarks of mass $m_{i}$ confined within a bag as denoted by the
$\Theta$ function. For three flavors $i=u,d,s$ and three colors $N_{c}=3$ of
quarks, the number and baryon densities, energy density, pressure and chemical
potentials in the bag model are given by

$\displaystyle n_{i}$
$\displaystyle=2N_{c}\int^{k_{Fi}}\frac{d^{3}k}{(2\pi)^{3}},\quad
n_{B}=\frac{1}{3}\sum_{i}n_{i}$ (17) $\displaystyle\epsilon_{Q}$
$\displaystyle=2N_{c}\sum_{i}\int^{k_{Fi}}\frac{d^{3}k}{(2\pi)^{3}}\sqrt{k^{2}+m_{i}^{2}}+\epsilon_{\mathrm{int}}+B$
(18) $\displaystyle p_{Q}$
$\displaystyle=\frac{2N_{c}}{3}\sum_{i}\int^{k_{Fi}}\frac{d^{3}k}{(2\pi)^{3}}\frac{k^{2}}{\sqrt{k^{2}+m_{i}^{2}}}+p_{\mathrm{int}}-B$
(19) $\displaystyle\mu_{i}$
$\displaystyle=\sqrt{k_{Fi}^{2}+m_{i}^{2}}+\mu_{\mathrm{int},i}$ (20)

The upper limit of the integrals $k_{Fi}$ is the Fermi wave number for each
species $i$, which, at zero temperature, appropriately terminates the
integration over $k$. The first terms in $\epsilon_{Q}$ and in $p_{Q}$ are
free Fermi gas (FG) contributions, $\epsilon_{\mathrm{FG}}$ and
$p_{\mathrm{FG}}$ respectively, the second terms are due to
$\mathcal{L}_{\text{int }}$ and $B$ is the bag constant that accounts for the
cost of confining the quarks inside a bag. The $m_{i}$ are quark masses,
generally taken to be current quark masses. The $u$ and $d$ quark masses are
commonly set to zero (because at high density, $k_{Fi}$ in these cases far
exceed $m_{i}$), whereas that of the $s$ quark is taken at its Particle Data
Group (PDG) value. The QCD perturbative calculations of $\epsilon_{\text{pert
}}$ and $p_{\text{pert }}$, and the ensuing results for the structure of NSs
containing quarks within the cores as well as self-bound strange quark stars
are discussed in [5]. At leading order of QCD corrections, the results are
qualitatively similar to what one obtains by simply using the FG results with
an appropriately chosen value of $B$. As results of perturbative calculations
are deemed to be valid only for $n_{B}\geq 40n_{s}$, they are dropped in the
vMIT model. The Lagrangian density from vector interactions

$\quad\mathcal{L}_{\text{vec
}}=-G_{v}\sum_{i}\bar{\psi}\gamma_{\mu}V^{\mu}\psi+\left(m_{V}^{2}/2\right)V_{\mu}V^{\mu}\,,$
(21)

where interactions among the quarks occur via the exchange of a vector-
isoscalar meson $V^{\mu}$ of mass $m_{V},$ is chosen in Ref. [54]. Explicitly,

$\displaystyle\epsilon_{Q}$
$\displaystyle=\sum_{i}\epsilon_{\mathrm{FG},\mathrm{i}}+\frac{1}{2}\left(\frac{G_{v}}{m_{V}}\right)^{2}n_{Q}^{2}+B$
(22) $\displaystyle p_{Q}$
$\displaystyle=\sum_{i}p_{\mathrm{FG},\mathrm{i}}+\frac{1}{2}\left(\frac{G_{v}}{m_{V}}\right)^{2}n_{Q}^{2}-B$
(23) $\displaystyle\mu_{i}$
$\displaystyle=\sqrt{k_{Fi}^{2}+m_{i}^{2}}+\left(\frac{G_{v}}{m_{V}}\right)^{2}n_{Q}\,,$
(24)

where $n_{Q}=\sum_{i}n_{i},$ and the bag constant $B$ is chosen appropriately
to enable a transition to matter containing quarks. Note that terms associated
with the vector interaction above are similar to those in hadronic models. We
studied model parameters in a wide range $B^{1/4}=(155-180)~{}\mathrm{MeV}$
and $a=\left(G_{v}/m_{V}\right)^{2}=(0.1-0.3)~{}\mathrm{fm}^{2}$ and report
results for specific values within this range.

##  IV Sound Speeds in the Pure and Mixed Phases

As the difference of the adiabatic and equilibrium sound speeds drives the
restoring force for $g$-modes, it is instructive to collect some general
expressions for these two sound speeds in the pure and mixed phases. For the
pure phase of $npe$ matter, these expressions are derived and applied in [20],
but given their central role in this work, and the fact that we also extend
the application to $npe\mu$ and quark matter, we detail their derivation below
for completeness. For the mixed phase, we derive expressions that have not, to
our knowledge, been previously reported in the literature. First, a point of
notation: the equilibrium squared sound speed is commonly defined in the
literature [20, 64, 6] by the symbol $c_{s}^{2}$, which we reserve here for
the squared adiabatic sound speed, as in [29]. The equilibrium sound speed is
defined by

$\displaystyle c_{e}^{2}=\frac{dp}{d\epsilon}\,,$ (25)

where $p$ and $\epsilon$ are the total pressure and energy density in matter

$\displaystyle\epsilon=\epsilon_{B}\,+\sum\limits_{l=e^{-},\,\mu^{-}}\epsilon_{l}\,,\quad
p=p_{B}\,+\sum\limits_{l=e^{-},\,\mu^{-}}p_{l}\,.$ (26)

In Eq.(26), the leptonic energies are the $T$=0 degenerate Fermi gas
expressions for massive leptons. Being a total derivative, the derivative is
taken along the curve satisfying both mechanical and chemical equilibrium,
i.e., $\beta$-equilibrium conditions hold. In NSM, when only $npe$ are present
in equilibrium, the composition at fixed baryon density ($n_{n}+n_{p}$) is
completely fixed once the proton fraction $x_{p}$ (=$x_{e}$ by charge
neutrality) is determined. In this case, the squared adiabatic sound speed is
defined as

$\displaystyle c_{s}^{2}=\left(\frac{\partial
p}{\partial\epsilon}\right)_{x}\,,$ (27)

where $x=x_{p}=x_{e}$. In the partial derivative, the composition is held
fixed, i.e., $\beta$-equilibrium conditions are imposed only after all
derivatives have been evaluated. The resulting distinction between these two
speeds plays an important role in determining the oscillation frequencies of
non-radial oscillations such as $g$-modes:

$\displaystyle\omega^{2}\propto\left(\frac{1}{c_{e}^{2}}-\frac{1}{c_{s}^{2}}\right)=\frac{(c_{s}^{2}-c_{e}^{2})}{c_{e}^{2}c_{s}^{2}}\,.$
(28)

Note that both the above speeds are dependent on density which varies over a
large range in NSM. Furthermore, an individual knowledge of both speeds is
required. In what follows, we apply Eqs. (25) and (27) to the case of a pure
and mixed phase.

###  IV.1 The Pure Phase

####  IV.1.1  Sound speeds in $npe$ matter

It is useful to recast the general expressions Eqs. (25) and (27) in terms of
derivatives of the individual chemical potentials with respect to density,
since such expressions are amenable to both analytical and numerical checks.
Without loss of generality, we have

$\displaystyle
c_{e}^{2}=\frac{dp}{d\epsilon}=\left(\frac{dp}{dn_{B}}\right)\bigg{/}\left(\frac{d\epsilon}{dn_{B}}\right)\,,$
(29)

Considering $npe$ matter as an example, differentiating the total energy
energy density inclusive of electrons

$\displaystyle\epsilon(n_{B},x)=n_{B}[M_{B}+E(n_{B},x)]$ (30)

with respect to $n_{B}$, we have

$\displaystyle\left(\frac{d\epsilon}{dn_{B}}\right)=\frac{\epsilon}{n_{B}}+n_{B}\left(\frac{dE}{dn_{B}}\right)\,,$
(31)

where $E(n_{B},x)$ is the energy per baryon. The second term on the right hand
side of Eq. (31) becomes

$\displaystyle\left(\frac{dE}{dn_{B}}\right)=\left(\frac{\partial E}{\partial
n_{B}}\right)+\left(\frac{\partial E}{\partial
x}\right)_{n_{B}}\left(\frac{dx}{dn_{B}}\right)\,.$ (32)

For the equilibrium sound speed, the $\beta$-equilibrium condition
$\left(\frac{\partial E}{\partial x}\right)_{n_{B}}=0$ yields

$\displaystyle\left(\frac{dE}{dn_{B}}\right)=\left(\frac{\partial E}{\partial
n_{B}}\right)\,.$ (33)

Thus,

$\displaystyle\left(\frac{d\epsilon}{dn_{B}}\right)=\frac{\epsilon}{n_{B}}+\frac{1}{n_{B}}\left[n_{B}^{2}\left(\frac{dE}{dn_{B}}\right)\right]=\frac{(\epsilon+p)}{n_{B}}\,.$
(34)

From the thermodynamic identity, using charge neutrality ($x=x_{e}$) and beta-
equilibrium,

$\displaystyle\epsilon+p$ $\displaystyle=$
$\displaystyle\mu_{n}n_{n}+\mu_{p}n_{p}+\mu_{e}n_{e}$ (35) $\displaystyle=$
$\displaystyle\mu_{n}n_{B}-(\mu_{n}-\mu_{p}-\mu_{e})n_{B}=\mu_{n}n_{B}\,,$

leading to the simple result

$\displaystyle\left(\frac{d\epsilon}{dn_{B}}\right)=\left(\frac{\partial\epsilon}{\partial
n_{B}}\right)_{x}=\mu_{n}$ (36)

This implies that

$\displaystyle c_{e}^{2}$ $\displaystyle\equiv$
$\displaystyle\frac{dp}{d\epsilon}=\frac{1}{\mu_{n}}\left(\frac{dp}{dn_{B}}\right)=\frac{1}{\mu_{n}}n_{B}\left(\frac{d\mu_{n}}{dn_{B}}\right)\,$
(37) $\displaystyle=$ $\displaystyle\left(\frac{d\ln\mu_{n}}{d\ln
n_{B}}\right)\,,$

where we have again taken advantage of the thermodynamic identity to relate
the required derivative of $p$ to that of $\mu_{n}$.

The adiabatic squared sound speed can be expressed as

$\displaystyle c_{s}^{2}$ $\displaystyle=$ $\displaystyle\left(\frac{\partial
p}{\partial\epsilon}\right)_{x}=\left(\frac{\partial p}{\partial
n_{B}}\right)_{x}\bigg{/}\left(\frac{\partial\epsilon}{\partial
n_{B}}\right)_{x}\,$ (38) $\displaystyle=$
$\displaystyle\frac{n_{B}}{\epsilon+p}\left(\frac{\partial p}{\partial
n_{B}}\right)_{x}=\frac{1}{\mu_{\rm avg}}\left(\frac{\partial p}{\partial
n_{B}}\right)_{x}\,,$ (39)

where we have used the equality in Eq. (34), which is valid at constant
composition even in the absence of $\beta$-equilibrium, and introduced an
average chemical potential $\mu_{\rm
avg}=(\sum_{i}\mu_{i}n_{i})/n_{B}=(\epsilon+p)/n_{B}$. Since $p=p_{B}+p_{e}$,

$c_{s}^{2}=\frac{1}{(\epsilon+p)}\left[\left(u\frac{\partial p_{B}}{\partial
u}\right)_{x}+\left(u\frac{\partial p_{e}}{\partial u}\right)_{x}\right]\,.$
(40)

The required derivatives are analytic:

$\displaystyle\left(u\frac{\partial p_{B}}{\partial u}\right)_{x}$
$\displaystyle=$
$\displaystyle\frac{2}{9\pi^{2}}\frac{k_{F}^{5}}{E_{F}}+n_{s}u\Big{[}2a_{0}u+b_{0}\gamma(\gamma+1)u^{\gamma}\Big{]}$
(41) $\displaystyle+$ $\displaystyle
n_{s}u(1-2x)^{2}\left\\{\frac{2k_{F}^{2}}{27E_{F}}\left(1-\frac{9k_{F}^{2}}{10E_{F}^{2}}+\frac{3k_{F}^{4}}{10E_{F}^{4}}\right)\right.$
$\displaystyle+$
$\displaystyle\left.\bigg{[}2u(a_{1}-a_{0})+b_{1}\gamma_{1}(\gamma_{1}+1)u^{\gamma_{1}}\right.$
$\displaystyle-$
$\displaystyle\left.b_{0}\gamma(\gamma+1)u^{\gamma}\bigg{]}\right\\}$
$\displaystyle\left(u\frac{\partial p_{e}}{\partial u}\right)_{x}$
$\displaystyle=$ $\displaystyle\frac{1}{3}n_{e}\mu_{e}\,.$ (42)

Thus, the difference of squared sound speeds becomes

$\displaystyle c_{s}^{2}-c_{e}^{2}=\frac{1}{\mu_{\rm avg}}\left(\frac{\partial
p}{\partial
n_{B}}\right)_{x}-\frac{1}{\mu_{n}}\left(\frac{dp}{dn_{B}}\right)\,.$ (43)

At this point all the necessary ingredients for the calculation of the speed-
of-sound difference are present. It is instructive, however, to obtain a
complementary expression in which its physical causes, namely
$\beta$-equilibrium and compositional gradients, are made explicit. To that
end, we proceed as follows: Noting that

$\displaystyle\frac{dp}{dn_{B}}=\left(\frac{\partial p}{\partial
n_{B}}\right)_{x}+\left(\frac{\partial p}{\partial
x}\right)_{n_{B}}\frac{dx}{dn_{B}}\,,$ (44)

Eq. (43) can be recast as

$\displaystyle c_{s}^{2}-c_{e}^{2}$ $\displaystyle=$
$\displaystyle\left(\frac{1}{\mu_{\rm
avg}}-\frac{1}{\mu_{n}}\right)\left(\frac{\partial p}{\partial
n_{B}}\right)_{x}-\frac{1}{\mu_{n}}\left(\frac{\partial p}{\partial
x}\right)_{n_{B}}\frac{dx}{dn_{B}}$ (45) $\displaystyle=$
$\displaystyle-~{}\frac{x\tilde{\mu}}{\mu_{\rm
avg}\mu_{n}}\left(\frac{\partial p}{\partial
n_{B}}\right)_{x}-\frac{1}{\mu_{n}}\left(\frac{\partial p}{\partial
x}\right)_{n_{B}}\frac{dx}{dn_{B}}\,$

where $\tilde{\mu}=\mu_{e}+\mu_{p}-\mu_{n}$. Anticipating that
$\beta$-equilibrium will be imposed at the end, we note that the first term
above vanishes as $\tilde{\mu}=0$, which leads to

$\displaystyle c_{s}^{2}-c_{e}^{2}=-\frac{1}{\mu_{n}}\left(\frac{\partial
p}{\partial x}\right)_{n_{B}}\frac{dx}{dn_{B}}\,.$ (46)

Utilizing $p=n_{B}^{2}\frac{\partial E}{\partial n_{B}}$ and interchanging the
order of derivatives

$\displaystyle
c_{s}^{2}-c_{e}^{2}=-\frac{1}{\mu_{n}}n_{B}^{2}\frac{\partial}{\partial
n_{B}}\left(\frac{\partial E}{\partial x}\right)_{n_{B}}\frac{dx}{dn_{B}}\,,$
(47)

which can be further rewritten as 666Observe the interesting relation
$\frac{\partial p}{\partial x}=n_{B}^{2}\frac{\partial\tilde{\mu}}{\partial
n_{B}}$, noted also in the context of bulk viscosity studies [65].

$\displaystyle
c_{s}^{2}-c_{e}^{2}=-\frac{1}{\mu_{n}}n_{B}^{2}\left(\frac{\partial\tilde{\mu}}{\partial
n_{B}}\right)_{x}\frac{dx}{dn_{B}}\,.$ (48)

It remains now to determine $\frac{dx}{dn_{B}}$. As

$\displaystyle d\tilde{\mu}=\left(\frac{\partial\tilde{\mu}}{\partial
n_{B}}\right)_{x}dn_{B}+\left(\frac{\partial\tilde{\mu}}{\partial
x}\right)_{n_{B}}dx=0\,,$ (49)
$\displaystyle\frac{dx}{dn_{B}}=-\left(\frac{\partial\tilde{\mu}}{\partial
n_{B}}\right)_{x}\bigg{/}\left(\frac{\partial\tilde{\mu}}{\partial
x}\right)_{n_{B}}\,.$ (50)

With this relation, Eq. (48) becomes

$\displaystyle c_{s}^{2}$ $\displaystyle=$ $\displaystyle
c_{e}^{2}+\frac{\left[n_{B}\left(\frac{\partial\tilde{\mu}}{\partial
n_{B}}\right)_{x}\right]^{2}}{\mu_{n}\left(\frac{\partial\tilde{\mu}}{\partial
x}\right)_{n_{B}}}\,,$ (51)

which illustrates the influence of the density and compositional gradients of
the two sound speeds. Thus far, we have simply retraced the steps originally
given in [20] (Sec. 4.2 of this reference). In $npe$ matter under the
constraint of charge neutrality, the independent variables chosen are $n_{B}$
and $x$, and thus a partial derivative of $\tilde{\mu}$ with respect to
$n_{B}$ ($x$) implies that $x$ ($n_{B}$) is fixed.

Casting the expressions for the sound speeds in terms of the chemical
potentials is expedient, as illustrated below for the case of $npe$ matter.
Note that the average chemical potential $\mu_{\rm avg}=\mu_{n}$ only in
$\beta$-equilibrium. At fixed $x$, with the relation
$\tilde{\mu}=\mu_{e}-\hat{\mu}$,

$\displaystyle n_{B}\frac{\partial\tilde{\mu}}{\partial n_{B}}$
$\displaystyle=$ $\displaystyle u\frac{\partial(\mu_{e}-\hat{\mu})}{\partial
u}\,.$ (52)

For the term in Eq. (52) involving electrons, we have

$\displaystyle\mu_{e}=\hbar c~{}(3\pi^{2}n_{s}u)^{1/3}x^{1/3}\quad{\rm
and}\quad u\frac{\partial\mu_{e}}{\partial u}=\frac{\mu_{e}}{3}\,,$ (53)

while for baryons,777Details of the derivatives of the kinetic part of the
symmetry energy are given in Appendix A.

$\displaystyle\hat{\mu}$ $\displaystyle=$
$\displaystyle\mu_{n}-\mu_{p}=4S_{2}(u)(1-2x)\quad{\rm with}$ $\displaystyle
S_{2}(u)$ $\displaystyle=$ $\displaystyle
S_{2k}+S_{2i}=\frac{k_{F}^{2}}{6E_{F}}$ $\displaystyle+$
$\displaystyle(a_{1}-a_{0})u+b_{1}u^{\gamma_{1}}-b_{0}u^{\gamma}\,,$
$\displaystyle uS_{2k}^{\prime}$ $\displaystyle=$
$\displaystyle\frac{1}{3}\cdot
2S_{2k}\left[1-18\left(\frac{S_{2k}}{k_{F}}\right)^{2}\right],$
$\displaystyle{\rm and}\quad uS_{2i}^{\prime}$ $\displaystyle=$
$\displaystyle(a_{1}-a_{0})u+\gamma_{1}b_{1}u^{\gamma_{1}}-\gamma
b_{0}u^{\gamma}\,,$ $\displaystyle uS_{2}^{\prime}$ $\displaystyle=$
$\displaystyle uS_{2k}^{\prime}+uS_{2i}^{\prime}.$ (54)

Putting together the above results, we have

$\displaystyle n_{B}\frac{\partial\tilde{\mu}}{\partial
n_{B}}=\frac{\mu_{e}}{3}-4(1-2x)~{}uS_{2}^{\prime}\,.$ (55)

Derivatives with respect to $x$ of $\tilde{\mu}$ at fixed density are also
straightforward. For

$\displaystyle\frac{\partial\tilde{\mu}}{\partial
x}=\frac{\partial(\mu_{e}-\hat{\mu})}{\partial x}\,,$ (56)

we note that

$\displaystyle\frac{\partial\mu_{e}}{\partial x}$ $\displaystyle=$
$\displaystyle\frac{1}{3}\frac{\mu_{e}}{x}\quad{\rm
and}\quad\frac{\partial\hat{\mu}}{\partial x}=-8S_{2}(u)\quad{\rm so~{}that}$
$\displaystyle\frac{\partial\tilde{\mu}}{\partial x}$ $\displaystyle=$
$\displaystyle\frac{1}{3}\frac{\mu_{e}}{x}+8S_{2}(u)\,.$ (57)

The equivalence of Eqs. (40) and (51) is established analytically in Appendix
A.5.

####  IV.1.2 Sound speeds in $npe\mu$ matter

Going beyond the results in [20], one way to include muons is by choosing
$n_{B}$, $x=x_{e}+x_{\mu}$ and $x_{\mu}\equiv y$ as the independent variables.
The formal expression for the squared adiabatic sound speed remains the same
as in $npe$ matter, i.e., Eq. (40) but now $\left(u~{}\partial p_{e}/\partial
u\right)_{x}$ [Eq. (42)] is replaced by

$\left(u\frac{\partial p_{lep}}{\partial
u}\right)_{x,x_{\mu}}=\frac{1}{3}n_{e}\mu_{e}+\frac{1}{3}n_{\mu}\left(\frac{\mu_{\mu}^{2}-m_{\mu}^{2}}{\mu_{\mu}}\right)\,.$
(58)

where $lep=e^{-},\mu^{-}$.

Furthermore, by retracing the steps leading to Eq. (51), its $npe\mu$
equivalent is obtained as

$c_{s}^{2}-c_{e}^{2}=-\frac{1}{\mu_{n}}\left(\left.\frac{\partial P}{\partial
x}\right|_{n_{B},y}\frac{dx}{dn_{B}}+\left.\frac{\partial P}{\partial
y}\right|_{n_{B},x}\frac{dy}{dn_{B}}\right)$ (59)

with 888The intermediate steps leading to Eqs. (60)-(61) are detailed in
Appendix. A.6

$\displaystyle\frac{dx}{dn_{B}}$ $\displaystyle=$
$\displaystyle\frac{\left.\frac{\partial\tilde{\mu}_{x}}{\partial
y}\right|_{n_{B},x}\left.\frac{\partial\tilde{\mu}_{y}}{\partial
n_{B}}\right|_{x,y}-\left.\frac{\partial\tilde{\mu}_{y}}{\partial
y}\right|_{n_{B},x}\left.\frac{\partial\tilde{\mu}_{x}}{\partial
n_{B}}\right|_{x,y}}{\left.\frac{\partial\tilde{\mu}_{x}}{\partial
y}\right|_{n_{B},x}\left.\frac{\partial\tilde{\mu}_{y}}{\partial
x}\right|_{n_{B},y}-\left.\frac{\partial\tilde{\mu}_{y}}{\partial
y}\right|_{n_{B},x}\left.\frac{\partial\tilde{\mu}_{x}}{\partial
x}\right|_{n_{B},y}}$ (60) $\displaystyle\frac{dy}{dn_{B}}$ $\displaystyle=$
$\displaystyle\frac{\left.\frac{\partial\tilde{\mu}_{x}}{\partial
x}\right|_{n_{B},y}\left.\frac{\partial\tilde{\mu}_{y}}{\partial
n_{B}}\right|_{x,y}-\left.\frac{\partial\tilde{\mu}_{y}}{\partial
x}\right|_{n_{B},y}\left.\frac{\partial\tilde{\mu}_{x}}{\partial
n_{B}}\right|_{x,y}}{\left.\frac{\partial\tilde{\mu}_{x}}{\partial
x}\right|_{n_{B},y}\left.\frac{\partial\tilde{\mu}_{y}}{\partial
y}\right|_{n_{B},x}-\left.\frac{\partial\tilde{\mu}_{y}}{\partial
x}\right|_{n_{B},y}\left.\frac{\partial\tilde{\mu}_{x}}{\partial
y}\right|_{n_{B},x}}~{}.$ (61)

The chemical potentials $\tilde{\mu}_{x}=\mu_{p}+\mu_{e}-\mu_{n}$ and
$\tilde{\mu}_{y}=\mu_{\mu}-\mu_{e}$ are zero in $\beta$-equilibrated matter.
Equations (59)-(61), while demonstrating that compositional gradients are at
the core of g-mode oscillations, are lengthy and computationally more involved
compared to the direct calculation of the adiabatic sound speed in $npe\mu$
matter using Eqs. (41) and (58). For the sake of completeness, we provide here
the explicit expressions for the adiabatic sound speed in $npe\mu$ matter
arising from (59)-(61) which are in excellent numerical agreement with the
more direct method:

$\displaystyle
c_{s}^{2}=c_{e}^{2}+\frac{1}{\mu_{n}}\Big{(}T_{1}+T_{2}+T_{3}+T_{4}\Big{)}$
(62)

where $T_{j}=N_{j}/D$ with

$\displaystyle N_{1}$
$\displaystyle=\left[\frac{\mu_{e}}{3}-4(1-2x)~{}uS_{2}^{\prime}\right]^{2}$
(63) $\displaystyle N_{2}$
$\displaystyle=\left[\frac{\mu_{e}}{3}-4(1-2x)~{}uS_{2}^{\prime}\right]8S_{2}x_{e}\left[\frac{k_{F_{\mu}}}{k_{F_{e}}}-\frac{x_{\mu}}{x_{e}}\right]$
$\displaystyle N_{3}$
$\displaystyle=\left[\frac{k_{F_{\mu}}^{2}}{3\mu_{e}}-4(1-2x)~{}uS_{2}^{\prime}\right]\left(\mu_{e}+8S_{2}x_{e}\right)\left[\frac{x_{\mu}}{x_{e}}-\frac{k_{F_{\mu}}}{k_{F_{e}}}\right]$
$\displaystyle N_{4}$
$\displaystyle=\left[\frac{k_{F_{\mu}}^{2}}{3\mu_{e}}-4(1-2x)~{}uS_{2}^{\prime}\right]\frac{k_{F_{\mu}}}{k_{F_{e}}}\left[\frac{\mu_{e}}{3}-4(1-2x)~{}uS_{2}^{\prime}\right]$
$\displaystyle D$
$\displaystyle=\left[\frac{\mu_{e}}{3x_{e}}+8S_{2}\left(1+\frac{k_{F_{\mu}}}{k_{Fe}}\right)\right]$

and $k_{F_{e}}=\mu_{e}$ (massless electrons) and
$k_{F_{\mu}}$=$\sqrt{{\mu_{\mu}}^{2}-m_{\mu}^{2}}$. These equations explicitly
display the connection to the nuclear symmetry energy $S_{2}$ and its density
derivative $S_{2}^{\prime}$.

At the muon threshold ($k_{F_{\mu}},x_{\mu}$=0 $\Rightarrow$ $x$=$x_{e}$), it
is easy to see that $N_{2},N_{3},N_{4}$=0, while $N_{1}(\equiv N_{1}^{npe}$)
recovers Eq. (51) for $npe$ matter. At extremely high baryon density, muons
are ultra-relativistic ($\mu_{\mu}$=$k_{F_{\mu}}$=$k_{F_{e}}$,
$x_{\mu}$=$x_{e}$=$x/2$) so that $N_{2},N_{3}$=0,
$N_{1}$=$N_{4}$=$N_{1}^{npe}/2$ and the total leptonic contribution to the
sound speed is equally divided between electrons and muons.

####  IV.1.3  Sound speeds in Quark Matter

We now move to a discussion of sound speeds in dense quark matter at zero
temperature. For the pure quark phase, the difference of the two sound speeds
has been computed to leading order in the quark mass [66, 36] using the non-
interacting 3-flavor FG model with massive quarks (see Sec.III.2). These
expressions reveal that for the non-interacting FG model, a non-zero quark
mass is necessary to support $g$-modes. This is because a system of massless
$uds$ quarks is charge neutral with equal numbers of each flavor at any
density; effectively, there is no change in composition with density to drive
the $g$-mode.

To leading order in the $s$-quark’s mass $m_{s}$, the Brunt-Väisälä frequency
is [36]

$\displaystyle N_{q}\simeq\left(\frac{g}{2\pi
c_{e}}\right)\left(\frac{m_{s}^{2}}{\sqrt{B}}\right)\,,$ (64)

where $c_{e}^{2}=dp_{q}/d\epsilon_{q}$ is the equilibrium squared sound speed
in QM999Numerically, $N_{q}\approx 100$ Hz for a current quark mass
$m_{s}\approx 100$ MeV, but the effect of interactions in addition to this
yields significantly lower values for $N_{q}$ [36].. It is possible to obtain
an exact expression for $c_{e}^{2}$ and $c_{s}^{2}$ in QM for the FG model,
and also for the vMIT model, as we show below.

The equilibrium sound speed may be simply calculated by numerically evaluating
$c_{e,vMIT}^{2}=dp/d\epsilon$ in the pure quark phase. However, additional
insight into its compositional structure is gained by expressing it in terms
of the various chemical potentials involved. Starting from the relation (valid
in $\beta$-equilibrium)

$\displaystyle\mu_{n}=2\mu_{d}+\mu_{u}=(2\mu_{d}^{\ast}+\mu_{u}^{\ast})+3an_{Q}\,,$
(65)

where $\mu_{f}^{\ast}=\sqrt{k_{F_{f}}^{2}+m_{f}^{2}}$ for a quark of flavor
$f$, and using $c_{e}^{2}=d\ln\mu_{n}/\ln n_{B}$,101010In charge neutral and
$\beta$-equilibrated matter,
$\mu_{B}=\sum_{f}x_{f}\mu_{f}+\sum\limits_{\ell=e,\mu}x_{\ell}\mu_{\ell}=\mu_{n}$
as in the nucleonic phase. we obtain

$\displaystyle c_{e,vMIT}^{2}=\frac{1}{\mu_{n}}$ $\displaystyle\bigg{[}$
$\displaystyle\frac{1}{3}\left\\{2\mu_{d}^{\ast}\left(1-\frac{m_{d}^{2}}{\mu_{d}^{\ast^{2}}}\right)\frac{d\ln
n_{d}}{d\ln n_{B}}\right.$ (66) $\displaystyle+$
$\displaystyle\left.\mu_{u}^{\ast}\left(1-\frac{m_{u}^{2}}{\mu_{u}^{\ast^{2}}}\right)\frac{d\ln
n_{u}}{d\ln n_{B}}\right\\}$ $\displaystyle+$ $\displaystyle
3an_{Q}\bigg{]}\,.$

Contributions from the leptons are implicitly included in the above
expression.

For the non-interacting FG model, the pressure $p=\sum_{f}p_{\rm
FG}(\mu_{f},\mu_{e})$. Introducing the partial fractions $x_{f}=n_{f}/n_{B}$,
where $n_{f}=(\mu_{f}^{2}-m_{f}^{2})^{3/2}/\pi^{2}$ and
$n_{e}=\mu_{e}^{3}/3\pi^{2}$, the partial derivative of pressure with respect
to baryon density in the definition of the adiabatic sound speed in Eq. (39)
can be re-expressed in terms of partial derivatives with respect to the
various chemical potentials, yielding

$\displaystyle c_{s,{\rm FG}}^{2}=\frac{1}{\mu_{\rm
avg}}\left[\sum_{f}\frac{1}{3}\mu_{f}x_{f}\left(1-\frac{m_{f}^{2}}{\mu_{f}^{2}}\right)+\frac{1}{3}\mu_{e}x_{e}\right]\,,$
(67)

where $\mu_{\rm avg}=(\sum\limits_{f=u,d,s,e}n_{f}\mu_{f})/n_{B}$. Note that
if all $m_{f}$ = 0 (i.e, one is in a charge neutral phase with $x_{e}$ = 0),
$c_{s,{\rm FG}}^{2}=c_{e,{\rm FG}}^{2}=1/3$ and there can be no $g$-modes.
Inclusion of ${\cal O}(\alpha_{s})$ corrections to this model does not change
the fact that a non-zero quark mass is necessary for $g$-modes.

In the vMIT model given by Eqs. (24),
$\mu_{f}^{\ast}=\sqrt{k_{F_{f}}^{2}+m_{f}^{2}}=\mu_{f}-an_{Q}$, and as was
done for the FG model, we compute partial derivatives with respect to
$\mu_{f}^{\ast}$, noting that $n_{f}(\mu_{f})$ = $n_{\rm FG}(\mu_{f}^{\ast})$.
The resulting expression for the adiabatic sound speed in the vMIT model is

$\displaystyle c_{s,vMIT}^{2}=\frac{1}{\mu_{\rm avg}}$ $\displaystyle\bigg{[}$
$\displaystyle\sum_{f}\frac{1}{3}\mu_{f}^{\ast}x_{f}\left(1-\frac{m_{f}^{2}}{\mu_{f}^{\ast^{2}}}\right)$
(68) $\displaystyle+$
$\displaystyle\frac{1}{3}\mu_{e}x_{e}+3an_{Q}\bigg{]}\,,$

where all quantities $\mu_{f},\mu_{e},x_{e},x_{f}$ are equilibrium
values111111Inclusion of muons is straightforward and adds a term,
$\frac{1}{3}\mu_{\mu}x_{\mu}\left(1-\frac{m_{\mu}^{2}}{\mu_{\mu}^{\ast^{2}}}\right)$,
on the right hand side of Eq. (67).. If we switch off interactions
($a\rightarrow$ 0), we recover results of the non-interacting FG model.
Interestingly, if we retain the interaction term, but set all quark masses
equal or to zero (implying that $x_{e}$ = 0), we find that
$c_{s,vMIT}^{2}=c_{e,vMIT}^{2}$ so stable $g$-modes are not supported in the
pure quark phase. Therefore, while both sound speeds are modified by
interactions, e.g.,

$\displaystyle c_{s,vMIT}^{2}=\frac{1}{\mu_{\rm
avg}}\left[\mu_{q}^{*}+3an_{Q}\right]\neq c_{s,{\rm FG}}^{2}\,,$ (69)

at asymptotically high density where quark masses are negligible, there can be
no $g$-modes in quark matter in the vMIT model 121212$g$-modes would still
exist in a mixed phase of vMIT quark matter and nucleons as the electron
fraction would vary from $\beta$-processes involving nucleons.

Note that when all chemical potentials and partial fractions are set to their
equilibrium values for $c_{s,vMIT}^{2}$ in the pure phase, $\mu_{\rm
avg}=\mu_{n}$. A comparison of Eqs. (66) and (68) reveals the differences
between the two sound speeds. While effects of interactions enter in the same
formal way for the two squared speeds, the occurrence of the logarithmic
derivatives of the quark densities distinguishes $c_{e,vMIT}^{2}$ from
$c_{s,vMIT}^{2}$ which features the partial fractions $x_{f}$. This difference
is the principal reason for the latter to become larger than the former. In
both cases, the $d$-quark contributions are larger than those of $u$ and $s$
quarks.

###  IV.2 Sound Speeds in the Mixed Phase

Once we have expressions for the sound speed in a pure phase of quarks or
nucleons, it is possible to compute the sound speed in the mixed phase of the
two, obtained from a Gibbs construction. The only information required, other
than the sound speeds in the pure phases, is the partial phase fraction of
quarks $\chi$ at any density. It is more convenient to begin with the
reciprocal relation

$\displaystyle\frac{1}{c_{e,{\rm mix}}^{2}}=\left(\frac{d\epsilon_{\rm
mix}}{dp_{\rm mix}}\right)\,.$ (70)

In a Gibbs mixed phase, the pressures in the two phases are equal, while the
energy density is a proportional mix of the quark ($q$) and nucleonic/hadronic
($h$) phases: $\epsilon_{\rm mix}=(1-\chi)\epsilon_{h}+\chi\epsilon_{q}$.
Substituting this in Eq.(70) gives

$\displaystyle\frac{1}{c_{e,{\rm
mix}}^{2}}=\frac{(1-\chi)}{c_{e,h}^{2}}+\frac{\chi}{c_{e,q}^{2}}+(\epsilon_{q}-\epsilon_{h})\frac{d\chi/dn_{B}}{dP/dn_{B}}\,.$
(71)

The derivatives in Eq.(71) must be computed numerically after solving for
$\chi$, hence afford no particular advantage over a direct numerical
computation of the sound speed from Eq.(70) itself. However, note that the
last term in Eq.(71) is always positive in the mixed phase.

As before, the general definition of the adiabatic sound speed applies to the
mixed phase

$\displaystyle c_{s,{\rm mix}}^{2}=\left(\frac{dp_{\rm mix}}{d\epsilon_{\rm
mix}}\right)_{x_{i}={\rm const.}=x_{i,{\rm eq.}}}\,,$ (72)

and the thermodynamic identity becomes $\epsilon_{\rm mix}+p_{\rm mix}$ =
$\sum_{i}n_{i}\mu_{i}$ = $n_{B}\mu_{\rm avg}$. Noting that the derivatives
$\partial\epsilon_{h}/\partial_{n_{B,h}}$ and
$\partial\epsilon_{q}/\partial_{n_{B,q}}$ are equal to the respective
$\mu_{\rm avg}$, it is once again more convenient to begin with the reciprocal
relation

$\displaystyle\frac{1}{c_{s,{\rm mix}}^{2}}=\left(\frac{d\epsilon_{\rm
mix}}{dp_{\rm mix}}\right)_{x_{i}={\rm const.}=x_{i,{\rm eq.}}}\,,$ (73)

and use the chain rule to compute derivatives with respect to density. This
leads to

$\displaystyle\left(\frac{d\epsilon_{\rm mix}}{dp_{\rm
mix}}\right)_{x_{i}=x_{i,{\rm eq.}}}=\frac{(1-\chi)\mu_{\rm
avg}}{\left(\frac{\partial p_{h}}{\partial n_{B},h}\right)}+\frac{\chi\mu_{\rm
avg}}{\left(\frac{\partial p_{q}}{\partial n_{B},q}\right)}\,,$ (74)

which, using Eq. (39), becomes

$\displaystyle\frac{1}{c_{s,{\rm
mix}}^{2}}=\frac{(1-\chi)}{c_{s,h}^{2}}+\frac{\chi}{c_{s,q}^{2}}\,.$ (75)

Comparing Eqs. (71) and (75), and to the extent that the two sound speeds in
the pure hadronic/quark phase are almost equal, we expect that the last term
in Eq. (71) which tracks the rapidly changing composition in the mixed phase,
is mainly responsible for $c_{s,\rm mix}^{2}>c_{e,\rm mix}^{2}$. The more
rapid the appearance of new chemical species and the softer the mixed phase,
the larger the Brunt-Väisälä frequency will be.

Furthermore, as will become evident from our results in Sec. V, the adiabatic
sound speed is continuous across the transition to and from the mixed phase,
while the equilibrium sound speed has a slight jump to accommodate the
derivative of $\chi$. The reciprocal relation for the adiabatic sound speeds
is reminiscent of the addition of resistors in a parallel circuit, with
voltage as pressure and current as energy density. Such impedance analogies
arise commonly in electrical engineering when modeling the behavior of
transducers.

##  V Results

###  V.1 Structural properties, sound speeds and the Brunt-Väisälä frequency

Figure 1:  Mass-radius curves for the ZL EOS without and with muons.
Configurations with muons are slightly more compact, but both cases support
$M_{\rm max}\simeq 2M_{\odot}$. Except for $L$, the EOS parameters are
$K_{0}=220$ MeV, $S_{v}=31$ MeV, and $\gamma_{1}=1.6$ for all curves.

Figure 1 shows $M$-$R$ curves for ZL EOSs with and without muons for the
indicated parameters in the caption. The radii of $\sim 1.4~{}M_{\odot}$
stars, $R_{1.4}$, for the different models shown lie within the bounds
inferred from available data. For example, data from X-ray observations have
yielded $R_{1.4}=9$-$14$ km for canonical masses of $\sim 1.4~{}M_{\odot}$
[67, 68, 69]. Measured tidal deformations from gravitational wave data in the
binary NS merger GW170817 give 8.9-13.2 km for binary masses in the range
1.36(1.17)-1.6(1.36) $M_{\odot}$ [70], whereas for the same masses Capano et
al. [71] report $11\pm 1$ km. X-ray pulse analysis of NICER data from PSR
J0030+0451 by Miller et al. (2019) [39] finds $13.02^{+1.14}_{-1.19}$ km for
$M=1.44\pm 0.15~{}M_{\odot}$, whereas for the same star Riley et al. (2019)
[38] obtain $12.71^{+1.14}_{-1.19}$ km and
$M=1.34^{+0.15}_{-0.16}~{}M_{\odot}$. The maximum masses ($\simeq 2M_{\odot}$)
of these EOSs131313By adjusting the constants of the ZL EOS to make the EOS
stiffer (yet causal) at supra-nuclear densities, masses larger than
$2M_{\odot}$ can be obtained; an example will be shown later. are also within
the uncertainties of high mass NSs which range from $1.908\pm
0.016~{}M_{\odot}$ to $2.27^{+0.17}_{-0.15}~{}M_{\odot}$ [72, 73, 74, 75, 76].
Although differences in $R_{1.4}$ with and without muons for a given EOS are
small, the appearance of muons in the star leads to distinct features in the
Brunt-Väisälä frequency (see below).

Figure 2:  Difference between the adiabatic and equilibrium squared sound
speeds (normalized to the squared speed of light) for the ZL EOS ($K_{0}$=220
MeV, $S_{v}=31$ MeV, $L=60$ MeV and $\gamma_{1}$=1.6) without and with muons.

In Fig. 2, differences in the two squared sound speeds are shown as a function
of $n_{B}$ with and without muons for the ZL EOS with $L=60$ MeV. The small
jump at $n_{B}\simeq 0.14~{}{\rm fm}^{-3}$, the density at which muons appear,
is caused by a sudden drop in the equilibrium sound speed. The differences at
large densities are due to the increasing concentration of muons.

Figure 3:  The Brunt-Väisälä frequency in the NS for the ZL EOS ($K_{0}$=220
MeV, $S_{v}=31$ MeV, $L=60$ MeV, and $\gamma_{1}$=1.6) without and with muons.
Figure 4:  EOS for the mixed phase of nucleons and quarks (middle curve) using
the Gibbs construction. For the ZL EOS without muons, $K_{0}$=220 MeV,
$S_{v}=31$ MeV, $L=60$ MeV, and $\gamma_{1}$=1.6. Parameters for the vMIT EOS
are: ($m_{u},m_{d},m_{s}$)=(5,7,150) MeV, $B^{1/4}$=180 MeV and $a$=0.1. The
circle indicates the central $p$ and $\epsilon$ of the maximum mass star
($n_{c,{\rm max}}/n_{s}$=7.63, for $M_{\rm max}/M_{\odot}$=1.82).

Figure 3 shows the Brunt-Väisälä frequency $N$ vs $r/R$ in the star. In the
results shown, the crust is assumed to be a homogeneous fluid for simplicity,
hence $N$ vanishes there. The location where muons begin to appear is signaled
by the small kink in the bottom panel. Overall, $N$ is slightly larger with
muons in the density range in the core of a 1.4$M_{\odot}$ star, consistent
with Fig. 2. This has a proportional impact on the $g$-mode frequency as shown
in the next section.

The EOS of the mixed phase following the Gibbs construction, and the ZL EOS
for the nucleonic sector and the vMIT EOS for the quark sector is shown in
Fig. 4. The ZL EOS does not include the small effect of muons. In the quark
sector, muons have not been included since their impact relative to quarks is
tiny.

Figure 5: Quark fraction vs $n_{B}$ corresponding to Fig. 4. The circle
indicates the central density of the maximum mass star ($n_{c,{\rm max}}$=1.22
${\rm fm}^{-3}$ for $M_{\rm max}/M_{\odot}$=1.82).

The compositional change in the mixed phase is indicated by the quark fraction
$\chi$ in Fig. 5. The steep rise of $\chi$ from the onset indicates the sort
of rapid compositional change that can impact the $g$-mode frequency. A
similar effect has been reported [30] in the context of the appearance of
strange baryons (e.g. hyperons), which is not a phase transition. Note, that
for the EOSs considered, the central density of the maximum mass star,
indicated by the filled circle on the curve, lies in the mixed phase so that
the pure quark phase is not realized.

Figure 6:  The two sound speeds (top panel) and their differences (bottom
panel) in the mixed phase for the EOS parameters corresponding to Fig. 4. The
pure quark phase is not achieved prior to the maximum mass in this case. The
termination at $n_{B}$=0.08 fm-3 demarcates the core-crust boundary, since we
assume $c_{s}$=$c_{e}$ in the core. Both sound speeds take much smaller values
in the crust than in the core. The circle indicates the central density of the
maximum mass star ($n_{c,{\rm max}}$=1.22 ${\rm fm}^{-3}$ for $M_{\rm
max}/M_{\odot}$=1.82).

Figure 6 shows results for the individual sound speeds and their differences
for the mixed phase. The two sound speeds in the mixed phase behave very
differently. Specifically, the equilibrium sound speed suddenly drops (rises)
at the onset (end) of the mixed phase, whereas the adiabatic sound speed
varies smoothly.

Figure 7: The Brunt-Väisälä frequency in a hybrid star of mass
$1.4~{}M_{\odot}$. The ZL EOS does not include muons and parameters for the
nuclear and quark EOS are as in Fig. 4. Quarks enter at $n_{B}\simeq$ 0.42
fm-3 corresponding to $r/R$ = 0.383, and the mixed phase extends beyond the
central density. The value of $N$ decreases towards the core due to the
decreasing value of $g$, even as the sound speed difference does not change
much.

The Brunt-Väisälä frequency of a $1.4~{}M_{\odot}$ hybrid star is shown in
Fig. 7. Note the broader width of the peak when quarks enter, and its location
in denser regions of the star, as compared to the nucleonic stars depicted in
Fig. 3. This explains why the $g$-mode, which is a long-wavelength global
oscillation, is strongly impacted by the mixed phase (see results in the next
section).

Figure 8:  The mass-radius curves for a hybrid star (Gibbs construction) with
EOS parameters chosen such that the mixed phase supports $M_{\rm
max}=2.05~{}M_{\odot}$. In the left panel, muons are not included, whereas the
right panel is with muons included.

Figure 8 shows $M$-$R$ curves for a hybrid star whose $M_{\rm
max}=2.05~{}M_{\odot}$. This value is obtained by increasing the
compressibility parameter of the ZL EOS from $K_{0}$ = 220 to $K_{0}$ = 240
MeV, and increasing $\gamma_{1}$ from 1.6 to 1.8 while maintaining causality.
Including muons pushes the onset of the mixed phase to slightly higher
densities, which causes the maximum mass of a hybrid star with muons to be
higher than for a hybrid star without muons. This is in contrast to the effect
of muons in an ordinary NS, where the softening results in a lower maximum
mass. The leftmost curves in these figures refer to a self-bound quark star,
and are shown here to provide contrast.

###  V.2 Boundary conditions for the $g$-mode oscillation

Having established the equilibrium structure and computed the sound speeds, we
have all the variables necessary to solve Eqs. (1) at hand, except for the
boundary conditions that determine the (real) eigenfrequencies. The boundary
conditions for Newtonian structure equations are obtained as a straightforward
limiting case of Eqs. (1), and are discussed at length in [20]. To summarize
those results, in the Newtonian non-relativistic case, regularity of
$\xi_{r},~{}\delta p/\rho$ can be checked by Taylor expansion around $r$ = 0.
The resulting condition is:

$\displaystyle
r^{2}\xi_{r}=\frac{l}{\omega^{2}}(Y_{0}+\phi_{0})r^{l+1},\quad\frac{\delta
p}{\rho}=Y_{0}r^{l}\,,$ (76)

where $Y_{0},~{}\phi_{0}$ are constants. For our purposes, $\phi_{0}$ = 0
since we ignore perturbations in the gravitational potential, as in [20].
$Y_{0}$ is an arbitrary normalization constant allowed by the linearity of
these equations. Effectively, this means that the overall scale of the
eigenfunctions is arbitrary. It must be determined by external factors, such
as the strength of the force (tidal effects in a merger, for example). The
normalization has no impact on the numerical value of the eigenfrequency. It
is therefore conventional to choose $Y_{0}$ = 1. We will make, for simplicity,
and without loss of generality, a slightly different choice:

$\frac{l}{\omega^{2}}Y_{0}=1$ (77)

so that (for $l$=2), $\xi_{r}\rightarrow r$ as $r\rightarrow 0$. In practice,
we start the integration slight off-center, so $\xi_{r}$ will be small but
non-zero. The other condition at the center becomes

$\frac{\delta p}{\rho}=\frac{\omega^{2}}{l}r^{l}{\rm e}^{-\nu_{0}}\,,$ (78)

again, with $l$ = 2 for our case. For the relativistic form of the oscillation
equations, the above conditions still apply with the change $\frac{\delta
p}{\rho}\rightarrow\frac{\delta p}{\epsilon+p}$. The boundary condition at the
surface is the vanishing of the Lagrangian pressure perturbation $\Delta
p=c_{s}^{2}\Delta\epsilon=0$. This projects out the radial component of
$\vec{\xi}$. In the non-relativistic case, $\nabla p=-\rho g$ while in the
relativistic case, $\nabla p=-gh$ with $h=(\epsilon+p)$ the enthalpy.

With some algebra, one can arrive at a simpler form of Eqs. (1):

$\displaystyle\frac{dU}{dr}$ $\displaystyle=$
$\displaystyle\frac{g}{c_{s}^{2}}U+{\rm e}^{\lambda/2}\left[\frac{l(l+1){\rm
e}^{\nu}}{\omega^{2}}-\frac{r^{2}}{c_{s}^{2}}\right]V$
$\displaystyle\frac{dV}{dr}$ $\displaystyle=$ $\displaystyle{\rm
e}^{\lambda/2-\nu}\frac{\omega^{2}-N^{2}}{r^{2}}U+g\Delta(c^{-2})V\,,$ (79)

where $U$ = $r^{2}{\rm e}^{\lambda/2}\xi_{r}$, $V$ = $\delta p/(\epsilon+p)$
and
$\Delta(c^{-2})=c_{e}^{-2}-c_{s}^{-2}$.

We employ a 4th-order Runge-Kutta scheme to find a global solution of the
linear perturbation equations, Eqs. (V.2), subject to the boundary conditions
for the relativistic case outlined above. Since the solution set comprises
overtones, we selected the lowest order $g$-mode (highest frequency) by
checking that the radial eigenfunction $\xi_{r}$ has only one node inside the
star. The corresponding eignefrequency is plotted in the figures that follow.
We examine the trends of the $g$-mode vs. mass for various parameter choices,
for the pure nuclear, self-bound and hybrid stars.

Figure 9: Contrasts of the $g$-mode frequencies vs mass of normal NSs for the
ZL EOS without and with muons. The two curves with different $L$’s in each
panel are for EOSs with $K_{0}=220$ MeV, $S_{v}=31$ MeV and $\gamma_{1}$=1.6.

Figure 9 contrasts the influence of varying the density dependence of the
symmetry energy, by changing the slope of the symmetry energy parameter $L$ at
$n_{s}$, of the underlying ZL EOS for normal neutron stars with fixed
$K_{0}=220$ MeV and $S_{v}=31$ MeV. For $L$ = 60 MeV as well as $L$=70 MeV,
the softening effect of muons leads to a noticeable increase in the $g$-mode
frequency at a given mass. Comparing $L$ = 60 MeV with $L$ = 70 MeV for a
fixed composition however, the $g$-mode frequency for
$M\buildrel>\over{\sim}~{}$0.5-0.6 $\mathrm{M}_{\odot}\;$is higher for the
stiffer EOS.

Figure 10: Contrasts of g-mode frequencies vs stellar mass in a hybrid star.
Parameters of the EOSs are as in the insets. In the left panel, muons are not
included, whereas the right panel is with muons included.

In Fig. 10, results contrasting the $g$-mode frequencies in normal, hybrid,
and self-bound stars are presented. The contents of this figure constitute the
principal result of this work, viz., the abrupt rise in the scale of the
$g$-mode frequency at the onset of the mixed phase in the hybrid star. For the
EOS parameters displayed in the figure, the jump occurs around 1.4
$M_{\odot}$, so that a hybrid star in a merger would have a distinctly higher
g-mode frequency than a normal NS. In the top panel, the ZL EOS does not
include muons, whereas in the bottom panel the ZL EOS includes muons. The
$g$-mode frequency in the mixed phase is again higher than in a pure phase,
but since the mixed phase appears at a higher density due to muons, the rise
in the $g$-mode is less dramatic compared to a hybrid star without muons.
Results for the self-bound star are shown here for comparison, and to
emphasize that its $g$-mode frequency is comparatively small (10-50 Hz).
Unlike the $f$-mode frequency for the hybrid star, which gradually
interpolates between those of the normal NS and self-bound star [50, 77] and
shows no dramatic effects of compositional changes, the $g$-mode frequency for
the hybrid star is the highest of all and is sensitive to the onset of quarks
- making it less subject to ambiguity. One does not need to know the mass of
the star to ascertain if it can be a hybrid star if the $g$-mode frequency can
be precisely determined.

The unusually large $g$-mode frequency for the hybrid star with a Gibbs mixed
phase may be understood in a qualitative sense using general thermodynamic
considerations without reference to details of the EOS. In general, the
equilibrium sound speed $c_{e,\mathrm{mix}}$ in a system with two conserved
charges ($\mu_{B}$ and $\mu_{Q}$) can be expressed as

$\displaystyle c_{e,\mathrm{mix}}^{2}$ $\displaystyle=$
$\displaystyle\frac{dp_{\mathrm{mix}}\left(\mu_{B},\mu_{Q}\right)}{d\epsilon_{\mathrm{mix}}}$
(80) $\displaystyle=$ $\displaystyle\frac{\partial
p_{\mathrm{mix}}}{\partial\mu_{B}}\left(\frac{d\mu_{B}}{d\epsilon_{\mathrm{mix}}}\right)+\frac{\partial
p_{\mathrm{mix}}}{\partial\mu_{Q}}\left(\frac{d\mu_{Q}}{d\epsilon_{\mathrm{mix}}}\right)$

where $\mu_{Q}$ is the charge chemical potential. Glendenning [61] showed that
in such a situation, while $\mu_{B}$ is smooth at the onset of the mixed
phase, $\mu_{Q}$ is not, as there is freedom to rearrange charges between the
two phases to achieve global charge neutrality and minimize the free energy.
In fact, the steady rise with density of $\mu_{Q}$ in the pure nuclear phase
changes abruptly to a decline in the mixed phase, tempering the equilibrium
sound speed as shown by our numerical results presented in Fig. 6 and
confirmed by other works [78] which use different EOS from ours for the
nucleon and quark sector. On the other hand, the adiabatic sound speed
$c_{s,\mathrm{mix}}$ is evaluated at fixed composition and shows no such
effect, hence the difference of the two sound speeds (usually small in a pure
phase) abruptly increases in the mixed phase. This is reflected as a positive
jump in the Brunt-Väisälä frequency and therefore of the $g$-mode in the mixed
phase.

##  VI $g$-mode Energy and Tidal Forcing

Unlike the Sun, where convection from within can drive oscillations, any
oscillations of an evolved NS likely require an external agent to excite and
sustain the perturbation beyond its normal damping time. A violent event such
as a NS merger is bound to produce oscillations in the pre-merger phase due to
tidal forcing or in the postmerger (ringdown) phase as the hypermassive
remnant relaxes to its stable rotating configuration. Here, we estimate the
impact of the $g$-mode on tidal phasing leading up to the merger, as the
$g$-mode spectrum in the postmerger remnant can be modified by thermal and
convective effects which are beyond the scope of the current work. We follow
[18] and assume spherically symmetric non-rotating stars, the Newtonian
approximation to orbital dynamics and quadrupolar gravitational wave emission.
These simplifying approximations allow for a first estimate of the excitation
energy and amplitude of the $g$-mode, as well as the phase difference due to
dynamic tides associated to the $g$-mode (not to be confused with the quasi-
static tides due to global deformation). Our estimates can be systematically
improved by going to the post-Newtonian approximation or numerical relativity.

The estimates are derived by modeling the NS as a forced simple harmonic
oscillator with a mass $M_{\ast}$, radius $R_{\ast}$ and a natural frequency
$\omega$=$\omega_{g}$, the angular frequency of the $g$-mode. The forcing
comes from the quadrupolar moment of the companion star’s gravitational force
(mass $M$), which couples to the $g$-mode. By following the analysis of [18],
we arrive at an expression for the accumulated phase error $\Delta\Phi(t)$
caused by the $g$-mode:

$\displaystyle\Delta\Phi(t)\approx\frac{3\pi\Gamma}{m}\left[\frac{\Omega_{e}(t)}{\Omega}-1\right]\left(-\frac{\Delta
E}{E}\right)\,,$ (81)

where $\triangle E$ is the energy pumped into the $g$-mode, $E$ the total
(kinetic plus potential) orbital energy of the system, $\Omega_{e}(t)$ the
time-dependent orbital frequency of the binary, and $\Omega$=$\omega_{g}/m$.
The quantity $m$ in Eq. (81) is the azimuthal mode index ($m$=2 in this case).
Finally, $\Gamma$ is a quantity that appears as a result of applying the
stationary phase approximation to the evaluation of the time to resonance
[18], and is quantified below. A $\Delta\Phi(t)$ of ${\cal O}(1)$ signifies a
large deviation from the point particle approximation to the gravitational
waveform from the merger. Explicitly, the quantity $\triangle E$ (for angular
quantum number $l$) is given by

$\displaystyle\Delta E$ $\displaystyle=$
$\displaystyle\left(\frac{5\pi}{384m}\right)\frac{M/M_{\ast}}{[1+M/M_{\ast}]^{(2l+1)/3}}\left(\frac{c^{2}R_{\ast}}{GM_{\ast}}\right)^{5/2}$
(82) $\displaystyle\times$
$\displaystyle\left(\frac{\Omega}{\Omega_{d}}\right)^{(4l-7)/3}\left(\frac{GM_{\ast}^{2}}{R_{\ast}}\right)S_{lm}^{2}\,.$

where $\Omega_{d}$=$(GM_{\ast}/R_{\ast}^{3})^{1/2}$ is a natural frequency
unit and $S_{lm}$ is proportional to the overlap integral between the mode
eigenstate $|lm\rangle$ and the vector spherical harmonic
$\left|P_{lm}\right\rangle=\nabla\left[r^{l}Y_{lm}(\theta,\phi)\right]$. The
total instantaneous orbital energy is $E=-GMM_{\ast}/2a$ with $a=a(t)$ the
instantaneous orbital separation. The evolution of the orbital frequency for a
circularized orbit using the formula for quadrupolar gravitational wave
emission gives

$\displaystyle\Omega_{e}(\tau)=\frac{1}{8}\left(\frac{aM_{c}}{c^{3}}\right)^{-5/8}\frac{1}{\tau^{3/8}}\,,$
(83)

where $M_{c}$ is the chirp mass of the binary system and $\tau$ is the time to
coalescence. All quantities on the right hand side in Eq. (81) can be
calculated, once the parameters of the binary ($M,M_{\ast},R_{\ast}$) and the
resonant $g$-mode frequency are fixed. We choose $M$=$M_{\ast}$=1.5
$M_{\odot}$ for neutron/hybrid stars and pure quark stars. The strongest tidal
coupling is likely to the $l$=$m$=2 $g$-mode whose characteristic frequency we
choose as

$\displaystyle\omega_{g}(NS)$ $\displaystyle\cong$ $\displaystyle
2\pi(200)\,{\rm Hz}\,,~{}~{}$ $\displaystyle\omega_{g}(HS)$
$\displaystyle\cong$ $\displaystyle 2\pi(300)\,{\rm
Hz}\,,~{}~{}\rm{and}~{}~{}$ $\displaystyle\omega_{g}(QS)$ $\displaystyle\cong$
$\displaystyle 2\pi(40)\,{\rm Hz}\,,~{}~{}$ (84)

based on the $g$-mode eigenfrequencies in the previous section. Even without
computing $S_{lm}$, one can estimate from Eq. (83) the time at which the
$g$-mode becomes resonant as

$\displaystyle\tau_{0}(NS)$ $\displaystyle\cong$ $\displaystyle 272\,{\rm
ms}\,,~{}~{}$ $\displaystyle\tau_{0}(HS)$ $\displaystyle\cong$ $\displaystyle
103\,{\rm ms}\,,~{}~{}\rm{and}~{}~{}$ $\displaystyle\tau_{0}(QS)$
$\displaystyle\cong$ $\displaystyle 22\,{\rm s}\,,$ (85)

where the the zero of time is the moment of coalescence. Assuming circularized
orbits, standard equations of binary orbit evolution $a(t)$ [79] give

$\displaystyle a_{0}(NS)$ $\displaystyle\cong$ $\displaystyle
111~{}\mathrm{km}\,,~{}~{}$ $\displaystyle a_{0}(HS)$ $\displaystyle\cong$
$\displaystyle 85~{}\mathrm{km}\,,~{}~{}{\rm and}~{}~{}$ $\displaystyle
a_{0}(QS)$ $\displaystyle\cong$ $\displaystyle 326~{}\mathrm{km}\,.$ (86)

We note that the $g$-mode for the hybrid star, which has a larger resonant
frequency than neutron or quark stars, is excited later in the merger and is
likely to be stronger in amplitude due to the close separation of the binary
since the forcing term is $\propto 1/a^{3}$ for $l$=2. Finally, from our
calculations for the $g$-mode eigenfunction and the associated density
perturbation $\delta\epsilon(r)$, we estimate

$\displaystyle S_{lm}^{NS}$ $\displaystyle\cong$ $\displaystyle\quad 4.5\times
10^{-3}\,,~{}~{}$ $\displaystyle S_{{lm}}^{HS}$ $\displaystyle\cong$
$\displaystyle\quad 6.2\times 10^{-3}\,~{}~{}{\rm and}~{}~{}$ $\displaystyle
S_{lm}^{QS}$ $\displaystyle\cong$ $\displaystyle\quad 9.9\times 10^{-6}\,$
(87)

using Eq. (40) for $S_{lm}$ from [32]. From these estimates, Eq. (82) can be
utilized to yield the estimated fractional orbital energy pumped into the
$g$-mode:

$\displaystyle\left|\frac{\Delta E}{E}\right|^{NS}$ $\displaystyle\cong$
$\displaystyle 2.3\times 10^{-3}\,,~{}~{}$ $\displaystyle\left|\frac{\Delta
E}{E}\right|^{HS}$ $\displaystyle\cong$ $\displaystyle 5.9\times
10^{-3}\,,~{}~{}~{}~{}{\rm and}~{}~{}$ $\displaystyle\left|\frac{\Delta
E}{E}\right|^{QS}$ $\displaystyle\cong$ $\displaystyle\quad 2\times
10^{-9}\,,$ (88)

and finally from Eq. (81), we obtain the phase error due to the resonant
excitation of the $g$-mode to be

$\displaystyle\triangle\phi^{NS}$ $\displaystyle\cong$ $\displaystyle
0.8\,,~{}~{}{\rm and}~{}~{}$ $\displaystyle\triangle\phi^{HS}$
$\displaystyle\cong$ $\displaystyle 0.45\,,~{}~{}$
$\displaystyle\triangle\phi^{QS}$ $\displaystyle\cong$ $\displaystyle 6\times
10^{-4}\,.$ (89)

Note that $\Delta\phi^{NS}$ and $\triangle\phi^{HS}$ are comparable. Despite
$\left(\frac{\Delta E}{E}\right)$ being larger for a hybrid star as expected,
its higher $g$-mode frequency means it is excited later in the merger, when
there is less time left for accumulating a phase error. These results are very
sensitive to the value of $S_{lm}$ ($\Delta\Phi\propto S_{lm}^{2}$), which
itself can vary by a factor of 2 or more depending on the EOS.

### Comparison with other works

Table 1: Comparison of characteristic $g$-mode frequencies (denoted by $\omega_{g}$ in the table) reported in a selection of the literature. As other works usually fix the stellar mass $M$, we include this information. The symbol $\Lambda$ is used here as a shorthand to denote hyperonic degrees of freedom and SF denotes superfluidity in the nucleonic sector. Values of $f_{g}$ that vary with the NS mass can be inferred from Figs. 9 and 10 of this work. The entries are representative, not exhaustive. Authors [Ref.] | Core | $M$ | $f_{g}=\omega_{g}/(2\pi)$
---|---|---|---
| Composition | [$M_{\odot}$] | [kHz]
Reisenegger & Goldreich [17] | $npe$ | 1.405 | 0.215
Lai [20] | $npe$ | 1.4 | 0.073
Kantor & Gusakov [29] | $npe$ | 1.4 | 0.13
Kantor & Gusakov [29] | $npe\mu$ | 1.4 | 0.19
Kantor & Gusakov [29] |  $npe\mu$(SF) | 1.4 | 0.46
Dommes & Gusakov [30] |  $npe\mu\Lambda$(SF) | 1.634 | 0.742
Yu & Weinberg [33] | $npe\mu$ | 1.4 | 0.13
Yu & Weinberg [33] |  $npe\mu$(SF) | 2.0 | 0.45
Rau & Wasserman [34] |  $npe\mu$(SF) | 2.0 | 0.45
Jaikumar et al. [this work] | $npe$ | 1.4 | 0.24
Jaikumar et al. [this work] | $npe\mu$ | 1.4 | 0.27
Jaikumar et al. [this work] | $npe\mu q$ | 2.0 | 0.58

Table 1 compares our results for zero-temperature core $g$-modes in the Gibbs
mixed phase of hadrons and quarks with other works, some of which also find an
enhancement of the frequency due to other compositional mixes or collective
fluid effects like superfluidity, although values in the table do not include
the effect of entrainment on the $g$-mode, which has also been studied.
Details about the different EOSs used, the effect of non-zero temperature and
entrainment can be found by perusing the respective reference. We confirm the
value of the $g$-mode frequency for $npe$ and $npe\mu$ non-superfluid matter
described by the Akmal-Pandharipande-Ravenhall (APR)-EOS as reported in [29],
which also serves as a test of our numerics. In comparison to [29] with the
APR-EOS or [33] with the SLy4 equation of state, the ZL-EOS yields a larger
value of the $g$-mode frequency as it is less stiff than either of those two.
While the EOS and the treatment of gravitational perturbations differ between
the cited works, the results for $npe\mu$ matter with superfluidity appear to
be in general agreement with each other. Note the considerably larger value of
the $g$-mode frequency for hyperonic stars with superfluidity compared to
hybrid stars or superfluid neutron stars. All of these, in turn, are larger
than non-superfluid neutron/hyperonic stars although the latter employ
Newtonian gravity. A study of $g$-mode frequencies and damping times in
superfluid hybrid stars is a future objective that would make this comparison
more complete.

##  VII Summary and Conclusions

The main objective of this work was to ascertain the characteristics of
$g$-mode oscillations of NSs containing QM in their cores. Toward this end,
the nucleon-to-quark phase transition was treated using Gibbs construction
which renders an otherwise sharp first order transition smooth. The cores of
such hybrid stars accommodate admixtures of nucleonic and quark matter, the
pure quark phase being never achieved. This feature, while allowing contrasts
between the structural properties (e.g., $M$ vs $R$) of normal and hybrid
stars to be made also permits comparisons of observables that depend on their
interior compositions, such as short- and long-term cooling, oscillation
modes, etc. Determining the composition of the star is essential to break the
degeneracy that exists in the masses and radii of normal and hybrid stars as
one may be masquerading as the other.

The nucleonic part of the EOS used in this work tallies with nuclear
systematics near and below $n_{s}$ in addition to being consistent with
results from modern chiral EFT calculations up to $2n_{s}$ for which
uncertainty quantifications have been made. The EOS employed in the quark
sector is sufficiently stiff, hence non-perturbative, to support $\sim
2M_{\odot}$ NSs required by recent observations. Furthermore, the overall EOS
gives radii of $\sim 1.4M_{\odot}$ stars that lie within the bounds of recent
determinations. The EOS is also consistent with the tidal deformation inferred
from gravitational wave detection in the event GW170817. Appendix A summarizes
the structural properties for the EOSs used and provides mathematical details
for the derivation of the sound speeds.

Unlike for $M$-$R$ curves for which only the pressure vs density relation
(EOS) is sufficient, the analysis of $g$-mode oscillations requires
simultaneous information about the equilibrium and adiabatic squared sound
speeds, $c_{e}^{2}=dp/d\epsilon$ and $c_{s}^{2}=\partial
p/\partial\epsilon|_{x}$, where $x$ is the local proton fraction. The
distinction between these two sound speeds plays a central role in determining
the Brunt-Väisälä frequencies $\omega^{2}\propto c_{e}^{-2}-c_{s}^{-2}$ of
non-radial $g$-mode oscillations. Thus, a future detection of $g$-modes would
take gravitational wave astronomy beyond the current capability of $M$-$R$
measurements to determine the composition of the star.

We find that the $g$-mode is sensitive to the presence of QM in NSs, where
quarks are part of a mixed phase with nuclear matter in the core. The
equilibrium sound speed drops sharply at the boundary of the mixed phase (Fig.
5), raising the local Brunt-Väisälä frequency and the fundamental $g$-mode
frequency of the star (Fig. 6). Contrasts of $g$-mode frequencies between
normal and hybrid stars containing quark matter (Fig. 9) form the principal
results of our work.

Our analysis leads to the conclusion that in binary mergers where one or both
components may be a hybrid star, the fraction of tidal energy pumped into the
resonant $g$-mode in hybrid stars can exceed that of a NS by a factor of 2-3,
although resonance occurs later in the inspiral. On the other hand, a self-
bound star has a much weaker tidal overlap with the g-mode. The cumulative
tidal phase error in hybrid stars $\Delta\phi\cong$ 0.5 is comparable to that
from tides in ordinary NSs. While this happenstance may present a challenge in
distinguishing between the two cases, should the g-mode be excited to
sufficient amplitude in a postmerger remnant, its frequency spectrum would be
a possible indication for the existence of non-nucleonic matter, including
quarks. The detection of such $g$-mode frequencies in binary mergers observed
by current gravitational wave detectors seems challenging, but possible with
next generation detectors.

The novel features of this work include (i) use of nucleonic EOSs that are
consistent with constraints from modern chiral EFT calculations coupled with
sufficiently stiff quark EOSs to calculate structural properties of hybrid
stars that lie within the bounds of astrophysical measurements, (ii) a first
calculation of the two sound speeds and the principal $g$-mode frequency of
hybrid stars employing Gibbs phase criteria, and (iii) a concomitant analysis
of tidal phase effects in a binary merger due to $g$-modes in hybrid stars. In
future work, we aim to report on $g$-mode frequencies in alternative
treatments of quark matter in NSs such as a first order nucleon-to-quark phase
transition and crossover transitions as in quarkyonic matter.

_Acknowledgments._ — We are grateful to the anonymous referee for meticulous
review of the equations. We acknowledge discussions with Thomas Klähn on the
non-perturbative EOS for quark matter. Thanks are due to Sophia Han for
remarks on observational constraints on the EOS. P.J. is supported by the U.S.
NSF Grant No. PHY-1608959 and PHY-1913693. The research of A.S. and M.P. was
supported by the U.S. Department of Energy, Grant No. DE-FG02-93ER-40756. C.C.
acknowledges support from the European Union’s Horizon 2020 research and
innovation programme under the Marie Skłodowska-Curie grant agreement No.
754496 (H2020-MSCA-COFUND-2016 FELLINI).

##  Appendix A Determination of EOS constants in SNM and PNM for the ZL EOS

###  A.1 SNM

The constants $a_{0},b_{0}$ and $\gamma$ in Eq. (10) for SNM are determined by
utilizing the empirical properties of SNM at $u=1$. Specifically, the values
used are $E_{1/2}=-B=-16$ MeV at $n_{s}=0.16~{}{\rm fm^{-3}}$,
$p_{1/2}/n_{s}=0$, and $K_{1/2}=220$ MeV. Manipulating the relations

$\displaystyle-B$ $\displaystyle=$ $\displaystyle T_{1/2}+a_{0}+b_{0}$ (90)
$\displaystyle 0$ $\displaystyle=$ $\displaystyle
T_{1/2}^{\prime}+a_{0}+\gamma b_{0}$ (91) $\displaystyle\frac{K_{1/2}}{9}$
$\displaystyle=$ $\displaystyle
T_{1/2}^{\prime\prime}+2T_{1/2}^{\prime}+2a_{0}+\gamma(\gamma+1)b_{0}\,,$ (92)

the constants are given by

$\displaystyle\gamma$ $\displaystyle=$
$\displaystyle\frac{K_{1/2}/9-T_{1/2}^{\prime\prime}}{T_{1/2}-T_{1/2}^{\prime}+B}\,,\quad$
$\displaystyle b_{0}$ $\displaystyle=$
$\displaystyle\frac{K_{1/2}/9-T_{1/2}^{\prime\prime}}{\gamma(\gamma-1)}\quad{\rm
and}$ $\displaystyle a_{0}$ $\displaystyle=$ $\displaystyle-B-
T_{1/2}-b_{0}\,,$ (93)

where $T_{1/2}^{\prime}=u\frac{dT_{1/2}}{du}$ and
$T_{1/2}^{\prime\prime}=u^{2}\frac{d^{2}T_{1/2}}{du^{2}}$. Explicit
expressions for these derivatives are

$\displaystyle T_{1/2}^{\prime}$ $\displaystyle=$
$\displaystyle\left.\frac{p_{1/2}^{\rm kin}}{n}\right|_{n_{s}}$
$\displaystyle=$
$\displaystyle\frac{1}{n_{s}}\cdot\frac{2}{12\pi^{2}}\bigg{[}k_{F}E_{F}\left(k_{F}^{2}-\frac{3}{2}M_{B}^{2}\right)$
$\displaystyle+$
$\displaystyle\frac{3}{2}M_{B}^{4}\ln\left(\frac{k_{F}+E_{F}}{M_{B}}\right)\bigg{]}_{k_{Fs}}\,,$
$\displaystyle T_{1/2}^{\prime\prime}$ $\displaystyle=$
$\displaystyle\frac{K_{1/2}^{\rm
kin}}{9}-2T_{/12}^{\prime}=\frac{k_{Fs}^{2}}{3E_{Fs}}-2T_{/12}^{\prime},$ (94)

where $k_{Fs}=(3\pi^{2}n_{s}/2)^{1/3}$. To obtain the first term in the
rightmost equality above, it is advantageous to use the thermodynamical
identity $p=n\mu-\epsilon$ for the kinetic parts, whence
$\frac{dp}{dn}=n\frac{d\mu}{dn}=\frac{d\mu}{dk_{F}}\frac{dk_{F}}{dn}$. The
result quoted above ensues from the relations
$\frac{dk_{F}}{dn}=\frac{k_{F}}{3n}$ and
$\mu=E_{F}=\sqrt{k_{F}^{2}+M_{B}^{2}}$ both evaluated at $n_{s}$ and $k_{Fs}$.

Numerical values of the derivatives and constants so derived are

$\displaystyle T_{1/2}$ $\displaystyle\simeq$ $\displaystyle 21.79~{}{\rm
MeV},~{}~{}~{}T_{1/2}^{\prime}\simeq 14.34~{}{\rm MeV,}$ $\displaystyle
T_{1/2}^{\prime\prime}$ $\displaystyle=$ $\displaystyle-5.030~{}{\rm
MeV,}~{}~{}~{}\gamma\simeq 1.256,$ $\displaystyle~{}~{}~{}a_{0}$
$\displaystyle\simeq$ $\displaystyle-129.3~{}{\rm MeV},\quad{\rm
and}~{}~{}~{}b_{0}\simeq 91.49~{}{\rm MeV}.$ (95)

as in ZL. For other permissible values of $K_{1/2}$ in the range $220\pm 30$
MeV, Eqs. (A.1, A.1) can be used to determine the corresponding constants.

###  A.2 PNM

In the PNM sector in which $x=0$, the constants in Eq. (15) to determined are
$a_{1},b_{1}$ and $\gamma_{1}$. As in SNM, $E_{0}$ and $T_{0}$ are relative to
the baryon mass $M_{B}$. Denoting the energy per baryon of PNM by $E_{0}$, its
various terms and the associated pressure are

$\displaystyle E_{0}$ $\displaystyle=$ $\displaystyle
T_{0}+V_{0}=T_{0}+a_{1}u+b_{1}u^{\gamma}_{1}\,,$ $\displaystyle p_{0}$
$\displaystyle=$ $\displaystyle n_{s}\left(u^{2}\frac{dE_{0}}{du}\right)$ (96)
$\displaystyle=$ $\displaystyle
n_{s}\left(u^{2}T_{0}^{\prime}+a_{1}u^{2}+\gamma_{1}b_{1}u^{\gamma_{1}+1}\right)\,.$

Evaluating the above equations at $u=1$ leads to

$\displaystyle E_{0}=S_{v}-B=T_{0}+a_{1}+b_{1}\,,$ (97) $\displaystyle
p_{0}=n_{s}(T_{0}^{\prime}+a_{1}+\gamma_{1}b_{1})\,,$ (98)

where $S_{v}=(E_{0}-E_{1/2})$ at $u=1$. The last equation above is generally
written as

$\displaystyle\frac{p_{0}}{n_{s}}$ $\displaystyle=$
$\displaystyle\frac{L}{3}\quad{\rm with}\quad
L=3\left[n\frac{dS_{v}}{dn}\right]_{n_{s}}=3~{}[uS_{v}^{\prime}]_{u=1}\quad{\rm
so~{}that}$ $\displaystyle\frac{L}{3}$ $\displaystyle=$ $\displaystyle
T_{0}^{\prime}+a_{1}+\gamma_{1}b_{1}\,,$ (99)

where $S_{v}^{\prime}=\frac{dS_{v}}{du}$. Manipulating Eqs. (97) and (99)
leads to the relations

$\displaystyle b_{1}$ $\displaystyle=$
$\displaystyle\frac{\frac{L}{3}+B-S_{v}+T_{0}-T_{0}^{\prime}}{\gamma_{1}-1}\quad{\rm
and}\quad$ $\displaystyle a_{1}$ $\displaystyle=$ $\displaystyle
S_{v}-B-T_{0}-b_{1}\,.$ (100)

Taking guidance from the empirical properties of isospin asymmetric nuclear
matter, we choose $S_{v}=31$ MeV, $L$ in the range ($30$-$70$) MeV, and
$\gamma_{1}=5/3$. The resulting values of the constants are

$\displaystyle a_{1}$ $\displaystyle\simeq$
$\displaystyle-\left(\frac{L}{2}+14.72\right)~{}{\rm MeV}\quad{\rm and}\quad$
$\displaystyle b_{1}$ $\displaystyle\simeq$
$\displaystyle\left(\frac{L}{2}-4.62\right)~{}{\rm MeV}.$ (101)

###  A.3 Sensitivity of the EOS constants

The EOS constants above depend on the input values of $B,~{}n_{s},~{}K_{0}$
and $S_{v},~{}L,~{}\gamma_{1}$ in SNM and PNM, respectively. Although we have
used representative values for these quantities at $u=1$, nuclear data permits
variations in them. Furthermore, one or more sets of these constants may be
correlated, as for example, $S_{V}$ and $L$. Additional constraints are to
support $\simeq 2M_{\odot}$ NSs and to maintain causality, at least within the
star. These points must be borne in mind when varying the input values,
particularly when correlated errors are present in theoretical evaluations of
these quantities.

###  A.4 NS properties with the ZL EOSs

The various properties shown in Table 2 below are for beta-equilibrated normal
NSs and correspond to variations in the characteristic properties of the ZL
EOSs.

Table 2: Structural properties of nucleonic NSs with $M=1.4\,M_{\odot}$ and $M_{\rm max}$ for the ZL EOSs. For each mass, the compactness parameter $\beta=(GM/c^{2}R)\simeq(1.475/R)(M/M_{\odot})$, $n_{c}$, $p_{c}$ and $y_{c}$ are the central values of the density, pressure and proton fraction, respectively. The corresponding equilibrium squared speeds of sound are denoted by $c_{e}^{2}$. The $\Lambda$’s denote tidal deformabilities.  Property | ZL-A | ZL-B | ZL-C | Units
---|---|---|---|---
$K_{0}$ | 220 | 220 | 240 | MeV
$S_{v}$ | 31 | 31 | 31 | MeV
$L$ | 50 | 70 | 60 | MeV
$\gamma_{1}$ | 1.6 | 1.6 | 1.8 |
$R_{1.4}$ | 11.77 | 12.69 | 12.61 | km
$\beta_{1.4}$ | 0.175 | 0.163 | 0.164 |
$n_{c,1.4}/n_{s}$ | 3.35 | 2.78 | 2.75 |
$p_{c,1.4}$ | 83.65 | 60.59 | 61.50 | ${\rm MeV~{}fm^{-3}}$
$(c_{e}^{2})_{c,1.4}$ | 0.385 | 0.345 | 0.363 | $c^{2}$
$\Lambda_{1.4}$ | 713.4 | 970.2 | 504.1 |
$R_{\rm max}$ | 10.01 | 10.68 | 10.8 | km
$M_{\rm max}$ | 1.997 | 2.02 | 2.13 | $M_{\odot}$
$\beta_{\rm max}$ | 0.294 | 0.279 | 0.291 |
$n_{\rm c,max}/n_{s}$ | 7.71 | 6.96 | 6.67 |
$p_{\rm c,max}$ | 798.89 | 602.04 | 646.7 | ${\rm MeV~{}fm^{-3}}$
$(c_{e}^{2})_{\rm c,max}$ | 0.874 | 0.777 | 0.866 | $c^{2}$
$\Lambda_{\rm max}$ | 9.11 | 9.6 | 6.39 |

Structural properties of hybrid stars are discussed and shown in various
figures in the text.

###  A.5 Proof of equivalence of Eqs. (40) and (51)

Here we establish the equivalence of the direct approach to computing
$c_{s}^{2}$ from Eq. (40) with that from Eq. (51) for a general EOS with a
parabolic dependence of the proton fraction $x$ in the case of $n,p,e$ matter.
Both approaches yield identical results, which we have verified numerically as
well.

The equilibrium squared sound speed is

$c_{e}^{2}=\left(\frac{dP}{d\epsilon}\right)_{\rm
eq}=\frac{n_{B}\frac{d}{dn_{B}}\left(P_{\rm
bar}+P_{e}\right)}{\left(\epsilon+P\right)}\,,$ (102)

where the total pressure $P=P_{\rm bar}+P_{e}$ is comprised of the pressure
from baryons and electrons. Writing

$\frac{d}{dn_{B}}=\frac{\partial}{\partial
n_{B}}+\frac{dx}{dn_{B}}\frac{\partial}{\partial x}\,,$ (103)

we get

$c_{e}^{2}=\frac{\left(n_{B}\frac{\partial P_{bar}}{\partial
n_{B}}+n_{B}\frac{\partial P_{e}}{\partial
n_{B}}\right)}{\epsilon+P}+\frac{n_{B}\left(\frac{\partial P_{\text{bar
}}}{\partial x}\frac{dx}{dn_{B}}+\frac{\partial P_{\text{e}}}{\partial
x}\frac{dx}{dn_{B}}\right)}{\epsilon+P}\,.$ (104)

Comparing with Eq. (39), the first term on the right hand side is simply
$c_{s}^{2}$. For the specific case of an EOS with parabolic dependence in $x$
(of which the APR-EOS in Kantor & Gusakov [29] and the ZL-EOS in Zhao and
Lattimer [43] are examples), we have

$E(u,x)=E_{0}(u)+(1-2x)^{2}S_{2}(u);\quad u=n_{B}/n_{0}\,.$ (105)

where $n_{0}$ is the saturation density, $E_{0}$ the energy per baryons
(neutrons and protons) and $S_{2}$ the symmetry energy. Computing the pressure
and its derivatives with respect to $n_{B}$ and $x$ for the EOS in Eq. (105),
we find

$\frac{\partial P_{\text{bar}}}{\partial x}=-4n_{B}(1-2x)uS_{2}^{\prime}\,,$
(106)

where the prime on $S_{2}$ is with respect to $u$. Since (assuming massless
electrons) $\frac{\partial P_{\text{e}}}{\partial x}=\frac{1}{3}n_{B}\mu_{e}$,
it follows that

$c_{e}^{2}=c_{s}^{2}+\frac{1}{\mu_{n}}\left(\frac{\mu_{e}}{3}-4(1-2x)uS_{2}^{\prime}\right)n_{B}\frac{dx}{dn_{B}}\,.$
(107)

From the $\beta$-equilibrium condition

$\mu_{e}=4S_{2}(1-2x)\Rightarrow\mu_{e}/(4S_{2})=(1-2x)\,,$ (108)

we get upon differentiation

$\frac{1}{4S_{2}}\left(\frac{d\mu_{e}}{dn_{B}}\right)-\left(\frac{\mu_{e}}{4S_{2}^{2}}\right)\left(\frac{dS_{2}}{dn_{B}}\right)=-2\left(\frac{dx}{dn_{B}}\right)\,.$
(109)

Using

$\left(\frac{d\mu_{e}}{dn_{B}}\right)=\frac{\mu_{e}}{3x}\left(\frac{dx}{dn_{B}}\right)+\frac{\mu_{e}}{3n_{B}}\,,$
(110)

and solving for $dx/dn_{B}$ from Eq.(109), a minor rearrangement yields

$n_{B}\left(\frac{dx}{dn_{B}}\right)=\frac{-\left(\frac{\mu_{e}}{3}-4(1-2x)uS_{2}^{\prime}\right)}{\left(\frac{\mu_{e}}{3x}+8S_{2}\right)}\,.$
(111)

Putting together Eq. (111) with Eq. (107), we get

$\displaystyle
c_{e}^{2}=c_{s}^{2}-\frac{1}{\mu_{n}}\frac{\left(\frac{\mu_{e}}{3}-4(1-2x)uS_{2}^{\prime}\right)^{2}}{\left(\frac{\mu_{e}}{3x}+8S_{2}\right)}\,.$
(112)

Comparing Eq.(112) with Eqs. (51), (55) and (57) in the text, namely,

$\displaystyle c_{s}^{2}$ $\displaystyle=$ $\displaystyle
c_{e}^{2}+\frac{\left[n_{B}\left(\frac{\partial\tilde{\mu}}{\partial
n_{B}}\right)_{x}\right]^{2}}{\mu_{n}\left(\frac{\partial\tilde{\mu}}{\partial
x}\right)_{n_{B}}}$ (113) $\displaystyle
n_{B}\frac{\partial\tilde{\mu}}{\partial n_{B}}$ $\displaystyle=$
$\displaystyle\frac{\mu_{e}}{3}-4(1-2x)~{}uS_{2}^{\prime}$ (114)
$\displaystyle\frac{\partial\tilde{\mu}}{\partial x}$ $\displaystyle=$
$\displaystyle\frac{1}{3}\frac{\mu_{e}}{x}+8S_{2}(u)\,,$ (115)

we see that Eq. (51) is consistent with the direct definition of $c_{s}^{2}$
from Eq. (40) and that Eq. (51) applies in general for any form of the EOS
with a parabolic dependence in $x$, although in the text we arrived at Eqs.
(55) and (57) in the context of the ZL-EOS.

###  A.6 Derivation of Eqs. (59)-(61)

In $npe\mu$ matter, we choose the independent variables to be the baryon
density $n_{B}$, the lepton fraction $x$, and the muon fraction $y\equiv
x_{\mu}$. The electron fraction $x_{e}$ is the difference $x-y$.

The starting point for the speed-of-sound difference is

$c_{s}^{2}-c_{e}^{2}=\frac{1}{\mu_{avg}}\left.\frac{\partial P}{\partial
n_{B}}\right|_{x,y}-\frac{1}{\mu_{n}}\frac{dP}{dn_{B}}~{}.$ (116)

The total derivative of the pressure $P(n_{B},x,y)$ with respect to $n_{B}$ is
given by

$\frac{dP}{dn_{B}}=\left.\frac{\partial P}{\partial
n_{B}}\right|_{x,y}+\left.\frac{\partial P}{\partial
x}\right|_{n_{B},y}\frac{dx}{dn_{B}}+\left.\frac{\partial P}{\partial
y}\right|_{n_{B},x}\frac{dy}{dn_{B}}$ (117)

and therefore

$\displaystyle c_{s}^{2}$ $\displaystyle-$ $\displaystyle
c_{e}^{2}=\left(\frac{1}{\mu_{avg}}-\frac{1}{\mu_{n}}\right)\left.\frac{\partial
P}{\partial n_{B}}\right|_{x,y}$ $\displaystyle-$
$\displaystyle\frac{1}{\mu_{n}}\left(\left.\frac{\partial P}{\partial
x}\right|_{n_{B},y}\frac{dx}{dn_{B}}+\left.\frac{\partial P}{\partial
y}\right|_{n_{B},x}\frac{dy}{dn_{B}}\right)$ $\displaystyle=$
$\displaystyle\left(\frac{\mu_{n}-\mu_{avg}}{\mu_{avg}~{}\mu_{n}}\right)\left.\frac{\partial
P}{\partial n_{B}}\right|_{x,y}$ $\displaystyle-$
$\displaystyle\frac{1}{\mu_{n}}\left(\left.\frac{\partial P}{\partial
x}\right|_{n_{B},y}\frac{dx}{dn_{B}}+\left.\frac{\partial P}{\partial
y}\right|_{n_{B},x}\frac{dy}{dn_{B}}\right).$ (119)

The average chemical potential $\mu_{avg}$ is

$\mu_{avg}=(1-x)\mu_{n}+x\mu_{p}+(x-y)\mu_{e}+y\mu_{y}$ (120)

which means that

$\displaystyle\mu_{n}-\mu_{avg}$ $\displaystyle=$ $\displaystyle
x(\mu_{n}-\mu_{p}-\mu_{e})+y(\mu_{e}-\mu_{\mu})$ (121) $\displaystyle\equiv$
$\displaystyle-x\tilde{\mu}_{x}-y\tilde{\mu}_{y}$

with the obvious definitions for $\tilde{\mu}_{x}$ and $\tilde{\mu}_{y}$.

In $\beta$-equilibrium $\mu_{n}=\mu_{avg}$ (as well as
$\tilde{\mu}_{x}=\tilde{\mu}_{y}=0$) and, correspondingly,

$c_{s}^{2}-c_{e}^{2}=-\frac{1}{\mu_{n}}\left(\left.\frac{\partial P}{\partial
x}\right|_{n_{B},y}\frac{dx}{dn_{B}}+\left.\frac{\partial P}{\partial
y}\right|_{n_{B},x}\frac{dy}{dn_{B}}\right)~{},$ (122)

which is Eq. (59) in the main text. Using $P=n_{B}^{2}\left.\frac{\partial
E}{\partial n_{B}}\right|_{x,y}$, the speed-of-sound difference can be
expressed as

$\displaystyle c_{s}^{2}-c_{e}^{2}$ $\displaystyle=$
$\displaystyle-\frac{n_{B}^{2}}{\mu_{n}}\frac{\partial}{\partial
n_{B}}\left.\left(\left.\frac{\partial E}{\partial
x}\right|_{n_{B},y}\frac{dx}{dn_{B}}+\left.\frac{\partial E}{\partial
y}\right|_{n_{B},x}\frac{dy}{dn_{B}}\right)\right|_{x,y}$ (123)
$\displaystyle=$
$\displaystyle-\frac{n_{B}^{2}}{\mu_{n}}\left(\left.\frac{\partial\tilde{\mu}_{x}}{\partial
n_{B}}\right|_{x,y}\frac{dx}{dn_{B}}+\left.\frac{\partial\tilde{\mu}_{y}}{\partial
n_{B}}\right|_{x,y}\frac{dy}{dn_{B}}\right)~{}.$

The calculation of $\frac{dx}{dn_{B}}$ and $\frac{dy}{dn_{B}}$ begins from the
total differentials of $\tilde{\mu}_{x}$ and $\tilde{\mu}_{y}$ which are

$d\tilde{\mu}_{x}=\left.\frac{\partial\tilde{\mu}_{x}}{\partial
n_{B}}\right|_{x,y}dn_{B}+\left.\frac{\partial\tilde{\mu}_{x}}{\partial
x}\right|_{n_{B},y}dx+\left.\frac{\partial\tilde{\mu}_{x}}{\partial
y}\right|_{n_{B},x}dy=0$ (124)

and

$d\tilde{\mu}_{y}=\left.\frac{\partial\tilde{\mu}_{y}}{\partial
n_{B}}\right|_{x,y}dn_{B}+\left.\frac{\partial\tilde{\mu}_{y}}{\partial
x}\right|_{n_{B},y}dx+\left.\frac{\partial\tilde{\mu}_{y}}{\partial
y}\right|_{n_{B},x}dy=0~{}.$ (125)

From the former differential, it follows that

$dy=\frac{-1}{\left.\frac{\partial\tilde{\mu}_{x}}{\partial
y}\right|_{n_{B},x}}\left(\left.\frac{\partial\tilde{\mu}_{x}}{\partial
n_{B}}\right|_{x,y}dn_{B}+\left.\frac{\partial\tilde{\mu}_{x}}{\partial
x}\right|_{n_{B},y}dx\right)$ (126)

which, when substituted in the latter, leads to

$\displaystyle 0$ $\displaystyle=$
$\displaystyle\left.\frac{\partial\tilde{\mu}_{y}}{\partial
n_{B}}\right|_{x,y}dn_{B}+\left.\frac{\partial\tilde{\mu}_{y}}{\partial
x}\right|_{n_{B},y}dx$ (127) $\displaystyle-$
$\displaystyle\frac{\left.\frac{\partial\tilde{\mu}_{y}}{\partial
y}\right|_{n_{B},x}}{\left.\frac{\partial\tilde{\mu}_{x}}{\partial
y}\right|_{n_{B},x}}\left(\left.\frac{\partial\tilde{\mu}_{x}}{\partial
n_{B}}\right|_{x,y}dn_{B}+\left.\frac{\partial\tilde{\mu}_{x}}{\partial
x}\right|_{n_{B},y}dx\right)~{}.$

One then collects terms proportional to $dn_{B}$ and $dx$

$\displaystyle\left(\left.\frac{\partial\tilde{\mu}_{y}}{\partial
n_{B}}\right|_{x,y}-\frac{\left.\frac{\partial\tilde{\mu}_{y}}{\partial
y}\right|_{n_{B},x}}{\left.\frac{\partial\tilde{\mu}_{x}}{\partial
y}\right|_{n_{B},x}}\left.\frac{\partial\tilde{\mu}_{x}}{\partial
n_{B}}\right|_{x,y}\right)dn_{B}$ (128) $\displaystyle=$
$\displaystyle-\left(\left.\frac{\partial\tilde{\mu}_{y}}{\partial
x}\right|_{n_{B},y}-\frac{\left.\frac{\partial\tilde{\mu}_{y}}{\partial
y}\right|_{n_{B},x}}{\left.\frac{\partial\tilde{\mu}_{x}}{\partial
y}\right|_{n_{B},x}}\left.\frac{\partial\tilde{\mu}_{x}}{\partial
x}\right|_{n_{B},y}\right)dx$

or, equivalently,

$\frac{dx}{dn_{B}}=\frac{\left.\frac{\partial\tilde{\mu}_{x}}{\partial
y}\right|_{n_{B},x}\left.\frac{\partial\tilde{\mu}_{y}}{\partial
n_{B}}\right|_{x,y}-\left.\frac{\partial\tilde{\mu}_{y}}{\partial
y}\right|_{n_{B},x}\left.\frac{\partial\tilde{\mu}_{x}}{\partial
n_{B}}\right|_{x,y}}{\left.\frac{\partial\tilde{\mu}_{x}}{\partial
y}\right|_{n_{B},x}\left.\frac{\partial\tilde{\mu}_{y}}{\partial
x}\right|_{n_{B},y}-\left.\frac{\partial\tilde{\mu}_{y}}{\partial
y}\right|_{n_{B},x}\left.\frac{\partial\tilde{\mu}_{x}}{\partial
x}\right|_{n_{B},y}}~{}.$ (129)

Similarly,

$\frac{dy}{dn_{B}}=\frac{\left.\frac{\partial\tilde{\mu}_{x}}{\partial
x}\right|_{n_{B},y}\left.\frac{\partial\tilde{\mu}_{y}}{\partial
n_{B}}\right|_{x,y}-\left.\frac{\partial\tilde{\mu}_{y}}{\partial
x}\right|_{n_{B},y}\left.\frac{\partial\tilde{\mu}_{x}}{\partial
n_{B}}\right|_{x,y}}{\left.\frac{\partial\tilde{\mu}_{x}}{\partial
x}\right|_{n_{B},y}\left.\frac{\partial\tilde{\mu}_{y}}{\partial
y}\right|_{n_{B},x}-\left.\frac{\partial\tilde{\mu}_{y}}{\partial
x}\right|_{n_{B},y}\left.\frac{\partial\tilde{\mu}_{x}}{\partial
y}\right|_{n_{B},x}}~{}.$ (130)

The speed-of-sound difference, as given by Eqs. (123), (129), and (130) is
physically transparent because $\beta$-equilibrium and compositional gradients
are brought to the forefront via $\tilde{\mu}_{x(y)}$ and $\partial/\partial
x(y)$, respectively (the latter two equations are Eqs. (60) and (61) in the
main text). However, this intuitive picture comes at the expense of
computational complexity.

## References

  * Collins and Perry [1975] J. C. Collins and M. J. Perry, Phys. Rev. Lett. 34, 1353 (1975).
  * Kurkela  _et al._ [2010] A. Kurkela, P. Romatschke, and A. Vuorinen, Phys. Rev. D 81, 105021 (2010).
  * Kurkela and Vuorinen [2016] A. Kurkela and A. Vuorinen, Phys. Rev. Lett. 117, 042501 (2016).
  * Gorda  _et al._ [2018] T. Gorda, A. Kurkela, P. Romatschke, M. Säppi, and A. Vuorinen, Phys. Rev. Lett. 121, 202701 (2018).
  * Annala  _et al._ [2020] E. Annala, T. Gorda, A. Kurkela, J. Nättilä, and A. Vuorinen, Nature Physics 16, 907 (2020).
  * Han  _et al._ [2019] S. Han, M. A. A. Mamun, S. Lalit, C. Constantinou, and M. Prakash, Phys. Rev. D 100, 103022 (2019).
  * Page  _et al._ [2015] D. Page, J. M. Lattimer, M. Prakash, and A. W. Steiner, Novel Superfluids, International Series of Monographs on Physics, edited by K. H. Bennemann, and J. B. Ketterson (Oxford University Press) , 505 (2015).
  * Potekhin  _et al._ [2015] A. Y. Potekhin, J. A. Pons, and D. Page, Space Sci. Rev. 191, 239 (2015).
  * Hinderer [2008] T. Hinderer, The Astrophysical Journal 677, 1216 (2008).
  * Damour and Nagar [2009] T. Damour and A. Nagar, Phys. Rev. D 80, 084035 (2009).
  * Postnikov  _et al._ [2010] S. Postnikov, M. Prakash, and J. M. Lattimer, Phys. Rev. D 82, 024016 (2010).
  * Abbott  _et al._ [2017] B. P. Abbott, others for the LIGO Scientific Collaboration, and Virgo Collaboration, Phys. Rev. Lett. 119, 161101 (2017).
  * Parisi  _et al._ [2020] A. Parisi, C. Vásquez Flores, C. H. Lenzi, C.-S. Chen, and G. Lugones, arXiv e-prints , arXiv:2009.14274 (2020).
  * Tews  _et al._ [2019] I. Tews, J. Margueron, and S. Reddy, European Physical Journal A 55, 97 (2019).
  * Kokkotas  _et al._ [2001] K. D. Kokkotas, T. A. Apostolatos, and N. Andersson, Mon. Not. Roy. Ast. Soc. 320, 307 (2001).
  * Alford  _et al._ [2019] M. G. Alford, S. Han, and K. Schwenzer, Journal of Physics G Nuclear Physics 46, 114001 (2019).
  * Reisenegger and Goldreich [1992] A. Reisenegger and P. Goldreich, Astrophys. J. 395, 240 (1992).
  * Reisenegger and Goldreich [1994] A. Reisenegger and P. Goldreich, Astrophys. J. 426, 688 (1994).
  * Cox [1980] J. P. Cox,  _Theory of Stellar Pulsation_ (Princeton University Press, 1980).
  * Lai [1994] D. Lai, Mon. Not. Roy. Astron. Soc. 270, 611 (1994).
  * Finn [1987] L. S. Finn, Mon. Not. Roy. Ast. Soc. 227, 265 (1987).
  * Strohmayer [1993] T. E. Strohmayer, Astrophys. J. 417, 273 (1993).
  * McDermott [1990] P. N. McDermott, Mon. Not. Roy. Ast. Soc. 245, 508 (1990).
  * Lee [1995] U. Lee, Astron. Astrophys. 303, 515 (1995).
  * Prix and Rieutord [2002] R. Prix and M. Rieutord, Astron. & Astrophys. 393, 949–963 (2002).
  * Andersson and Comer [2001] N. Andersson and G. L. Comer, Mon. Not. Roy. Ast. Soc. 328, 1129 (2001).
  * Gusakov and Kantor [2013] M. E. Gusakov and E. M. Kantor, Phys. Rev. D 88, 101302 (2013).
  * Gualtieri  _et al._ [2014] L. Gualtieri, E. M. Kantor, M. E. Gusakov, and A. I. Chugunov, Phys. Rev. D 90, 024010 (2014).
  * Kantor and Gusakov [2014] E. M. Kantor and M. E. Gusakov, Mon. Not. Roy. Astron. Soc. 442, 90 (2014).
  * Dommes and Gusakov [2016] V. A. Dommes and M. E. Gusakov, Mon. Not. Roy. Astron. Soc. 455, 2852 (2016).
  * Passamonti  _et al._ [2016] A. Passamonti, N. Andersson, and W. C. G. Ho, Mon. Not. Roy. Astron. Soc. 455, 1489 (2016).
  * Yu and Weinberg [2017] H. Yu and N. N. Weinberg, Mon. Not. Roy. Astron. Soc. 470, 350 (2017).
  * Yu and Weinberg [2017] H. Yu and N. N. Weinberg, Mon. Not. Roy. Ast. Soc. 464, 2622 (2017).
  * Rau and Wasserman [2018] P. B. Rau and I. Wasserman, Mon. Not. Roy. Ast. Soc. 481, 4427 (2018).
  * Gusakov  _et al._ [2014] M. E. Gusakov, P. Haensel, and E. M. Kantor, Monthly Notices of the Royal Astronomical Society 439, 318 (2014).
  * Wei  _et al._ [2020] W. Wei, M. Salinas, T. Klähn, P. Jaikumar, and M. Barry, The Astrophysical Journal 904, 187 (2020).
  * Abbott  _et al._ [2018] B. P. Abbott _et al._ , Phys. Rev. Lett. 121, 161101 (2018).
  * Riley  _et al._ [2019] T. E. Riley _et al._ , The Astrophysical Journal Letters 887, L21 (2019).
  * Miller  _et al._ [2019] M. C. Miller _et al._ , The Astrophysical Journal Letters 887, L24 (2019).
  * Bedaque and Steiner [2015] P. Bedaque and A. W. Steiner, Phys. Rev. Lett. 114, 031103 (2015).
  * Kanakis-Pegios  _et al._ [2020] A. Kanakis-Pegios, P. S. Koliogiannis, and C. C. Moustakidis, Phys. Rev. C 102, 055801 (2020).
  * McLerran and Reddy [2019] L. McLerran and S. Reddy, Phys. Rev. Lett. 122, 122701 (2019).
  * Zhao and Lattimer [2020] T. Zhao and J. M. Lattimer, Phys. Rev. D 102, 023021 (2020).
  * Brillante and Mishustin [2014] A. Brillante and I. N. Mishustin, EPL 105, 39001 (2014).
  * Lau and Yagi [2020] S. Y. Lau and K. Yagi, arXiv e-prints (2020), arXiv:2012.13000 [astro-ph.HE] .
  * Jaikumar  _et al._ [2008] P. Jaikumar, G. Rupak, and A. W. Steiner, Phys. Rev. D 78, 123007 (2008).
  * Sa’d [2008] B. A. Sa’d, arXiv e-prints (2008), arXiv:0806.3359 [astro-ph] .
  * Rupak and Jaikumar [2010] G. Rupak and P. Jaikumar, Phys. Rev. C 82, 055806 (2010).
  * Sotani  _et al._ [2013] H. Sotani, T. Maruyama, and T. Tatsumi, Nucl. Phys. A906, 37 (2013).
  * Flores and Lugones [2014] C. V. Flores and G. Lugones, Class. Quant. Grav. 31, 155002 (2014).
  * Ranea-Sandoval  _et al._ [2018] I. F. Ranea-Sandoval, O. M. Guilera, M. Mariani, and M. G. Orsaria, Journal of Cosmology and Astroparticle Physics 2018 (12), 031.
  * Heiselberg [1995] H. Heiselberg, in  _Strangeness and Quark Matter_ , edited by G. Vassiliadis, A. D. Panagiotou, S. Kumar, and J. Madsen (1995) pp. 298–307.
  * Xia  _et al._ [2018] C.-J. Xia, G.-X. Peng, T.-T. Sun, W.-L. Guo, D.-H. Lu, and P. Jaikumar, Phys. Rev. D 98, 034031 (2018).
  * Gregorian [2015] P. Gregorian, Masters Thesis, Eberhard Karls Universität Tübingen (2015), uRL:https://dspace.library.uu.nl/handle/1874/306758.
  * McDermott  _et al._ [1983] P. N. McDermott, H. M. van Horn, and J. F. Scholl, Astrophys. J. 268, 837 (1983).
  * Chandrasekhar and Ferrari [1991] S. Chandrasekhar and V. Ferrari, Proceedings of the Royal Society of London. Series A: Mathematical and Physical Sciences 432, 247 (1991).
  * Lai [1999] D. Lai, Mon. Not. Roy. Astron. Soc. 307, 1001 (1999).
  * Drischler  _et al._ [2016] C. Drischler, A. Carbone, K. Hebeler, and A. Schwenk, Phys. Rev. C 94, 054307 (2016).
  * Lattimer and Prakash [2011] J. M. Lattimer and M. Prakash, From Nuclei to Stars, ed. S. Lee, Singapore: WorldScientific , 275 (2011).
  * Gomes  _et al._ [2019] R. O. Gomes, P. Char, and S. Schramm, The Astrophysical Journal 877, 139 (2019).
  * Glendenning [1992] N. K. Glendenning, Phys. Rev. D 46, 1274 (1992).
  * Baym and Chin [1976] G. Baym and S. A. Chin, Physics Letters B 62, 241 (1976).
  * Klähn and Fischer [2015] T. Klähn and T. Fischer, Astrophys. J. 810, 134 (2015).
  * Tews  _et al._ [2018] I. Tews, J. Carlson, S. Gandolfi, and S. Reddy, The Astrophysical Journal 860, 149 (2018).
  * Haensel  _et al._ [2000] P. Haensel, K. P. Levenfish, and D. G. Yakovlev, Astron. & Astrophys. 357, 1157 (2000), arXiv:astro-ph/0004183 [astro-ph] .
  * Abney  _et al._ [1996] M. Abney, R. I. Epstein, and A. V. Olinto, Astrophys. J. Lett. 466, L91 (1996).
  * Lattimer [2012] J. M. Lattimer, Annual Review of Nuclear and Particle Science 62, 485 (2012).
  * Özel and Freire [2016] F. Özel and P. Freire, Ann. Rev. Astron. & Astrophys. 54, 401 (2016).
  * Özel  _et al._ [2016] F. Özel, D. Psaltis, T. Güver, G. Baym, C. Heinke, and S. Guillot, The Astrophysical Journal 820, 28 (2016).
  * De  _et al._ [2018] S. De, D. Finstad, J. M. Lattimer, D. A. Brown, E. Berger, and C. M. Biwer, Phys. Rev. Lett. 121, 091102 (2018).
  * Capano  _et al._ [2020] C. D. Capano, I. Tews, S. M. Brown, B. Margalit, S. De, S. Kumar, D. A. Brown, B. Krishnan, and S. Reddy, Nature Astronomy 4, 625 (2020).
  * Demorest  _et al._ [2010] P. B. Demorest, T. Pennucci, S. M. Ransom, M. S. E. Roberts, and J. W. T. Hessels, Nature (London) 467, 1081 (2010).
  * Fonseca  _et al._ [2016] E. Fonseca, T. T. Pennucci, J. A. Ellis, _et al._ , Astrophys. J. 832, 167 (2016).
  * Arzoumanian  _et al._ [2018] Z. Arzoumanian, A. Brazier, S. Burke-Spolaor, _et al._ , The Astrophysical Journal Supplement Series 235, 37 (2018).
  * Cromartie  _et al._ [2020] H. T. Cromartie _et al._ , Nature Astronomy 4, 72 (2020).
  * Linares  _et al._ [2018] M. Linares, T. Shahbaz, and J. Casares, The Astrophysical Journal 859, 54 (2018).
  * Flores  _et al._ [2017] C. V. Flores, Z. B. Hall, II, and P. Jaikumar, Phys. Rev. C96, 065803 (2017).
  * Wu and Shen [2019] X. H. Wu and H. Shen, Phys. Rev. C 99, 065802 (2019).
  * Maggiore [2007] M. Maggiore,  _Gravitational Waves Vol 1: Theory and Experiments_ (Oxford University Press, 2007).