# Cosmic twins are likely rare: exploring pathways to hadron-quark stars with
the same mass as neutron stars

Mahdi Naseri<EMAIL_ADDRESS>Department of Astronomy, University of
Arizona, Tucson, AZ, USA  Gabriele Bozzola<EMAIL_ADDRESS>Department of
Astronomy, University of Arizona, Tucson, AZ, USA
Division of Geological and Planetary Sciences, California Institute of
Technology, Pasadena, CA, USA  Vasileios Paschalidis<EMAIL_ADDRESS>Departments of Astronomy and Physics, University of Arizona, Tucson, AZ, USA

(August 31, 2024)

###### Abstract

A viable model for the dense matter equation of state above the nuclear
saturation density is a hadron-to-quark phase transition at densities relevant
to compact objects. In this case, stable hybrid hadron-quark stars can arise.
An even more interesting scenario is one where the hadron-to-quark phase
transition results in the emergence of a third branch of stable compact
objects (in addition to white dwarfs and neutron stars). Inherent to the
presence of a third family of compact stars is the existence of twin stars –
hybrid stars with the same mass as the corresponding neutron stars, but with
smaller radii. Interestingly, the neutron star-twin star scenario is
consistent with GW170817. If twin stars exist in nature, it raises a question
about the mechanism that leads to their formation. Here, we explore
gravitational collapse as a pathway to the formation of low-mass twin stars.
We perform fully general relativistic simulations of the collapse of a stellar
iron core, modeled as a cold degenerate gas, to investigate whether the end
product is a neutron star or a twin star. Our simulations show that even with
unrealistically large perturbations in the initial conditions, the core
bounces well below the hadron-to-quark phase transition density, if the
initial total rest mass is in the twin star range. Following cooling, these
configurations produce neutron stars. We find that twin stars can potentially
form due to mass loss, e.g., through winds, from a slightly more massive
hybrid star that was initially produced in the collapse of a more massive core
or if the maximum NS mass is below the Chandrasekhar mass limit. The challenge
in producing twin stars in gravitational collapse, in conjunction with the
fine tuning required because of their narrow mass range, suggest the rarity of
twin stars in nature.

††preprint: APS/123-QED

##  I Introduction

One of the most important open questions in nuclear physics research is the
state of matter above the nuclear saturation density, $\rho_{s}\approx
2.7\times 10^{14}$ g cm-3. The equation of state (EOS) of cold neutron-rich
matter is reasonably well understood at densities below $\sim 2\,\rho_{s}$ [1,
2, 3] (see also [4, 5, 6, 7, 8, 9, 10, 11, 12, 13]). Uncertainties exist even
for finite temperature dense nuclear matter, see e.g., [14, 15, 16].
Theoretical understanding at these densities is aided by laboratory
experiments involving heavy ion collisions, see e.g., [17, 18, 19, 20, 21, 22,
23]. However, the dense-matter cold EOS remains highly uncertain at densities
of $\sim 2-40\,\rho_{s}$ because of limitations in theoretical or lattice
quantum chromodynamic (QCD) approaches, see, e.g., [24, 25, 26, 27, 3]. As a
result, a plethora of theoretical possibilities for how matter behaves at
densities higher than $2\,\rho_{s}$ is currently allowed. In this work we
focus on the possibility that a quark deconfinement phase transition can take
place at densities encountered inside compact objects.

One of the fundamental consequences of QCD is color confinement. In the
absence of strong medium effects, quarks and gluons are confined inside
hadrons. Quark confinement can be violated at high temperatures or high
densities, leading to the formation of a new state of matter known as quark-
gluon plasma [28, 29]. It is currently unclear whether quark deconfinement can
take place in the deep interiors of compact objects, but this possibility has
been considered by a number of works (see, e.g., [30, 31, 32, 33, 34, 35, 36,
37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 38, 37, 47, 48, 44, 49, 50].

Figure 1:  Mass-central density plot for a sequence of TOV stars based on the
EOS detailed in Table 1. Stars with a fixed rest mass $M_{0}=M_{0,s}\sim
M_{ch}$ are designated by labels A, B, C, D, and E. The inset zooms in on the
twin star part of the curve.

If quark deconfinement takes place at astrophysically relevant densities, an
interesting scenario arises when the phase transition occurs such that a
“third family” of compact objects (in addition to white dwarfs and neutron
stars) emerges [51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63]. Compact
objects in this third family are hybrid hadron-quark stars or hybrid stars –
configurations with a quark core enveloped by a hadronic shell. The stable
third family of compact stars is separated from the branch of stable neutron
stars (NSs) by an unstable branch of compact objects, just like stable NSs are
separated from stable white dwarfs (WDs) by an unstable branch. Fig. 1 depicts
these branches on a mass versus central density plot (further details on this
plot are discussed in Sec. II). As a result, there exists a set of stable
hybrid stars (HSs) that have the same mass as NSs but smaller radii. These HSs
are typically referred to as twin stars (TSs) [54, 64, 65, 66, 67, 68, 69, 70,
71]. Recent works suggest that current experimental and multimessenger
observations of NSs are either consistent with or provide evidence in favor of
the existence of hadron-to-quark phase transitions in their cores [72, 73, 74,
75, 76, 77], but see [78, 79]. The third family of compact objects with
potentially low-mass TSs has been shown to be consistent with a number of
experiments and observations [34, 80, 81, 67, 82, 83, 84, 85]. Twin stars are
particularly interesting because proof of their existence would provide
valuable insights into the quark-gluon phase diagram.

The third family can arise when the EOS features a strong hadron-quark phase
transition (i.e. characterized by a sufficiently large jump in energy density
at a roughly fixed pressure) which is first-order (meaning that the first
derivative of free energy with respect to a thermodynamic variable is
discontinuous) [53]. Alternatively, hybrid EOSs with mixed pasta phases can
also lead to the formation of the third family, e.g., [86]. What determines
the nature of the phase transition is the surface tension between the quark
and hadron phases [87, 88]. If the tension between these phases is low, a
mixed phase of quark and nucleonic matter forms in-between purely nuclear and
quark matter phases. By contrast, if the tension is high, a sharp transition
boundary is favorable. Both possibilities are theoretically allowed, because
the surface tension is not known accurately.

Given the large number of works suggesting that hybrid and twin stars are
compatible with observations, an important question arises: if these stars
exist, how can they form? We begin to tackle this question in this work.

Stellar gravitational collapse appears to be a natural pathway for the
formation of HSs. One might expect this, given that WDs and NSs are endpoints
of stellar evolution, see e.g., [89, 90]. However, forming TSs through
gravitational collapse faces an important challenge; if one imagines the
collapse of stellar core at approximately constant rest mass, following
cooling, the configuration would contract and encounter first the branch of
stable NSs, and then the branch of TSs, assuming further contraction. Thus, it
is unclear whether there exists a natural pathway during stellar collapse that
would prefer to form TSs over NSs. In fact, the work of [91] treated the
evolution of configurations in the TS mass range that are initially on the
unstable branch, and found that these stars naturally transition to the stable
NS branch. Thus, Ref. [91] concluded that it would be challenging to form TSs
through gravitational collapse because in a quasi-adiabatic, constant rest-
mass contraction of a proto-NS that is on its way to forming a TS, it would
have to go through the unstable branch first, which appears to favor NSs.

The primary objective of this paper is to explore whether HSs, and in
particular, TSs, can form as a result of the gravitational collapse of a
stellar core. To investigate this pathway toward HSs formation, we perform
hydrodynamic simulations of unstable WDs in full general relativity by
considering several types of initial perturbations differing in the degree of
violence of the collapse. The advantage of considering a WD for our initial
conditions is that it models two distinct channels to forming neutron/hybrid
stars: i) the collapse of a degenerate iron core of a massive star, and ii)
the accretion-induced collapse of a WD [92]. Here, we do not treat detailed
microphysics or neutrino effects. Instead, we focus on treating properly
relativistic gravitation, which is critical when compact objects arise.

Our simulations demonstrate that for a type of EOS which gives rise to a third
family of compact objects, the gravitational collapse of a WD in the twin star
mass regime prefers the formation of NSs and not TSs, even under extreme
initial perturbations. The only pathway to forming TSs that we were able to
identify is as follows: first, a massive core should collapse to a slightly
more massive stable HS, which may subsequently lose a small amount of mass,
for example, in the form of winds or due to rotation or a “grazing” collision
with a black hole, to ultimately settle as a TS. The results of our
simulations, in conjunction with the narrow mass range over which twin stars
exist suggest that twin stars should be rare. Thus, if a HS star was involved
in GW170817, then it likely was not a twin star.

The paper is structured as follows. In Sec. II, we construct the EOS that we
employ in this work based on a realistic EOS with a hadron-to-quark phase
transition. In Sec. III, we discuss our initial data and evolution methods
along with the numerical scheme employed in the simulations. Results of
various simulations exploring the formation of HSs under different conditions
are presented and discussed in Sec. IV. We conclude in Sec. V, where we
summarize our main findings. Unless specified otherwise, throughout this paper
we adopt geometrized units where $c=G=1$, with $c$ being the speed of light in
vacuum and $G$ representing the gravitational constant.

##  II Equations of State

There exists a broad range of possibilities for the EOS in the deep interior
of an NS. Numerous studies have explored different hadronic models, quark
models and hybrid hadron-quark models. For comprehensive reviews readers are
referred to [27, 93, 94, 95]. Apart from microscopic, so-called realistic
EOSs, the EOS can also be treated phenomenologically, e.g., by giving the
pressure or sound speed as a function of rest-mass or energy density, see,
e.g., [96, 97, 72, 98]. Here, we restrict our discussion to the astrophysical
implications of EOSs that are based on realistic models but treated
phenomenologically. The EOSs we consider give rise to a third family of stable
compact stars. However, to simulate the collapse of a stellar iron core or a
WD, the EOSs we construct must encompass a wide range of densities, spanning
the entire range from below neutron drip to the supranuclear regime. Here we
describe how we construct the phenomenological EOSs we adopt.

The EOS below the nuclear saturation density down to WD densities is
reasonably well understood. Our base phenomenological EOS selected for this
density regime is derived from the 6-parameter piecewise polytropic EOS
introduced in [99]. For the high-density regime, the EOS needs to incorporate
a phase transition, so that a stable HS branch exists. The high-density EOS
that we adopt here is based on the T9 EOS of [91], which is a piecewise
polytropic representation of the T9 EOS used in [36], which was in turn based
on the ACS-II with $\xi=0.90$ EOS of [34]. This EOS describes zero-temperature
matter in beta-equilibrium. Its hadronic part is derived from [100], the quark
phase builts upon the MIT bag model [101, 102, 103, 104], and the low-density
crust component is added from the models in [105, 106]. The quark phase of the
original EOS adopts the constant sound speed parametrization

$P(\epsilon)=\begin{cases}P_{tr}&\epsilon_{q}\leq\epsilon\leq\epsilon_{tr}\\\
P_{tr}+c_{s}^{2}(\epsilon-\epsilon_{tr})&\epsilon\geq\epsilon_{tr}\;,\end{cases}$
(1)

where $\epsilon$ is the energy density,
$P_{tr}=P(\epsilon_{q})=P(\epsilon_{tr})$ is the pressure of the hadronic
matter in the transition region $\epsilon\in[\epsilon_{q},\epsilon_{tr}]$, and
$c_{s}$ denotes the sound speed.

Given that the original T9 EOS is provided in tabulated form, we represent the
high-density EOS, similar to the low-density regime, as a piecewise polytrope
of the form

$P=k_{i}\rho_{0}^{\Gamma_{i}}\;,$ (2)

where $\rho_{0}$ is rest-mass density, and
$\rho_{0,i}\leq\rho_{0}\leq\rho_{0,{i+1}}$ is the density range of each
polytropic segment. Piecewise polytropes are frequently employed in the
context of the NS EOS, as many proposed tabulated, realistic nuclear EOSs can
be well-fitted by them, see, e.g., [96, 7, 107, 98].

The next step is to merge these two EOSs into one that can describe the entire
range of densities relevant to compact objects. The low-density and high-
density EOSs represented as $P=P(\rho_{0})$ intersect at a matching density.
The combined EOS should initially follow the low-density EOS and transition to
the high-density EOS as the density increases beyond the matching density. To
ensure a smooth transition below and above the matching density, in the
density range of approximately $10^{11}-10^{12}\textrm{g cm}^{-3}$, we
linearly interpolate between the two regimes using logarithmic pressure and
density. This step yields a base EOS that covers a wide range of densities,
extending from the crust of a WD to the dense core of a HS.

After we join the low- and high-density EOSs, we solve the Tolman-Oppenheimer-
Volkoff (TOV) equations [89] for determining general relativistic hydrostatic
equilibrium stellar configurations. We then make small modifications to the
EOSs to satisfy the following set of conditions:

  1. 1.

A third family of compact objects is present.

  2. 2.

A $2M_{\odot}$ [108, 109, 110, 111, 112] compact object should be allowed by
the EOS.

  3. 3.

The Chandrasekhar mass $M_{ch}=1.44M_{\odot}$ is the maximum WD mass, assuming
a mean molecular weight of electrons $\mu_{e}=2$.

  4. 4.

For the range of densities in our simulations, the sound speed should be
subluminal. The sound speed for a piecewise polytropic EOS is calculated as
[96]

$c_{s}=\sqrt{\frac{dP}{d\epsilon}}=\sqrt{\frac{\Gamma_{i}P}{\epsilon+P}}\;.$
(3)

  5. 5.

In the phase transition region, the sound speed should not be zero because the
equations of hydrodynamics would be only weakly hyperbolic (see [91] and
discussion therein).

  6. 6.

The dominant energy condition should be satisfied for all densities. In the
case of a perfect fluid, this implies [113]

$\epsilon\geq|P|\;.$ (4)

Following small modifications to the piecewise polytropic parameters of our
EOS such that all aforementioned conditions are met we constructed the EOSs
presented in Table 1. The EOS designated as “EOS I” in the table is used
throughout most of this study. “EOS II” on the other hand is used only in one
our simulations of gravitational collapse to a hybrid star above the twin star
mass range. In EOS II, the maximum mass of the NS branch is slightly lower
than $M_{ch}$, while the other characteristics are almost identical to those
of EOS I. For the remainder of this paper, we always refer to EOS I of Table 1
unless EOS II is explicitly specified. A plot of these EOSs is provided in
App. A.

The choice of piecewise polytropic parameters is not unique, but there is not
much room for changing these parameters while still meeting all aforementioned
conditions. Our goal here is not to explore all possible equations of state,
because this is not feasible. However, by adopting the resulting
phenomenological EOSs, we can perform our point-of-principle calculations.

In geometrized units111The conversion factor between cgs units and geometrized
units with $M_{\odot}$ is given by
$\rho_{\text{cgs}}=G^{-3}\,M_{\odot}^{-2}\,c^{6}\,\rho_{G}$ for density, and
by $R_{\text{cgs}}=G\,M_{\odot}\,c^{-2}\,R_{G}$ for length, where $\rho_{G}$
and $R_{G}$ stand for density and radius in geometrized units, respectively.
with $M_{\odot}=1$, the first coefficient of the EOS in Eq (2) is
$k_{1}=20.7$. The value of $k_{i}$ for $i>1$, is determined by continuity
condition at the boundary of each two neighboring segments, expressed as
$P_{i}(\rho_{0,i+1})=P_{i+1}(\rho_{0,i+1})$. This condition leads to

$k_{i+1}=k_{i}\,\rho_{0,i+1}^{\Gamma_{i}-\Gamma_{i+1}}\;.$ (5)

Table 1:  The two polytropic EOSs that are adopted in this work, presented as 8-branch piecewise polytropes with 16 free parameters ($\Gamma_{i}$, $\rho_{0,i}$, $k_{1}$). Setting $M_{\odot}=1$, the first coefficient is $k_{1}=20.7$, and the values of $k_{i}$ for the other segments are determined by continuity at $\rho_{0,i}$, as given by Eq. (5).  | EOS I | EOS II
---|---|---
$i$ | $\Gamma_{i}$ | $\log_{10}\rho_{0,i}[M_{\odot}^{-2}]$ | $\Gamma_{i}$ | $\log_{10}\rho_{0,i}[M_{\odot}^{-2}]$
1 | 1.5000 | — | 1.5000 | —
2 | 1.3350 | -9.8833 | 1.3350 | -9.8833
3 | 1.1386 | -7.4573 | 1.1286 | -7.4661
4 | 2.3544 | -3.8475 | 2.3144 | -3.8475
5 | 3.3458 | -3.3199 | 3.3858 | -3.3109
6 | 0.2576 | -3.0081 | 2.4957 | -3.0126
7 | 5.1878 | -2.7019 | 4.1878 | -2.7290
8 | 7.6102 | -2.5405 | 8.1102 | -2.5405

At sufficiently high densities, relevant for the most massive stellar
configurations, the sound speed calculated by Eq. (3) becomes superluminal.
Superluminality for ultra-dense matter with realistic EOSs is not uncommon.
For example, it arises in the Akmal-Pandharipande-Ravenhall (APR) EOS [114],
and efforts have been made to explain this behavior of $c_{s}$ in high
densities, e.g. in [115]. However, we checked that the sound speed never
becomes superluminal in our dynamical spacetime simulations.

During the process of adjusting the EOS, we find that in order to have heavy
HSs with a maximum mass of at least $2M_{\odot}$, the sound speed must be
large at densities corresponding to pure quark matter. This behavior aligns
closely with the findings of [116, 117, 118, 25] on the necessity of sound
speed becoming close to the speed of light in densities of a few times
$\rho_{s}$. Physically, this implies that pure quark matter at intermediate
densities is strongly coupled, and violates $c_{s}^{2}=1/3$ predicted for
asymptotically free quarks when $\rho_{0}\gtrsim 40\,\rho_{s}$ [119].

With the EOS available, we solve the TOV equations for a wide range of central
rest-mass densities $\rho_{c}$ between $10^{8}-10^{16}\textrm{g cm}^{-3}$. We
show the resulting mass-central density plot in Fig. 1. The plot displays the
total gravitational mass $M$, also known as the Arnowitt-Deser-Misner (ADM)
mass, and the total rest mass $M_{0}$ of static stars as functions of their
central rest-mass density, with red and blue curves, respectively. The
branches of stable ($dM/d\rho_{c}>0$) and unstable ($dM/d\rho_{c}<0$) WDs, NSs
and HSs are explicitly identified based on the turning point theorem [89].

Figure 2:  A 2-dimensional view (top panels) of the matter distribution inside
the stable NS configuration C (left) and the HS configuration E (right). For
the HS, the distinct quark core is easily distinguishable by its brighter
color, indicative of the sudden increase in density. The density profile of
each star is also displayed in the bottom panels as a function of the radial
coordinate $r$. Red color in the background stands for the pure hadronic
matter, blue represents the pure quark matter, and the narrow purple region in
the plot of the HS density profile denotes the zone containing a mixture of
quarks and hadrons.

Our initial data for the collapse simulations are based on stellar
configuration B in Fig. 1, with a rest mass $M_{0,s}$ very close to $M_{ch}$
but on the unstable WD branch to help accelerate the collapse. The mass of
this star, $1.44\,M_{\odot}$, falls right in the mass range of TSs and allows
us to perform our point-of-principle calculations to test for twin star
formation. The TOV solutions with this choice of mass are labeled on this plot
as $\\{$A, B, C, D, E$\\}$. Configurations A and B are a stable and an
unstable WD, respectively. Configurations C and E represent a stable NS and a
stable HS in the twin star mass regime, although, strictly speaking, TSs have
the same gravitational mass and not the same rest mass. The internal structure
of these two compact stars is illustrated in Fig. 2.

The largest ADM mass is attained in our EOSs is by a HS. Stars on the unstable
branch that separates stable NSs and stable HSs, such as configuration D, (see
inset in Fig. 1), which we will refer to as unstable TSs, are entirely in the
twin star mass range. Fig. 1 can help visualize the basic question in our
work. During the gravitational collapse of a stellar core/WD along a constant
rest-mass path (dashed green horizontal line in Fig. 1), as the remnant cools
it encounters first the stable NS branch and then the stable TS branch as it
contracts quasiadiabatically. This begs the following question: if the NS
branch is encountered first, then how can stable twin stars form?

We point out that during the evolution the remnant entropy changes due to
heating, therefore the actual evolutionary track that we describe above is not
on the plane of Fig. 1, because that corresponds to zero-entropy
configurations. However, core bounce occurs above the density of the NS
configuration B, and as remnant cools and contracts the evolutionary path
approaches the one shown by the horizontal dashed line in Fig. 1.

Figure 3:  Mass-radius relationship based on the EOS used in this work. Colors
and points indicate the same quantities and configurations as in Fig. 1.
Solid, dashed, and dotted curves stand for stable HSs, unstable TSs, and
stable NSs branches, respectively. The green shaded area represents the mass
range of TSs.

In addition to mass, the radius is another important macroscopic property of
compact objects. For any given EOS, solving the TOV equations also determines
the radius of each configuration. Consequently, every EOS has a unique mass-
radius (M-R) relation that can be probed through multimessenger observations,
and subsequently, impose constraints on the EOS of dense nuclear matter (see,
e.g., [120, 121, 19, 122, 123, 27] for reviews). The M-R relation based on our
EOS is displayed in Fig. 3. The three compact stars labeled in Fig. 1 are also
shown in this plot with the same labels. We show these stars on the rest-mass
vs radius curve, because $M_{0}$ is conserved during collapse.

##  III  Initial data and evolution Methods

In the $M-\rho_{c}$ plane of Fig. 1, the unstable WD configuration B will
collapse at constant rest mass, and following cooling, it must ultimately
migrate towards one of the two stable configurations with the same rest-mass,
i.e., at points C or E. During its evolution, the star moves through
configurations that cannot be accurately described by cold equilibrium models,
so a priori, it is unclear which of these stable states the star will end up
to. To determine the ultimate fate of the unstable WD, we perform multiple
simulations of gravitational collapse by perturbing configuration B. As the
initial configuration contracts, it heats adiabatically, and via shocks, if
the core bounces. Therefore, to end up on a cold degenerate configuration, the
star has to cool. We adopt a covariant local effective cooling approach, the
details of which we present in App. B. Our cooling model is that of [124], and
is characterized by a single parameter–the cooling timescale $\tau_{c}$.

Except for the code that builds our tabulated EOS in Python, and our TOV
solver, our computational framework is based on the Einstein Toolkit [125,
126] – Cactus [127] and Simfactory [128], employing Carpet for mesh refinement
[129]. Post-processing of the simulation data is performed using the kuibit
Python package [130].

###  III.1 Initial data

The initial data are provided by the Einstein Toolkit thorn RNSID which is
based on the RNS code [131, 132]. This code builds (rotating) isolated stellar
configurations assuming a zero-temperature EOS in either tabulated or
polytropic forms. As our EOS is represented by a piecewise polytrope 1, we
convert it to a highly sampled table before passing it to RNSID. For
simplicity, we ignore rotation in this work.

To vary the degree of violence of the collapse of our initial configuration we
consider a set of different initial perturbations, that involve non-zero
velocity and/or pressure perturbations. We also consider mass depletion
perturbations. We describe these below. After we impose the initial
perturbations we do not resolve the Einstein constraints, but we check in all
our cases that the constraint violations are always small during the
evolution.

####  III.1.1 Pressure perturbation

The pressure perturbations are modeled as

$P\rightarrow(1+\xi_{p})P=P+\delta P\;,$ (6)

for every point in the configuration. Here, the value of $\xi_{p}={\delta
P}/{P}$ controls the strength of the perturbation. The pressure perturbation
could be either negative or positive. A negative $\delta P$ corresponds to
pressure depletion and consequently, speeds up the collapse. On the contrary,
a positive pressure perturbation corresponds to heating the stellar interior
that must be radiated away to maintain the equilibrium. We use the latter to
test our implementation of effective cooling in App. C.

####  III.1.2 Velocity perturbation

Our model for velocity perturbation is given by the following two-parameter
function

$\vec{v}\rightarrow\xi_{v}\,{\left(\frac{r}{R}\right)}^{\kappa}\;\hat{r}\;,$
(7)

where $\vec{v}$ is the coordinate 3-velocity, $\hat{r}$ the radial unit
vector, $r$ is the coordinate radius, $R$ the coordinate radius of the star,
and $\xi_{v}$ determines the amplitude of the velocity perturbation, and could
be either negative (collapsing star) or positive (expanding star). The
exponent $\kappa$ controls the radial profile of the perturbation. The case
$\kappa=0$ is a particular model with a velocity perturbation independent of
position. The choice of a power-law perturbation model is made for simplicity
and is not unique. We choose $\kappa\neq 0$ so the perturbation goes smoothly
to zero close to the center of the star. This is because for
$\xi_{v}=\mathcal{O}(-0.1)$ and $\kappa=0$ the star bounces almost immediately
and explodes.

####  III.1.3 Mass perturbation

To simulate a configuration that experiences a small amount of mass loss,
e.g., lost due to winds, the initial data are perturbed such that a low-
density shell of the stellar structure is set to the tenuous atmospheric
density we maintain in our simulations. This perturbation essentially changes
the total rest mass as

$M_{0}\rightarrow(1+\xi_{m})M_{0}=M_{0}+\delta M_{0}\;.$ (8)

Here, $\xi_{m}={\delta M_{0}}/{M_{0}}$ specifies the fraction of rest mass
removed or added from/to the original star, and hence, takes a negative value
in the case of mass loss. In practice, this is done by defining a threshold
rest-mass density below which the rest-mass density of the stars takes the
value of the tenuous atmospheric density until the desired amount of mass
$\delta M_{0}$ is depleted.

The parameters characterizing all simulations performed in this work are
summarized in Table 2.

Table 2: The set of simulations conducted in this study using EOS I is detailed below. For each case, the corresponding configuration in the $M-\rho_{c}$ plot for the initial data are listed, along with the perturbation parameters. The last column describes how cooling is performed during the simulation, if enabled. Further details are provided in Sec. IV. Configuration F represents a stable HS with $M_{0}=1.49M_{\odot}$, and is discussed in Sec. IV.3. The simulations based on the last two rows are presented in App. D.  Description | Initial TOV solution | $\xi_{p}$ | $\xi_{m}$ | $\xi_{v}$ | $\kappa$ | Cooling
---|---|---|---|---|---|---
Unstable WD collapse | B | -0.01 | 0.00 | 0.00 | 0.0 | SI
Unstable WD collapse | B | -0.01 | 0.00 | 0.00 | 0.0 | SII
Unstable WD collapse — perturbed pressure | B | -0.90 | 0.00 | 0.00 | 0.0 | SI
Unstable WD collapse — perturbed pressure | B | -0.90 | 0.00 | 0.00 | 0.0 | SII
Unstable WD collapse — perturbed velocity | B | -0.01 | 0.00 | -0.10 | 0.0 | —
Unstable WD collapse — perturbed velocity | B | -0.01 | 0.00 | -0.10 | 0.5 | SII
Unstable WD collapse — perturbed velocity | B | -0.01 | 0.00 | -0.10 | 1.0 | SII
Unstable WD collapse — perturbed velocity | B | -0.01 | 0.00 | -0.10 | 5.0 | SII
Unstable WD collapse — highly perturbed | B | -0.90 | 0.00 | -0.10 | 0.5 | SII
Unstable WD collapse — highly perturbed | B | -0.90 | 0.00 | -0.10 | 1.0 | SII
Unstable WD collapse — highly perturbed | B | -0.90 | 0.00 | -0.10 | 5.0 | SII
Unstable WD collapse — highly perturbed | B | -0.90 | 0.00 | -0.30 | 1.0 | SII
Unstable WD collapse — highly perturbed | B | -0.90 | 0.00 | -0.40 | 1.0 | SII
Stable massive HS | F | 0.00 | 0.00 | 0.00 | 0.0 | —
Massive HS — mass loss | F | 0.00 | -0.03 | 0.00 | 0.0 | —
Stable NS | C | 0.00 | 0.00 | 0.00 | 0.0 | —
Stable HS | E | 0.00 | 0.00 | 0.00 | 0.0 | —

###  III.2 Evolution

The spacetime initial data are evolved in the Baumgarte-Shapiro-Shibata-
Nakamura (BSSN) formulation [133, 134] as implemented in the public Lean code
[135]. The spacetime gauge choice adopted here includes the “1+log” and the
“$\Gamma$ – driver” conditions [136, 137]. The hydrodynamic initial data are
evolved with the publicly available code IllinoisGRMHD [138, 139]. We have
implemented in IllinoisGRMHD the cooling of [124] that we describe in App. B .
We also test our implementation of cooling in App. C. The EOS for the
evolution is hybrid with a cold and thermal component

$P=P_{\rm cold}+P_{\rm th},$ (9)

where the cold pressure $P_{\rm cold}$ is given by EOS I or II, and $P_{\rm
th}=(\Gamma_{\rm th}-1)\epsilon_{\rm th}$, with $\epsilon_{\rm th}$ being the
thermal energy density. We employ $\Gamma_{\rm th}=2$.

In the WD collapse simulations, since the radius of the initial configuration
(WD) is about 100 times larger than that of the final configuration (HS/NS),
it is computationally inefficient to perform the simulation with the highest
resolution necessary to resolve a HS/NS from the start. Thus, we employ
adaptive mesh refinement. The initial grid structure consists of 2 refinement
levels, with the finer resolution set to 100 points across the WD. However, as
the radius decreases over time (and the central density increases), the
initial resolution becomes inadequate, and we add progressively higher
resolution refinement levels. Based on an estimate of the final central
density as $\rho_{c,f}\sim 10^{-3}\,M_{\odot}^{-2}$ using the equilibrium cold
configurations (see Fig. 1) we construct 6 refinement levels that are
initially inactive. As the central rest-mass density rises, these refinement
levels are subsequently activated when the density reaches certain values.
When all the refinement levels are activated, the highest grid resolution is
$166\,\textrm{m}$. The density range between the initial central density,
$\rho_{c,i}$, and $\rho_{c,f}$ is divided into 6 equidistant logarithmic
segments. Consequently, there are 7 values of density, separating these 6
segments, that we denote $\\{\rho_{1},\rho_{2},\ldots,\rho_{7}\\}$, where
$\rho_{1}=\rho_{c,i}$ and $\rho_{7}=\rho_{c,f}$. Every time the maximum
density of the star exceeds
$\rho_{k}={\left(\frac{\rho_{c,f}}{\rho_{c,i}}\right)}^{1/6}\times\rho_{k-1}$
for $k=\\{2,3,\ldots,7\\}$, a new refinement level is activated. In the
collapse simulations, the outer boundary is set to $1.7R_{\rm WD}\approx
498M_{\rm WD}$ ($R_{\rm WD}=1058.55\,\textrm{km}$ is the initial WD radius and
$M_{\rm WD}$ its mass.). The half-side length of each new refinement level is
half that of the preceding level and its resolution twice that of the
preceding level. For the other simulations in this work (e.g. starting with a
more massive HS), the grid has the same structure as described for the stable
HS in App. D.

As the star collapses, it heats adiabatically, eventually resulting in a core
bounce once the EOS has stiffened and significant heat has been generated.
Without taking cooling into account, the collapse would stall and the star
would not continue to contract to reach either of the equilibrium
configurations on the TOV sequence. This excess heat is naturally radiated
away in the form of neutrinos.

In the absence of a physically accurate 6+1 dimensional neutrino code, we
crudely model cooling by locally removing any excess heat, as described in
App. B. The only free parameter that needs to be specified in this model is
the cooling timescale $\tau_{c}$. The value of $\tau_{c}$ has to be chosen
such that we respect the hierarchy of timescales in our problem, while also
considering the duration over which the simulations can be completed. If
$\tau_{c}$ is much larger than the dynamical timescale $t_{dyn}$, the
computations will take a very long time to complete. However, $\tau_{c}$ is
much smaller than $t_{dyn}$, this would result in rapid cooling, causing
significant perturbations instead of a smoother transition to the zero-
temperature remnant. As a star collapses, the density increases. Thus, the
dynamical timescale of the star, defined as
$t_{dyn}\equiv\frac{1}{\sqrt{\rho_{c}(t)}}$, is not constant over time.

We follow two different strategies to activate cooling in our simulations in
order to ensure that our final results are invariant with the cooling
strategy. In the first strategy (SI), the collapse begins without cooling, and
subsequently comes to a halt due to core bounce and heat generation. When the
star settles in this state, we activate cooling with $\tau_{c}=3\,t_{dyn}$,
where $t_{dyn}$ is determined based on the central density of the settled
configuration. In the second strategy (SII), cooling is activated from the
beginning with $\tau_{c}=3\,t_{dyn}$. However, as the collapse continues,
$\tau_{c}$ is updated to track the changing dynamical timescale as the maximum
rest-mass density increases. We update $\tau_{c}$ every time a new refinement
level is added to the configuration with the same relation
$\tau_{c}=3\,t_{dyn}$, where $t_{dyn}$ is determined based on the values of
$\\{\rho_{1},\rho_{2},\ldots,\rho_{7}\\}$ discussed above.

##  IV Results

In this section, we present the results of our simulations. We begin with the
simulations of WD collapse, both without and with strong initial
perturbations. Subsequently, we present the outcomes of the simulation
involving a HS with mass loss. Furthermore, we perform a simulation to explore
the formation of a stable HS with a mass exceeding the mass range of TSs,
arising from the collapse of a WD. We present how well the constraints are
satisfied during the evolution in App. E. We also perform two simulations
involving the stable NS and HS shown in Fig. 2 to demonstrate the code can
reliably evolve stable stars. The results of these simulations are discussed
in App. D.

###  IV.1 White dwarf collapse

A potential pathway to forming TSs is the gravitational collapse of the iron
core of a massive star or the accretion-induced collapse of a WD. An unstable
WD near the Chandrasekhar limit is an acceptable model for either scenario.
The initial central density and radius of the unstable WD at configuration B
are $\rho_{c}=2.14\times 10^{10}\textrm{g}\,\textrm{cm}^{-3}$ and $R_{\rm
WD}=1058.55\,\textrm{km}$. Following cooling, this configuration is expected
to collapse and settle into configuration C or E, since these are the only
stable configurations that preserve the total rest mass. To speed up the
collapse, we initially deplete $1\%$ of the pressure at every point inside the
WD, that is, $\xi_{p}=-0.01$.

Figure 4:  Snapshots of matter distribution within the star at four different
times during the collapse in the xy plane. The leftmost panel shows the
unstable WD configuration B, and the rightmost panel shows the stable NS
configuration C. The matter left as a low-density atmosphere around the stable
NS gradually falls onto the dense core during the final stages of the collapse
but its total mass is negligible. Thus, it cannot raise the maximum density,
as evidenced in Fig. 5.

Fig. 4 shows the outcome of the simulation adopting the SII cooling strategy
with four snapshots of the matter distribution. The dense core depicted in the
final snapshot of this figure shows the ultimate stable compact star. The mass
in the low-density region surrounding the core is negligible compared to that
of the central compact object. Fig. 5 presents the time evolution of the
maximum rest-mass density inside the star. The horizontal orange and blue
dashed lines display the central rest-mass density of the stable NS and the
stable HS with the same rest mass as the initial WD, i.e. configurations C and
E, respectively. The black dashed (green) curve shows the evolution of the
maximum rest-mass density in the SI (SII) cooling strategy.

In both scenarios, the initial configuration slowly contracts. After $t\sim
6\,t_{dyn}$ (here and in Fig. 5 $t_{dyn}$ corresponds to the dynamical
timescale of the initial WD), the collapse accelerates. Around $t\sim
6.70\,t_{dyn}$, a bounce in the density evolution occurs, halting the rapid
collapse. We note that the evolution until core bounce is insensitive to
whether cooling is active or not. After the bounce, the evolution of the
maximum rest-mass density depends on the cooling strategy. In the case where
cooling is initially inactive, the star settles at a density below both stable
NS and stable HS. Note that this configuration does not lie on the TOV
sequence with the underlying cold EOS, because it no longer has zero entropy.
As soon as cooling is activated, this intermediate configuration undergoes
additional collapse, and asymptotes to the stable NS configuration C (orange
dashed line). In the SII case, the star continues its collapse after the
bounce, since cooling already had reduced the thermal pressure. However, the
asymptotic configuration, when most of the heat has been removed, is again the
stable NS configuration C. Thus, the outcome of the collapse is independent of
the cooling strategy adopted here and the final cold configuration has a
maximum density that is below the quark-hadron phase transition region.
Therefore, the final configuration is a NS.

Figure 5:  Evolution of the maximum density for an unstable WD under 1%
pressure depletion. Here, $t_{\rm dyn}=26.46\,\text{ms}$ is the dynamical
timescale of the initial WD. The orange and blue dashed lines denote the rest-
mass density of the stable compact configurations C and E of Fig. 1,
respectively. The black dashed (green) curve illustrates the evolutionary path
in the SI (SII) cooling strategy. The inset zooms in on the bounce.

We note that in the part of the plot in Fig. 5 where the green and black
curves plateau at the stable NS threshold, every unit of $t/t_{dyn}$
corresponds to $\sim 170$ dynamical timescales of that stable NS. Hence, the
simulations were long enough to ensure that the evolution has reached the
stable cold configuration.

We also conducted a simulation with $\tau_{c}=1.5\,t_{dyn}$ as opposed to
$\tau_{c}=3\,t_{dyn}$, to check if the cooling timescale can affect the
results. We confirmed that faster cooling only accelerates the collapse, while
leaving the final outcome unaffected. Thus, we conclude that the final product
of the collapse of an unstable WD is a stable NS with the same total rest
mass; this initial data does not result in a stable HS, regardless of the
cooling approaches tested.

Figure 6:  Same as Fig. 5, but for the case with 90 percent of pressure
depletion in the initial data. The line styles and colors are coded as in Fig.
5, and similarly, $t_{\rm dyn}$ is the dynamical timescale of the initial WD

.

###  IV.2 Strongly perturbed initial data

Here we explore more violent perturbations to investigate whether the final
remnant can be a TS if the strong bounce can be avoided or happen at a higher
density past the hadron-to-quark phase transition. To cause the bounce to
occur at higher densities, one can potentially think of more extreme initial
conditions for the unstable WD.

Figure 7:  Maximum density vs time for the simulations of WD implosion with
initial perturbations in both velocity and pressure. Time is scaled with the
dynamical timescale of the initial WD. The colors maroon, green, and magenta
designate the three values of $\kappa=\\{0.5,1,5\\}$, respectively. The
velocity perturbation is fixed at $\xi_{v}=-0.1$. Dotted lines correspond to
$\xi_{p}=-0.01$, while solid lines to $\xi_{p}=-0.90$. All six simulations
employ the SII cooling strategy. The dashed lines designate the central
density of the stable NS/HS as in Fig. 5.

####  IV.2.1 Large pressure depletion

One of the parameters that could potentially facilitate the collapse is the
strength of the pressure depletion characterized by $\xi_{p}$. A very large
pressure depletion might result in a more rapid collapse. We conduct a set of
simulations with SI and SII cooling beginning with configuration B but
$\xi_{p}=-0.90$. While this level of pressure depletion is unlikely in nature,
we use it as an extreme case to determine if collapse alone can form a TS.
Since the WD is basically a Newtonian object, such a large pressure depletion
has a negligible impact both on the total energy budget of the star and on the
constraints.

The results of these simulations are displayed in Fig. 6. By comparing Fig. 6
and Fig. 5, it can be inferred that an exceptionally large pressure depletion
accelerates the collapse substantially. In this case, the bounce occurs at
$t\sim 2.42\,t_{dyn}$, almost 3 times earlier. Despite that, the core still
bounces at approximately the same density and halts the collapse. Following
cooling, the remnant eventually settles into a stable NS configuration.
Similar to the previous set of simulations, the end product is independent of
whether the SI or SII cooling strategy is adopted.

####  IV.2.2 Velocity perturbation

Given that the binding energy of the TS is larger than that of the NS with the
same rest-mass, it is possible that we must inject energy into the system to
overcome the bounce. One way to induce a more violent implosion is by
introducing an inward velocity perturbation. This perturbation is modeled with
Eq. (7), which consists of two free parameters $\xi_{v}$ and $\kappa$ to be
specified. We fix $\xi_{v}=-0.1$, which corresponds to an inward velocity at
the surface at 10% of the speed of light, and perform simulations with three
different values of $\kappa$, namely $\kappa=\\{0.5,1,5\\}$. Fig. 7 shows the
result of these simulations in the evolution plot of the maximum density,
represented by maroon, green, and magenta colors, respectively.

To explore a more diverse range of initial conditions, each of these
simulations is performed in two cases: $\xi_{p}=-0.01$ (dotted curves) and
$\xi_{p}=-0.90$ (solid curves). In all these simulations, cooling is activated
from the beginning, that is, we employ the SII cooling strategy. All these
cases yield the same final result as the previous simulation of the collapsing
WD; the original bounce near the last stages of the collapse is still
unavoidable, and the remaining core is still a stable NS, regardless of the
value of $\kappa$. Different values of $\kappa$ can only change the timescale
of the collapse. Among these three choices, $\kappa=0.5$ and $\kappa=5$ result
in the fastest and the slowest processes, respectively. Moreover, the
combination of the velocity perturbation with a larger pressure depletion in
the initial data leads to a faster implosion in all cases.

As the last set of simulations for the collapse of an unstable WD, we also
explore velocity perturbations with larger values of $\xi_{v}$. Fig. 8
presents these simulations with $\xi_{v}=-0.3$ (dotted curve) and
$\xi_{v}=-0.4$ (solid curve). Given that all choices of $\kappa$ led to the
same ultimate result, we only run these simulations for the case of
$\kappa=1$. Both simulations have 90 percent of pressure depletion,
$\xi_{p}=-0.90$, and are conducted in the SII. These represent the most
extreme initial conditions examined in this study. However, the result is once
again the stable NS at configuration C. The bounce still manifests after the
exponential collapse, although the process is even faster than in the previous
cases. The general outcome of these simulations suggests that even extreme
pressure and velocity perturbations cannot change the eventual fate of the
collapsing WD. The outcome is always a NS.

Figure 8:  The evolution of the maximum density of the collapsing WD with the
most extreme initial conditions in pressure and velocity. The dotted curve
shows the simulation with $\xi_{v}=-0.3$, and the solid curve represents the
simulation with $\xi_{v}=-0.4$. Both the cases are performed in SII, with
$\kappa=1$ and $\xi_{p}=-0.90$. Similar to Fig. 5, the dashed lines show the
maximum density of the stable NS/HS, and $t_{\rm dyn}$ is the dynamical
timescale of the initial WD.

###  IV.3 Mass loss

In this section, we investigate the scenario where the NS maximum mass is
somewhat lower than the Chandrasekhar limit so that a stellar core or a WD at
the Chandrasekhar limit collapses to the stable hybrid star branch above the
twin star mass regime. If this initial hybrid star experiences some mass loss,
e.g., due to strong winds, it could settle into a TS instead of an NS.

To model this scenario, we take the hybrid star configuration F in Fig. 9, and
deplete a small fraction of its mass as described in Eq. (8). The initial rest
mass of this configuration is $M_{0}=1.49\,M_{\odot}$, and exceeds $M_{ch}$.
The initial central rest-mass density is $\rho_{c}=1.63\times
10^{15}\,\textrm{g}\,\textrm{cm}^{-3}$ and the radius is
$R=11.02\,\textrm{km}$. By setting $\xi_{m}=-0.033$, the rest mass of the star
becomes $M_{0}\approx M_{0,s}$, aligning with the rest masses of stable stars
A, C, and E of Fig. 1. This initial perturbation is fairly strong compared to
a slow loss of mass through winds, but it allows us to get a glimpse into
whether this pathway to forming TSs is viable. We checked the constraint
violations after the mass loss remain small throughout the simulations we
performed (see discussion in App. E).

Figure 9:  A zoomed-in view of the NS/HS part of the Fig. 1. Configuration F
is a stable HS with a mass higher than the other three configurations
indicated. The green shaded area roughly denotes the mass range for TSs. The
red (blue) curve corresponds the ADM (rest) mass. Figure 10:  Time evolution
of the maximum density of configuration F without any perturbation (violet)
and with a 3.3 percent mass loss (black). The stable NS and HS configurations
C and E are shown in orange and blue, respectively. For each stable star, the
solid line denotes the actual simulation, while the dashed line represents the
initial central density. Furthermore, the shaded area indicates the density
range of the phase transition. In this plot, $t_{\rm dyn}$ is the dynamical
timescale of the initial HS. Figure 11:  Rest mass as a function of central
density for a sequence of TOV stars with EOS I (dotted) and EOS II (solid).
The dashed line represents the Chandrasekhar mass. Notably, there is no TS in
EOS II with a rest mass close to $M_{ch}$. In the density range relevant to
WDs, the two curves exhibit almost complete overlap. Thus, initial unstable WD
configuration (B) shares similar characteristics in both EOSs.

Fig. 10 illustrates the results. Four simulations are shown in this figure;
the stable configurations C, E, and F (without perturbation) are displayed by
the orange, blue, and violet solid lines, respectively. The black solid curve
represents the evolution of the maximum density of a star beginning at
configuration F but perturbed by $\xi_{m}=-0.033$. Interestingly, this plot
demonstrates that the last scenario successfully ends up with the formation of
a stable HS in the twin star regime. As a result of this large initial
perturbation, two types of oscillations can be seen around the maximum density
of the stable HS shown in blue. The first type exhibits oscillations with a
period of $\sim t_{dyn}$, which gradually disappears after a few dynamical
timescales. The second type of oscillation has a larger domain with $\sim\pm
10\,\%$ of deviations around the mean value of the central rest-mass density
of the TS E. The latter type of oscillation has been seen before in
simulations with EOSs that exhibit a strong phase transition, and is most
likely due to a significant portion of the mass of the star transitioning in
and out of the density range in the hadron-to-quark phase transition, where
the stiffness of the EOS changes rapidly leading to collapse (EOS softens) and
bounce (EOS stiffens). We expect that a gradual loss of rest-mass would make
such oscillations disappear, but the main result of the final configuration
being a TS with the same rest mass should hold. The important finding here is
that the maximum central density remains above the phase transition region as
shown in Fig. 10, and its mean is that of the TS with the same rest mass.

###  IV.4 Hybrid star formation

In the previous subsection, we demonstrated that TSs can form via mass loss
from a more massive HS. However, this process is plausible only if a massive
HS can form in the first place. In this subsection, we explore the formation
of a HS with a higher mass than TSs through gravitational collapse.

The simulation we perform here is the only one in this work based on EOS II
listed in Table 1. Fig. 11 compares the $M_{0}-\rho_{c}$ curves resulting from
these two EOSs. The TOV sequence resulting from EOS II has an NS peak lower
than $M_{ch}$. Therefore, there are no TSs with rest mass $M_{0,s}$.

Starting from an unstable WD with rest mass $M_{0,s}$ as initial data, there
should be only one stable equilibrium, cold configuration at higher density.
This equilibrium should correspond to a HS with a central density of
$\rho_{c}=1.75\times 10^{15}\,\textrm{g}\,\textrm{cm}^{-3}$ and a radius
$R=10.35\,\textrm{km}$. The results of this simulation are displayed in Fig.
12, where we show the evolution of the maximum rest-mass density. The blue
line denotes the described stable HS, and the shaded area shows the range of
densities in the phase transition. As expected, the remnant of the collapse is
this HS, as evidenced by the convergence of the maximum density to the central
density of this star.

Figure 12:  The evolution of the maximum density of a star beginning as an
unstable WD with a rest mass $M_{0}=M_{0,s}$. The blue line represents the
central density of the stable HS with this rest mass, and the phase transition
range of densities is indicated by the shaded gray area. Here $t_{\rm dyn}$
represents dynamical timescale of the initial WD. The inset zooms in on the
final stages of collapse, indicating the convergence of maximum density to the
central density of the stable HS.

The bounce observed in the previous simulations persists in this simulation as
well, due to the stiff EOS. SII cooling is employed in this simulation,
starting at $\tau_{c}=3\,t_{dyn}$. When the maximum density approaches the
phase transition segment of the EOS, the evolution of the maximum density
becomes very slow. Since we have already confirmed that the cooling timescale
does not affect the outcomes, we change the cooling timescale to
$\tau_{c}=1.5\,t_{dyn}$ at $t/t_{dyn}\sim 6.7$ to accelerate the transition.
As the maximum density enters the segment of the EOS that softens, the
configuration undergoes fast collapse to the HS formation. The only
perturbation introduced in this simulation is $\xi_{p}=-0.01$, applied solely
to speed up the initial collapse.

While not surprising, our simulations here demonstrate explicitly that forming
HSs with masses higher than TSs is viable.

##  V Conclusions

The exact form of the cold EOS beyond the nuclear saturation density is still
quite uncertain. Cores of NSs are a common place in nature where neutron-rich
matter can be found in this regime. The properties of compact stars are
intricately linked to the form of the EOS at these densities. If a strong
hadron-to-quark phase transition occurs at these compact star densities, the
EOS can potentially give rise to a third branch of stable compact objects that
consist of a quark core surrounded by a hadronic shell. These configurations
are known as hybrid hadron-quark stars or hybrid stars. Associated with a
third family of compact stars is a range of masses where for each NS there is
a corresponding HS with the same mass but a smaller radius. These stars are
referred to as twin stars.

In this work, we explored the formation of TSs by conducting hydrodynamical
simulations in full general relativity. We constructed piecewise polytropic
EOSs that approximate realistic EOSs with a hadron-to-quark phase transition,
and we extended them down to the EOS of a white dwarf. Using our EOSs, we
performed multiple simulations of the collapse of an unstable WD under
different initial conditions with varying degrees of extremality of initial
perturbations in pressure and/or velocity. Given that total rest mass is
conserved, the ultimate product of the collapse could either be a stable NS or
a stable HS. The general finding from all our simulations is that the unstable
WD always collapses into a stable NS, regardless of how extreme the initial
conditions are. This overall result remains the same even under extreme
initial perturbations.

If TSs exist, there should be at least one pathway for them to form. Following
the standard theory of NS formation, a natural path could be the formation of
a more massive HS that subsequently loses a small amount of mass in the form
of winds. Another potential avenue for mass loss would be a “grazing”
collision of a massive HS with a black hole where the massive HS undergoes an
episode of mass loss and then flies away. Such a scenario might take place in
a dense stellar cluster. Our simulations show that forming HSs heavier than in
the twin star mass regime is possible, and that forming a TS through a heavier
HS experiencing mass loss is a viable path. However, even this scenario
appears that it would require some fine-tuning for TSs to form, which would
place strong limits on the abundance of TSs. The fine tuning for the formation
of twin stars results because of the limited range of masses over which they
are predicted to reside by existing EOSs with hadron-to-quark phase
transitions, and hence they should be fairly rare objects. This conclusion and
the challenge in producing TSs as demonstrated by our simulations suggest that
if a HS star was involved in GW170817, then it likely was not a twin star.

We point out that our work is idealized in several ways and has a number of
caveats. First, we do not treat realistic neutrino effects or magnetic fields,
and the perturbations studied here represent a simplified version of processes
that may occur inside the dense core of massive stars during the last stages
of their evolution. In particular, scenarios involving mass loss due to winds
and mass gain resulting from the infall of matter from outer shells of the
initial progenitor are expected to happen in nature. The results of our
simulations indicate that these perturbations play a notable role in the
potential formation of stable TSs, assuming that they exist in nature. Thus,
it is important to treat these self-consistently. Nevertheless, our conclusion
that NSs tend to be the preferred outcome in the twin star mass range should
hold, not only because this is what our simulations demonstrate, but also
because of the very small range of masses over which twin stars are predicted
to exist by current EOSs, which in turn requires a delicate balance of mass
loss and mass gain.

To develop a more comprehensive model for the formation of TSs, further cases
should be studied with a wider range of parameters. Our simulations were
performed using static configurations. Rotation is another aspect that can be
added with different models to future simulations to study the same problem.
Additionally, these simulations were conducted considering only one type of
EOS. Neutrino effects, nuclear reaction networks, and a more realistic
progenitor are other ways to increase the realism of these simulations.

Another potential pathway to forming HSs involves mergers of binary NSs or WD-
NS [140, 124]. Exploring all these scenarios and including more realistic
physics is left for future work. Most work in the area of TSs revolves around
either observational constraints on or properties of the EOS. Our work here
points out that there is another important problem surrounding TSs: their
formation and abundance in nature.

###### Acknowledgements.

We thank Thomas Baumgarte for useful discussions. We are grateful to Zachariah
Etienne and Leo Werneck for releasing their work on the IllinoisGRMHD code
which supports piecewise polytropic EOS, and to all the Einstein Toolkit
developers and maintainers. We also thank Pedro Espino for sharing his EOS
tables, and D. Alvarez-Castillo, D. Blanschke and A. Sedrakian for giving us
permission to use the ACS-II equations of state developed in [34]. Our Python
TOV solver for tabulated equations of state was developed starting from a code
originally created by Nikolaos Stergioulas. Simulations were performed on the
Puma cluster at the University of Arizona, and on ACCESS resources under
allocation PHY190020. This work was supported in part by NSF grant
PHY-2145421, and NASA grant 80NSSC24K0771 to the University of Arizona. G. B.
was in part supported by the Frontera Fellowship by the Texas Advanced
Computing Center (TACC). Frontera is funded by NSF grant No. OAC-1818253.

##  Appendix A Extended EOS information

The EOSs used in this work consist of 8 segments and must satisfy the list of
conditions in Sec. II. Empirically we have found that each $\Gamma_{i}$
controls specific characteristics of the TOV sequence shown in Fig. 1.
$\Gamma_{1}$ primarily affects the maximum WD mass, $\Gamma_{2}$ alters the
density and mass of the WD peak, $\Gamma_{i}$ for $i=\\{3,4,5\\}$ control the
maximum NS mass, and $\Gamma_{6}$ changes the twin star mass range. Finally,
$\Gamma_{i}$ for $i=\\{7,8\\}$ affects the maximum HS mass. However, it is
important to emphasize that each of the free parameters affects more than just
one characteristic.

The two EOSs listed in Table 1 exhibit only small differences, which are
mostly noticeable at higher densities. Fig. 13 illustrates the pressure as a
function of energy density $\epsilon$ for our EOSs. The shaded red region
designates the density range where matter exists in the hadronic phase, and
the blue region shows the density range where it exists purely in quark form.
The narrower purple range, where the two regimes overlap ($i=6$ branch of the
EOS), is where the phase transition occurs. In this regime, the original
polytropic exponent was $\Gamma_{6}=0$, corresponding to a 1st-order phase
transition, but we modified it to $\Gamma_{6}=0.2576$ to avoid a zero sound
speed. This implies that, instead of a sharp boundary from the hadron to the
quark phase, there exists a mixture of hadrons and quarks in this narrow
region of our EOS. This modification makes the sound speed of this region
$\simeq 15\,\%$ of the speed of light, which allows for stable numerical
integration of the general relativistic hydrodynamic equations.

Figure 13: The relationship between the pressure and the energy density based
on the EOSs constructed in Table 1. EOS I is represented by the solid line,
while the dotted line denotes EOS II. The shaded red, blue, and purple areas
correspond to the pure hadronic, pure quark, and the mixed quark-hadron phases
for EOS I, respectively. The inset zooms in on the phase transition region.

##  Appendix B  Cooling formalism

In this appendix, we summarize the cooling formalism adopted in our work.

In the standard 3+1 decomposition formalism to numerically solve the Einstein
equation, the spacetime metric $g_{\mu\nu}$ is written as follows by [141]

$ds^{2}=-\alpha^{2}dt^{2}+\gamma_{ij}(dx^{i}+\beta^{i}dt)(dx^{j}+\beta^{j}dt)\;,$
(10)

where $\alpha$, $\beta^{i}$, and $\gamma_{ij}$ are the lapse function, shift
vector, and spatial metric ($i$ and $j$ are spatial indices). For a perfect
fluid with 4-velocity $u^{\mu}$, the stress-energy tensor is written as

$T^{\mu\nu}=\rho_{0}hu^{\mu}u^{\nu}+Pg^{\mu\nu}\;.$ (11)

Here, the specific enthalpy $h$ is related to the specific internal energy
$e=\frac{\epsilon}{\rho_{0}}-1$ through

$h=1+e+\frac{P}{\rho_{0}}\;.$ (12)

The evolution equations for general relativistic hydrodynamics (GRHD), are
given by [142]

$\begin{gathered}\partial_{t}\rho_{\ast}+\partial_{j}(\rho_{\ast}v^{j})=0\;,\\\\[10.00002pt]
\partial_{t}\tilde{\tau}+\partial_{i}(\alpha^{2}\sqrt{\gamma}\;T^{0i}-\rho_{\ast}v^{i})=s\;,\\\\[10.00002pt]
\partial_{t}\tilde{S}_{i}+\partial_{j}(\alpha\sqrt{\gamma}\;T^{j}_{\;\;i})=\frac{1}{2}\alpha\sqrt{\gamma}\;T^{\mu\nu}\;\partial_{i}g_{\mu\nu}\;,\\\\[10.00002pt]
\end{gathered}$ (13)

In these equations, $v^{i}=u^{i}/u^{0}$ is the coordinate three velocity, and
the conservative variables are defined as

$\begin{gathered}\rho_{\ast}=-\sqrt{\gamma}\;\rho_{0}n_{\mu}u^{\mu}\;,\\\\[10.00002pt]
\tilde{S}_{i}=-\sqrt{\gamma}\;T_{\mu\nu}n^{\mu}\gamma^{\nu}_{\;\;i}\;,\\\\[10.00002pt]
\tilde{\tau}=\sqrt{\gamma}\;T_{\mu\nu}n^{\mu}n^{\nu}\rho_{\ast}\;,\end{gathered}$
(14)

where $n^{\mu}=\frac{1}{\alpha}(1,-\beta^{i})$ is the timelike unit vector
normal to $t=$ constant slices.

In Eqs. (LABEL:GRMHD), the first equation corresponds to the baryon number
conservation. The time component of energy-momentum conservation
$\nabla_{\nu}T_{\mu}^{\;\;\nu}$ is expressed in the second equation, where
$s=-\alpha\sqrt{\gamma}\;T^{\mu\nu}\nabla_{\nu}n_{\mu}$ is a source term. The
spatial components of the energy-momentum conservation are encoded in the
third equation.

The set of Eqs. (LABEL:GRMHD) also requires an EOS to close. The general form
of the EOS that is assumed to solve this set of equations is given by Eq. (9)
[143], where $P_{cold}$ stands for the pressure of the zero-temperature
matter, and has the polytropic form introduced in Sec. II. Any additional
pressure due to heating, e.g., by shock heating, is encapsulated in $P_{th}$.
This term is considered as a $\Gamma$-law,

$P_{th}=(\Gamma_{th}-1)\rho_{0}e_{th}\;,$ (15)

where $e_{th}=e-e_{cold}$, and $e_{cold}$ is the specific internal energy
associated with the cold pressure, which can be found from the first law of
thermodynamics as

$e_{cold}(\rho_{0})=-\int
P_{cold}(\rho_{0})d\left(\frac{1}{\rho_{0}}\right)\;.$ (16)

For the specific case of Eq. (2), this becomes

$e_{cold}=\frac{k_{i}}{\Gamma_{i}-1}\rho_{0}^{\Gamma_{i}-1}+a_{i}\;,$ (17)

where $a_{i}$ is a constant of integration. Considering $\epsilon_{cold}$ as
the energy density when there is no heat, this constant is determined as [96]

$a_{i}=\frac{\epsilon_{cold}(\rho_{0,i})}{\rho_{0,i}}-1-\frac{k_{i}}{\Gamma_{i}-1}\rho_{0,i}^{\Gamma_{i-1}}\;.$
(18)

In the presence of radiation, the conservation of energy-momentum must be
modified to incorporate the radiation stress-energy tensor $R^{\mu\nu}$,

$\nabla_{\mu}(T^{\mu\nu}+R^{\mu\nu})=0\;.$ (19)

The dynamics of the radiation stress-energy tensor is described by
$\nabla_{\mu}R^{\mu\nu}=-G^{\nu}$, where $G^{\mu}$ is the radiation 4-force
density. Assuming local and isotropic (in a frame comoving with the fluid)
cooling, the resulting energy and momentum equations are then given by [124]

$\begin{gathered}\partial_{t}\tilde{\tau}+\partial_{i}(\alpha^{2}\sqrt{\gamma}\;T^{0i}-\rho_{\ast}v^{i})=s-\alpha^{2}\sqrt{\gamma}\;u^{0}\Lambda\;,\\\\[10.00002pt]
\partial_{t}\tilde{S}_{i}+\partial_{j}(\alpha\sqrt{\gamma}\;T^{j}_{\;\;i})=\frac{1}{2}\alpha\sqrt{\gamma}\;T^{\mu\nu}\;\partial_{i}g_{\mu\nu}-\alpha\sqrt{\gamma}\;u_{i}\Lambda\;,\end{gathered}$
(20)

where we choose the same emissivity $\Lambda$ as in [124]

$\Lambda=\frac{\rho_{0}}{\tau_{c}}e_{th}\;.$ (21)

Here $\tau_{c}$ is the cooling timescale.

Cooling changes the internal thermal energy over time. The governing equation
in a comoving frame can be written as [124]

$\frac{de_{th}}{d\tau}=\left[\frac{(\Gamma_{th}-1)}{\rho_{0}}\frac{d\rho_{0}}{d\tau}-\frac{1}{\tau_{c}}\right]e_{th}\;,$
(22)

where $\tau$ is the proper time. When no adiabatic contraction or expansion
takes place, then $\rho_{0}$ is constant, and the last equation describes an
exponential evolution of $e_{th}$ as a function of proper time, that is

$e_{th}\propto\exp{\left(-\frac{\tau}{\tau_{c}}\right)}\;.$ (23)

##  Appendix C Cooling tests

To validate our implementation of cooling we introduce a uniform perturbation
in pressure everywhere in a test simulation involving a stable NS. We adopt a
cold polytropic function EOS, $P=k\rho_{0}^{\Gamma}$, where $\Gamma=2$ and
$k=100\,M_{\odot}^{2}$. Equation (9) implies that any positive perturbation in
pressure adds some $P_{th}$ to the initial cold pressure. We conduct three
different simulations for a stable NS with the central density
$\rho_{c}=0.001\,M_{\odot}^{-2}$. In the first simulation, no pressure
perturbation is applied, and cooling remains inactive. In the second
simulation, we impose a small pressure perturbation with $\xi_{p}=0.01$, which
consequently, produces heat, but cooling is not active. Finally, we conduct a
simulation that includes the same pressure perturbation as the second
simulation, while cooling is active with $\tau_{c}=3.16\,t_{dyn}$. We continue
the simulations for $\sim 30\,t_{dyn}$.

Figure 14:  The value of $|e_{th}|$ is shown in 2-dimensional snapshots of
three simulations involving a stable NS governed by a $\Gamma=2$ polytropic
EOS with constant $k=100\,M_{\odot}^{2}$, and $\rho_{c}=0.001\,M_{\odot}^{-2}$
. The evolution of each simulation is illustrated in a separate row. The first
row describes an NS without any initial perturbation in pressure and without
cooling. The second row corresponds to the case with a small pressure
perturbation in the initial data and no cooling. The bottom row has the same
initial data as the second row but cooling is active. The plot demonstrates
that cooling is effective at removing excess heat in the bulk of the stellar
matter. The stellar surface discontinuity generates heat at every timestep and
would require more aggressive cooling to completely cool. However, the amount
of mass in that hotter region is negligible compared to that of the cold bulk
of the star.

Fig. 14 shows the result of these simulations for the absolute value of
$e_{th}$ at four different times. The three simulations are presented in the
three rows, respectively. Even in the first simulation, where there is no
additional thermal pressure in the initial data, small numerical errors result
in the generation of some heat that grows over time. These numerical errors
can also deplete heat, leading to negative $e_{th}$ at some points. This is
why we show the absolute value of this quantity. The initial data for this
simulation contains no initial heat, as there is no heat source of $P_{th}$ at
the beginning. In the second simulation, the perturbation $\xi_{p}=0.01$
results in heat in the initial data. This excess heat can accumulate alongside
the contribution from the numerical errors over time, resulting in a
relatively hotter configuration. The third simulation starts with the same
amount of initial heat as the second simulation, but cooling removes the extra
heat from the star. As can be seen in the bottom row of plots in Fig. 14, the
initial heat and a significant portion of the heat arising from numerical
errors are effectively removed after a few dynamical times.

We can also examine the evolution of $P_{th}$ at the center of the star in
these three simulation, as depicted in Fig. 15. In the first simulation,
$P_{th}$ increases over time as a result of numerical error. The value of
$P_{th}$ remains constant over time in the second simulation, as expected. In
the last simulation, the cooling mechanism gradually decreases the value of
$P_{th}$. After $\sim 18\,t_{dyn}$, the central value of $P_{th}$ in this
simulation drops below its corresponding value in both the other simulations,
revealing the successful removal of heat by cooling. The positive $\delta P$
also leads to a decrease in density at every point. This behavior can be seen
in Fig. 16. The positive perturbation induces more pronounced oscillations in
quantities such as density, similar to the oscillations that will be discussed
in App. D for stable stars. These oscillations remain strong even after many
$t_{dyn}$, while cooling has eliminated the majority of the heat. In this,
cooling drives the maximum density to converge to its initial value and
oscillate around it.

Figure 15:  The time evolution of $P_{th}/P$ at the center of the star is
shown for three different simulations. The blue curve corresponds to a
simulation without any perturbation or cooling. The orange curve has no
cooling, but incorporates a pressure perturbation with $\xi_{p}=0.01$. The
dark red curve corresponds to the simulation with $\xi_{p}=0.01$ but cooling
is active with $\tau_{c}=3.16\,t_{dyn}$. Figure 16:  Time evolution of the
central density for the three simulations shown in Fig. 15, with colors having
the same meaning.

The results of these simulations indicate that cooling is successfully
removing excess heat, restoring the EOS to its cold state. Having verified the
functionality of the code, we test its accuracy. In other words, cooling
should operate as modeled in App. B. To design an actual (quantitative)
cooling test, we may focus on the evolution of the quantity ${e_{th}}/{e}$ at
the center over the very early times, for a simulation with $\xi_{p}=0.01$ and
$\tau_{c}=3.16\,t_{dyn}$. If $\rho_{0}$ is fixed, the value of $e_{th}$ must
be decreasing as an exponential function with the proper time, as expressed by
Eq. (23). Unlike the previous plot, the evolving quantity here is shown as a
function of proper time $\tau$ rather than coordinate time $t$. This
adjustment is made due to the exponential behavior being dependent on proper
time. To convert the coordinate time to proper time, the value of the lapse
function is obtained at each timestep, and the conversion is achieved by

$d\tau=\alpha\,dt\;.$ (24)

Then, the result can be fitted to an exponential function characterized by a
damping timescale $\tau_{damp}$ as

$e_{th}\propto\exp{(-\frac{\tau}{\tau_{damp}})}\;.$ (25)

The solid curve in Fig. 17 shows the simulation results, and the dots
represent the fitted model based on Eq. (25) with $\tau_{damp}=3.20\,t_{dyn}$
to this curve at early times. Note that $\tau_{c}=\tau_{damp}$ only if
$d\rho_{0}/d\tau=0$. Although the density changes only little in our
simulations, it is not constant. Therefore we expect that the measured
$\tau_{damp}$ will deviate slightly from the expected value of $\tau_{c}$. In
our simulation, the difference between $\tau_{c}$ and $\tau_{damp}$ is only
$\sim 1\%$.

Figure 17:  The evolution of the quantity $e_{th}/e$ at the center of an NS
during the early times. The green curve represents the simulation results with
$\xi_{p}=0.01$ and $\tau_{c}=3.16\,t_{dyn}$. The fitted model described by Eq.
(25) is shown by dots, with the damping timescale set to
$\tau_{damp}=3.20\,t_{dyn}$.

##  Appendix D Stable stars

In this appendix, we simulate the two stable configurations C and E in Fig. 1,
representing a NS and a HS, respectively. We do this to demonstrate that our
code can evolve stable stars, and since we expect these two configurations to
be the endpoints for a collapsing unstable WD, i.e. configuration B.
Configuration C is the stable NS with a central density of
$\rho_{c}=6.04\times 10^{14}\,\textrm{g}\,\textrm{cm}^{-3}$ and a radius
$R=13.78\,\textrm{km}$. Configuration E, on the other hand, is a stable HS
with a central density of $\rho_{c}=1.57\times
10^{15}\,\textrm{g}\,\textrm{cm}^{-3}$ and a radius $R=11.40\,\textrm{km}$.

Figure 18:  Evolution of the central rest-mass density for the stable NS
configuration C (orange) and for the stable HS configuration E (blue). The
coordinate time $t$ is scaled with the dynamical timescale $t_{dyn}$, which
depends on the initial maximum density $\rho_{c}(0)$ and is different for the
two stars. The plot demonstrates that we can accurately evolve stable compact
stars with the EOSs we constructed.

For these simulations, the grid structure consists of 8 and 7 refinement
levels, respectively. In both configurations, the finest resolution is set to
100 points across the initial radius of the stars. The half-side length of the
coarsest box is $134.4R$ for the stable NS simulation and $67.2R$ for the
stable HS simulation.

Fig. 18 shows the relative deviation of the central density over time, with
respect to its initial value. This figure indicates that over $8\,t_{dyn}$,
the stable NS exhibits density fluctuations of $\sim 0.01\%$ at maximum, while
the largest amplitude of fluctuations for the HS is $\sim 0.1\%$. Due to the
transition to the more complex part of the EOS, the density evolution of the
stable HS follows a less predictable pattern.

The primary conclusion drawn from these simulations is that the stable
configurations of the TSs are dynamically stable indeed, and our code can
accurately evolve them.

##  Appendix E Constraints

To ensure the accuracy of the evolution of the Einstein equations we monitor
the Hamiltonian and momentum constraint equations (see [144] for a detailed
discussion on the 3+1 decomposition of Einstein’s equations). This is
important because after perturbing our initial data, we do not resolve the
constraints, so any constraint violations should remain small. In this
appendix, we demonstrate that this is the case.

The absolute value of the constraint violation averaged over volume as a
function of time is shown in Fig. 19 for both a stable NS (solid curves) and a
collapsing WD with cooling activated from the beginning (dotted curves). The
constraint violations are denoted by $\Delta_{k}$. Red corresponds to the
violation of the Hamiltonian constraint $\mathcal{H}$, while the momentum
constraints violations, $\mathcal{M}_{j}$ with $j=\\{x,y,z\\}$, correspond to
orange, green, and blue. $\Delta_{k}$ should remain close to zero at all
timesteps. As the star undergoes collapse, the computational errors increase,
leading to larger values of $\Delta_{k}$. Since $\Delta_{k}$ is a dimensionful
quantity, we also show $\Delta_{k}$ for the stable NS simulation discussed in
App. D for comparison. The plot demonstrates that the constraints remain well
satisfied throughout the collapse. The maximum value of $\Delta_{k}$ remains
small even in the simulations with more aggressive initial perturbations
presented in this work, and it similarly occurs during the late stages of the
collapse. The worst case in our simulations involved the HS with mass loss,
where constraint violations are still $\Delta_{k}\lesssim
10^{-9}\,M_{\odot}^{-2}$.

Figure 19:  The constraints as a function of time scaled by the initial
dynamical timescale for the Hamiltonian constraint (represented in red), and
x, y, and z components of the momentum constraint (shown in orange, green, and
blue, respectively). Solid curves stand for the simulation of the stable NS
configuration C, while dotted curves denote the collapsing WD configuration B.
The sudden rise in $\Delta_{k}$ when the initial WD enters rapid collapse is
unavoidable. However, the value of $\Delta_{k}$ remains sufficiently small
throughout the simulation.

## References

  * Lattimer [2012] J. M. Lattimer, The nuclear equation of state and neutron star masses, Ann. Rev. Nucl. Part. Sci. 62, 485 (2012), arXiv:1305.3510 [nucl-th] .
  * Drischler  _et al._ [2020] C. Drischler, R. J. Furnstahl, J. A. Melendez, and D. R. Phillips, How Well Do We Know the Neutron-Matter Equation of State at the Densities Inside Neutron Stars? A Bayesian Approach with Correlated Uncertainties, Phys. Rev. Lett. 125, 202702 (2020), arXiv:2004.07232 [nucl-th] .
  * Somasundaram  _et al._ [2023a] R. Somasundaram, I. Tews, and J. Margueron, Perturbative QCD and the neutron star equation of state, Phys. Rev. C 107, L052801 (2023a), arXiv:2204.14039 [nucl-th] .
  * Epelbaum  _et al._ [2009] E. Epelbaum, H.-W. Hammer, and U.-G. Meissner, Modern Theory of Nuclear Forces, Rev. Mod. Phys. 81, 1773 (2009), arXiv:0811.1338 [nucl-th] .
  * Machleidt and Entem [2011] R. Machleidt and D. R. Entem, Chiral effective field theory and nuclear forces, Phys. Rept. 503, 1 (2011), arXiv:1105.2919 [nucl-th] .
  * Tews  _et al._ [2013] I. Tews, T. Krüger, K. Hebeler, and A. Schwenk, Neutron matter at next-to-next-to-next-to-leading order in chiral effective field theory, Phys. Rev. Lett. 110, 032504 (2013), arXiv:1206.0025 [nucl-th] .
  * Hebeler  _et al._ [2013] K. Hebeler, J. M. Lattimer, C. J. Pethick, and A. Schwenk, Equation of state and neutron star properties constrained by nuclear physics and observation, Astrophys. J. 773, 11 (2013), arXiv:1303.4662 [astro-ph.SR] .
  * Drischler  _et al._ [2014] C. Drischler, V. Soma, and A. Schwenk, Microscopic calculations and energy expansions for neutron-rich matter, Phys. Rev. C 89, 025806 (2014), arXiv:1310.5627 [nucl-th] .
  * Hebeler  _et al._ [2015] K. Hebeler, J. D. Holt, J. Menendez, and A. Schwenk, Nuclear forces and their impact on neutron-rich nuclei and neutron-rich matter, Ann. Rev. Nucl. Part. Sci. 65, 457 (2015), arXiv:1508.06893 [nucl-th] .
  * Holt and Kaiser [2017] J. W. Holt and N. Kaiser, Equation of state of nuclear and neutron matter at third-order in perturbation theory from chiral effective field theory, Phys. Rev. C 95, 034326 (2017), arXiv:1612.04309 [nucl-th] .
  * Lonardoni  _et al._ [2020] D. Lonardoni, I. Tews, S. Gandolfi, and J. Carlson, Nuclear and neutron-star matter from local chiral interactions, Phys. Rev. Res. 2, 022033 (2020), arXiv:1912.09411 [nucl-th] .
  * Drischler  _et al._ [2021] C. Drischler, J. W. Holt, and C. Wellenhofer, Chiral Effective Field Theory and the High-Density Nuclear Equation of State, Ann. Rev. Nucl. Part. Sci. 71, 403 (2021), arXiv:2101.01709 [nucl-th] .
  * Alford  _et al._ [2022] M. G. Alford, L. Brodie, A. Haber, and I. Tews, Relativistic mean-field theories for neutron-star physics based on chiral effective field theory, Phys. Rev. C 106, 055804 (2022), arXiv:2205.10283 [nucl-th] .
  * Raithel  _et al._ [2021] C. Raithel, V. Paschalidis, and F. Özel, Realistic finite-temperature effects in neutron star merger simulations, Phys. Rev. D 104, 063016 (2021), arXiv:2104.07226 [astro-ph.HE] .
  * Raithel and Paschalidis [2023a] C. A. Raithel and V. Paschalidis, Influence of stellar compactness on finite-temperature effects in neutron star merger simulations, Phys. Rev. D 108, 083029 (2023a), arXiv:2306.13144 [astro-ph.HE] .
  * Raithel and Paschalidis [2023b] C. A. Raithel and V. Paschalidis, Detectability of Finite-Temperature Effects From Neutron Star Mergers with Next-Generation Gravitational Wave Detectors (2023b), arXiv:2312.14046 [astro-ph.HE] .
  * Li  _et al._ [1997] G.-Q. Li, C. H. Lee, and G. E. Brown, Kaon production in heavy ion collisions and maximum mass of neutron stars, Phys. Rev. Lett. 79, 5214 (1997), arXiv:nucl-th/9711002 .
  * Tsang  _et al._ [2012] M. B. Tsang _et al._ , Constraints on the symmetry energy and neutron skins from experiments and theory, Phys. Rev. C 86, 015803 (2012), arXiv:1204.0466 [nucl-ex] .
  * Horowitz [2019] C. J. Horowitz, Neutron rich matter in the laboratory and in the heavens after GW170817, Annals Phys. 411, 167992 (2019), arXiv:1911.00411 [astro-ph.HE] .
  * Huth  _et al._ [2022] S. Huth _et al._ , Constraining Neutron-Star Matter with Microscopic and Macroscopic Collisions, Nature 606, 276 (2022), arXiv:2107.06229 [nucl-th] .
  * Most  _et al._ [2023] E. R. Most, A. Motornenko, J. Steinheimer, V. Dexheimer, M. Hanauske, L. Rezzolla, and H. Stoecker, Probing neutron-star matter in the lab: Similarities and differences between binary mergers and heavy-ion collisions, Phys. Rev. D 107, 043034 (2023), arXiv:2201.13150 [nucl-th] .
  * Li  _et al._ [2023a] A. Li, G.-C. Yong, and Y.-X. Zhang, Testing the phase transition parameters inside neutron stars with the production of protons and lambdas in relativistic heavy-ion collisions, Phys. Rev. D 107, 043005 (2023a), arXiv:2211.04978 [nucl-th] .
  * Sorensen  _et al._ [2024] A. Sorensen _et al._ , Dense nuclear matter equation of state from heavy-ion collisions, Prog. Part. Nucl. Phys. 134, 104080 (2024), arXiv:2301.13253 [nucl-th] .
  * Machleidt and Sammarruca [2016] R. Machleidt and F. Sammarruca, Chiral EFT based nuclear forces: Achievements and challenges, Phys. Scripta 91, 083007 (2016), arXiv:1608.05978 [nucl-th] .
  * Tews  _et al._ [2018] I. Tews, J. Carlson, S. Gandolfi, and S. Reddy, Constraining the speed of sound inside neutron stars with chiral effective field theory interactions and observations, Astrophys. J. 860, 149 (2018), arXiv:1801.01923 [nucl-th] .
  * Nagata [2022] K. Nagata, Finite-density lattice QCD and sign problem: Current status and open problems, Prog. Part. Nucl. Phys. 127, 103991 (2022), arXiv:2108.12423 [hep-lat] .
  * Lattimer [2021] J. M. Lattimer, Neutron Stars and the Nuclear Matter Equation of State, Ann. Rev. Nucl. Part. Sci. 71, 433 (2021).
  * Shuryak [1980] E. V. Shuryak, Quantum Chromodynamics and the Theory of Superdense Matter, Phys. Rept. 61, 71 (1980).
  * McLerran [1986] L. D. McLerran, The Physics of the Quark - Gluon Plasma, Rev. Mod. Phys. 58, 1021 (1986).
  * Collins and Perry [1975] J. C. Collins and M. J. Perry, Superdense Matter: Neutrons Or Asymptotically Free Quarks?, Phys. Rev. Lett. 34, 1353 (1975).
  * Heiselberg and Hjorth-Jensen [2000] H. Heiselberg and M. Hjorth-Jensen, Phases of dense matter in neutron stars, Phys. Rept. 328, 237 (2000), arXiv:nucl-th/9902033 .
  * Blaschke  _et al._ [2013] D. Blaschke, D. E. Alvarez-Castillo, and S. Benic, Mass-radius constraints for compact stars and a critical endpoint, PoS CPOD2013, 063 (2013), arXiv:1310.3803 [nucl-th] .
  * Benic  _et al._ [2015] S. Benic, D. Blaschke, D. E. Alvarez-Castillo, T. Fischer, and S. Typel, A new quark-hadron hybrid equation of state for astrophysics - I. High-mass twin compact stars, Astron. Astrophys. 577, A40 (2015), arXiv:1411.2856 [astro-ph.HE] .
  * Paschalidis  _et al._ [2018] V. Paschalidis, K. Yagi, D. Alvarez-Castillo, D. B. Blaschke, and A. Sedrakian, Implications from GW170817 and I-Love-Q relations for relativistic hybrid stars, Phys. Rev. D 97, 084038 (2018), arXiv:1712.00451 [astro-ph.HE] .
  * Ranea-Sandoval  _et al._ [2019] I. F. Ranea-Sandoval, M. G. Orsaria, G. Malfatti, D. Curin, M. Mariani, G. A. Contrera, and O. M. Guilera, Effects of Hadron-Quark Phase Transitions in Hybrid Stars within the NJL Model, Symmetry 11, 425 (2019), arXiv:1903.11974 [nucl-th] .
  * Bozzola  _et al._ [2019] G. Bozzola, P. L. Espino, C. D. Lewin, and V. Paschalidis, Maximum mass and universal relations of rotating relativistic hybrid hadron-quark stars, Eur. Phys. J. A 55, 149 (2019), arXiv:1905.00028 [astro-ph.HE] .
  * Most  _et al._ [2019] E. R. Most, L. J. Papenfort, V. Dexheimer, M. Hanauske, S. Schramm, H. Stöcker, and L. Rezzolla, Signatures of quark-hadron phase transitions in general-relativistic neutron-star mergers, Phys. Rev. Lett. 122, 061101 (2019), arXiv:1807.03684 [astro-ph.HE] .
  * Bauswein  _et al._ [2019] A. Bauswein, N.-U. F. Bastian, D. B. Blaschke, K. Chatziioannou, J. A. Clark, T. Fischer, and M. Oertel, Identifying a first-order phase transition in neutron star mergers through gravitational waves, Phys. Rev. Lett. 122, 061102 (2019), arXiv:1809.01116 [astro-ph.HE] .
  * Chatziioannou and Han [2020] K. Chatziioannou and S. Han, Studying strong phase transitions in neutron stars with gravitational waves, Phys. Rev. D 101, 044019 (2020), arXiv:1911.07091 [gr-qc] .
  * Otto  _et al._ [2020] K. Otto, M. Oertel, and B.-J. Schaefer, Hybrid and quark star matter based on a nonperturbative equation of state, Phys. Rev. D 101, 103021 (2020), arXiv:1910.11929 [hep-ph] .
  * Gieg  _et al._ [2019] H. Gieg, T. Dietrich, and M. Ujevic, Simulating Binary Neutron Stars with Hybrid Equation of States: Gravitational Waves, Electromagnetic Signatures, and Challenges for Numerical Relativity, Particles 2, 365 (2019), arXiv:1908.03135 [gr-qc] .
  * Xie and Li [2021] W.-J. Xie and B.-A. Li, Bayesian inference of the dense-matter equation of state encapsulating a first-order hadron-quark phase transition from observables of canonical neutron stars, Phys. Rev. C 103, 035802 (2021), arXiv:2009.13653 [nucl-th] .
  * Bauswein  _et al._ [2020] A. Bauswein, S. Blacker, V. Vijayan, N. Stergioulas, K. Chatziioannou, J. A. Clark, N.-U. F. Bastian, D. B. Blaschke, M. Cierniak, and T. Fischer, Equation of state constraints from the threshold binary mass for prompt collapse of neutron star mergers, Phys. Rev. Lett. 125, 141103 (2020), arXiv:2004.00846 [astro-ph.HE] .
  * Prakash  _et al._ [2021] A. Prakash, D. Radice, D. Logoteta, A. Perego, V. Nedora, I. Bombaci, R. Kashyap, S. Bernuzzi, and A. Endrizzi, Signatures of deconfined quark phases in binary neutron star mergers, Phys. Rev. D 104, 083029 (2021), arXiv:2106.07885 [astro-ph.HE] .
  * Tan  _et al._ [2022] H. Tan, T. Dore, V. Dexheimer, J. Noronha-Hostler, and N. Yunes, Extreme matter meets extreme gravity: Ultraheavy neutron stars with phase transitions, Phys. Rev. D 105, 023018 (2022), arXiv:2106.03890 [astro-ph.HE] .
  * Takatsy  _et al._ [2023] J. Takatsy, P. Kovacs, G. Wolf, and J. Schaffner-Bielich, What neutron stars tell about the hadron-quark phase transition: A Bayesian study, Phys. Rev. D 108, 043002 (2023), arXiv:2303.00013 [astro-ph.HE] .
  * Weih  _et al._ [2020] L. R. Weih, M. Hanauske, and L. Rezzolla, Postmerger Gravitational-Wave Signatures of Phase Transitions in Binary Mergers, Phys. Rev. Lett. 124, 171103 (2020), arXiv:1912.09340 [gr-qc] .
  * Hanauske  _et al._ [2021] M. Hanauske, L. R. Weih, H. Stöcker, and L. Rezzolla, Metastable hypermassive hybrid stars as neutron-star merger remnants: A case study, Eur. Phys. J. ST 230, 543 (2021).
  * Haque  _et al._ [2024] S. Haque, R. Mallick, and S. K. Thakur, Effects of Onset of Phase Transition on Binary Neutron Star Mergers, Mon. Not. Roy. Astron. Soc. 527, 11575 (2024), arXiv:2207.14485 [astro-ph.HE] .
  * Blacker  _et al._ [2023] S. Blacker, A. Bauswein, and S. Typel, Exploring thermal effects of the hadron-quark matter transition in neutron star mergers, Phys. Rev. D 108, 063032 (2023), arXiv:2304.01971 [astro-ph.HE] .
  * Gerlach [1968] U. H. Gerlach, Equation of state at supranuclear densities and the existence of a third family of superdense stars, Phys. Rev. 172, 1325 (1968).
  * Kampfer [1981] B. Kampfer, On the Possibility of Stable Quark and Pion Condensed Stars, J. Phys. A 14, L471 (1981).
  * Glendenning and Kettner [2000] N. K. Glendenning and C. Kettner, Nonidentical neutron star twins, Astron. Astrophys. 353, L9 (2000), arXiv:astro-ph/9807155 .
  * Schertler  _et al._ [2000] K. Schertler, C. Greiner, J. Schaffner-Bielich, and M. H. Thoma, Quark phases in neutron stars and a ’third family’ of compact stars as a signature for phase transitions, Nucl. Phys. A 677, 463 (2000), arXiv:astro-ph/0001467 .
  * Ayvazyan  _et al._ [2013] N. S. Ayvazyan, G. Colucci, D. H. Rischke, and A. Sedrakian, Rotating hybrid compact stars, Astron. Astrophys. 559, A118 (2013), arXiv:1308.3053 [astro-ph.SR] .
  * Zacchi  _et al._ [2016] A. Zacchi, M. Hanauske, and J. Schaffner-Bielich, Stable hybrid stars within a SU(3) Quark-Meson-Model, Phys. Rev. D 93, 065011 (2016), arXiv:1510.00180 [nucl-th] .
  * Bejger  _et al._ [2017] M. Bejger, D. Blaschke, P. Haensel, J. L. Zdunik, and M. Fortin, Consequences of a strong phase transition in the dense matter equation of state for the rotational evolution of neutron stars, Astron. Astrophys. 600, A39 (2017), arXiv:1608.07049 [astro-ph.HE] .
  * Kaltenborn  _et al._ [2017] M. A. R. Kaltenborn, N.-U. F. Bastian, and D. B. Blaschke, Quark-nuclear hybrid star equation of state with excluded volume effects, Phys. Rev. D 96, 056024 (2017), arXiv:1701.04400 [astro-ph.HE] .
  * Alford and Sedrakian [2017] M. G. Alford and A. Sedrakian, Compact stars with sequential QCD phase transitions, Phys. Rev. Lett. 119, 161104 (2017), arXiv:1706.01592 [astro-ph.HE] .
  * Blaschke and Chamel [2018] D. Blaschke and N. Chamel, Phases of dense matter in compact stars, Astrophys. Space Sci. Libr. 457, 337 (2018), arXiv:1803.01836 [nucl-th] .
  * Alvarez-Castillo  _et al._ [2019] D. E. Alvarez-Castillo, D. B. Blaschke, A. G. Grunfeld, and V. P. Pagura, Third family of compact stars within a nonlocal chiral quark model equation of state, Phys. Rev. D 99, 063010 (2019), arXiv:1805.04105 [hep-ph] .
  * Blaschke  _et al._ [2020] D. Blaschke, A. Ayriyan, D. E. Alvarez-Castillo, and H. Grigorian, Was GW170817 a Canonical Neutron Star Merger? Bayesian Analysis with a Third Family of Compact Stars, Universe 6, 81 (2020), arXiv:2005.02759 [astro-ph.HE] .
  * Pfaff  _et al._ [2022] A. Pfaff, H. Hansen, and F. Gulminelli, Bayesian analysis of the properties of hybrid stars with the Nambu–Jona-Lasinio model, Phys. Rev. C 105, 035802 (2022), arXiv:2112.09595 [astro-ph.HE] .
  * Zacchi  _et al._ [2017] A. Zacchi, L. Tolos, and J. Schaffner-Bielich, Twin Stars within the SU(3) Chiral Quark-Meson Model, Phys. Rev. D 95, 103008 (2017), arXiv:1612.06167 [astro-ph.HE] .
  * Alvarez-Castillo  _et al._ [2016] D. Alvarez-Castillo, A. Ayriyan, S. Benic, D. Blaschke, H. Grigorian, and S. Typel, New class of hybrid EoS and Bayesian M-R data analysis, Eur. Phys. J. A 52, 69 (2016), arXiv:1603.03457 [nucl-th] .
  * Christian  _et al._ [2018] J.-E. Christian, A. Zacchi, and J. Schaffner-Bielich, Classifications of Twin Star Solutions for a Constant Speed of Sound Parameterized Equation of State, Eur. Phys. J. A 54, 28 (2018), arXiv:1707.07524 [astro-ph.HE] .
  * Montana  _et al._ [2019] G. Montana, L. Tolos, M. Hanauske, and L. Rezzolla, Constraining twin stars with GW170817, Phys. Rev. D 99, 103009 (2019), arXiv:1811.10929 [astro-ph.HE] .
  * Sieniawska  _et al._ [2019] M. Sieniawska, W. Turczanski, M. Bejger, and J. L. Zdunik, Tidal deformability and other global parameters of compact stars with strong phase transitions, Astron. Astrophys. 622, A174 (2019), arXiv:1807.11581 [astro-ph.HE] .
  * Sharifi  _et al._ [2021] Z. Sharifi, M. Bigdeli, and D. Alvarez-Castillo, Studying VLOCV twin compact stars with binary mergers, Phys. Rev. D 103, 103011 (2021), arXiv:2106.11219 [astro-ph.HE] .
  * Sen  _et al._ [2022] D. Sen, N. Alam, and G. Chaudhuri, Detailed analysis of the special points on M-R solutions of hybrid/twin stars, Phys. Rev. D 106, 083008 (2022), arXiv:2209.13890 [nucl-th] .
  * Tsaloukidis  _et al._ [2023] L. Tsaloukidis, P. S. Koliogiannis, A. Kanakis-Pegios, and C. C. Moustakidis, Twin stars as probes of the nuclear equation of state: Effects of rotation through the PSR J0952-0607 pulsar and constraints via the tidal deformability from the GW170817 event, Phys. Rev. D 107, 023012 (2023), arXiv:2210.15644 [astro-ph.HE] .
  * Annala  _et al._ [2020] E. Annala, T. Gorda, A. Kurkela, J. Nättilä, and A. Vuorinen, Evidence for quark-matter cores in massive neutron stars, Nature Phys. 16, 907 (2020), arXiv:1903.09121 [astro-ph.HE] .
  * Annala  _et al._ [2022] E. Annala, T. Gorda, E. Katerini, A. Kurkela, J. Nättilä, V. Paschalidis, and A. Vuorinen, Multimessenger Constraints for Ultradense Matter, Phys. Rev. X 12, 011058 (2022), arXiv:2105.05132 [astro-ph.HE] .
  * Li  _et al._ [2023b] J. J. Li, A. Sedrakian, and M. Alford, Ultracompact hybrid stars consistent with multimessenger astrophysics, Phys. Rev. D 107, 023018 (2023b), arXiv:2207.09798 [astro-ph.HE] .
  * Essick  _et al._ [2023] R. Essick, I. Legred, K. Chatziioannou, S. Han, and P. Landry, Phase transition phenomenology with nonparametric representations of the neutron star equation of state, Phys. Rev. D 108, 043013 (2023), arXiv:2305.07411 [astro-ph.HE] .
  * Yamamoto  _et al._ [2023] Y. Yamamoto, N. Yasutake, and T. A. Rijken, Quark phases in neutron stars consistent with implications from NICER observations, Phys. Rev. C 108, 035811 (2023), arXiv:2309.10233 [nucl-th] .
  * Lin and Steiner [2023] Z. Lin and A. Steiner, Indication of Sharp and Strong Phase-Transitions from NICER Observations (2023), arXiv:2310.01619 [astro-ph.HE] .
  * Somasundaram  _et al._ [2023b] R. Somasundaram, I. Tews, and J. Margueron, Investigating signatures of phase transitions in neutron-star cores, Phys. Rev. C 107, 025801 (2023b), arXiv:2112.08157 [nucl-th] .
  * Brandes  _et al._ [2023] L. Brandes, W. Weise, and N. Kaiser, Evidence against a strong first-order phase transition in neutron star cores: Impact of new data, Phys. Rev. D 108, 094014 (2023), arXiv:2306.06218 [nucl-th] .
  * Drago and Pagliara [2018] A. Drago and G. Pagliara, Merger of two neutron stars: predictions from the two-families scenario, Astrophys. J. Lett. 852, L32 (2018), arXiv:1710.02003 [astro-ph.HE] .
  * Burgio  _et al._ [2018] G. F. Burgio, A. Drago, G. Pagliara, H. J. Schulze, and J. B. Wei, Are Small Radii of Compact Stars Ruled out by GW170817/AT2017gfo?, Astrophys. J. 860, 139 (2018), arXiv:1803.09696 [astro-ph.HE] .
  * Essick  _et al._ [2020a] R. Essick, P. Landry, and D. E. Holz, Nonparametric Inference of Neutron Star Composition, Equation of State, and Maximum Mass with GW170817, Phys. Rev. D 101, 063007 (2020a), arXiv:1910.09740 [astro-ph.HE] .
  * De Pietri  _et al._ [2019] R. De Pietri, A. Drago, A. Feo, G. Pagliara, M. Pasquali, S. Traversi, and G. Wiktorowicz, Merger of compact stars in the two-families scenario, Astrophys. J. 881, 122 (2019), arXiv:1904.01545 [astro-ph.HE] .
  * Gorda  _et al._ [2023a] T. Gorda, K. Hebeler, A. Kurkela, A. Schwenk, and A. Vuorinen, Constraints on Strong Phase Transitions in Neutron Stars, Astrophys. J. 955, 100 (2023a), arXiv:2212.10576 [astro-ph.HE] .
  * Jie Li  _et al._ [2024] J. Jie Li, A. Sedrakian, and M. Alford, Hybrid star models in the light of new multi-messenger data, arXiv e-prints , arXiv:2401.02198 (2024), arXiv:2401.02198 [astro-ph.HE] .
  * Maslov  _et al._ [2019] K. Maslov, N. Yasutake, A. Ayriyan, D. Blaschke, H. Grigorian, T. Maruyama, T. Tatsumi, and D. N. Voskresensky, Hybrid equation of state with pasta phases and third family of compact stars, Phys. Rev. C 100, 025802 (2019), arXiv:1812.11889 [nucl-th] .
  * Glendenning [1997] N. K. Glendenning,  _Compact Stars: Nuclear Physics, Particle Physics, and General Relativity_, Astronomy and astrophysics library (Springer, 1997).
  * Weber [1999] F. Weber,  _Pulsars as astrophysical laboratories for nuclear and particle physics_ (Institute of Physics, Bristol, 1999).
  * Shapiro and Teukolsky [1983] S. L. Shapiro and S. A. Teukolsky,  _Black holes, white dwarfs, and neutron stars: The physics of compact objects_ (1983).
  * Kippenhahn  _et al._ [2012] R. Kippenhahn, A. Weigert, and A. Weiss,  _Stellar structure and evolution_, Astronomy and Astrophysics Library (Springer, 2012).
  * Espino and Paschalidis [2022] P. L. Espino and V. Paschalidis, Fate of twin stars on the unstable branch: Implications for the formation of twin stars, Phys. Rev. D 105, 043014 (2022), arXiv:2105.05269 [astro-ph.HE] .
  * Colgate  _et al._ [1980] S. A. Colgate, A. G. Petschek, and J. T. Kriese, The luminosity of type I supernovae, ApJ 237, L81 (1980).
  * Burgio  _et al._ [2021] G. F. Burgio, H. J. Schulze, I. Vidana, and J. B. Wei, Neutron stars and the nuclear equation of state, Prog. Part. Nucl. Phys. 120, 103879 (2021), arXiv:2105.03747 [nucl-th] .
  * Kojo [2021] T. Kojo, QCD equations of state and speed of sound in neutron stars, AAPPS Bull. 31, 11 (2021), arXiv:2011.10940 [nucl-th] .
  * Baym  _et al._ [2018] G. Baym, T. Hatsuda, T. Kojo, P. D. Powell, Y. Song, and T. Takatsuka, From hadrons to quarks in neutron stars: a review, Rept. Prog. Phys. 81, 056902 (2018), arXiv:1707.04966 [astro-ph.HE] .
  * Read  _et al._ [2009] J. S. Read, B. D. Lackey, B. J. Owen, and J. L. Friedman, Constraints on a phenomenologically parameterized neutron-star equation of state, Phys. Rev. D 79, 124032 (2009), arXiv:0812.2163 [astro-ph] .
  * Lindblom [2010] L. Lindblom, Spectral Representations of Neutron-Star Equations of State, Phys. Rev. D 82, 103011 (2010), arXiv:1009.0738 [astro-ph.HE] .
  * O’Boyle  _et al._ [2020] M. F. O’Boyle, C. Markakis, N. Stergioulas, and J. S. Read, Parametrized equation of state for neutron star matter with continuous sound speed, Phys. Rev. D 102, 083027 (2020), arXiv:2008.03342 [astro-ph.HE] .
  * Paschalidis  _et al._ [2011a] V. Paschalidis, Z. Etienne, Y. T. Liu, and S. L. Shapiro, Head-on collisions of binary white dwarf–neutron stars: Simulations in full general relativity, Phys. Rev. D 83, 064002 (2011a), arXiv:1009.4932 [astro-ph.HE] .
  * Colucci and Sedrakian [2013] G. Colucci and A. Sedrakian, Equation of state of hypernuclear matter: impact of hyperon–scalar-meson couplings, Phys. Rev. C 87, 055806 (2013), arXiv:1302.6925 [nucl-th] .
  * Witten [1984] E. Witten, Cosmic separation of phases, Phys. Rev. D 30, 272 (1984).
  * Haensel  _et al._ [1986] P. Haensel, J. L. Zdunik, and R. Schaefer, Strange quark stars, A&A 160, 121 (1986).
  * Alcock  _et al._ [1986] C. Alcock, E. Farhi, and A. Olinto, Strange Stars, ApJ 310, 261 (1986).
  * Zdunik [2000] J. L. Zdunik, Strange stars: Linear approximation of the EOS and maximum QPO frequency, Astron. Astrophys. 359, 311 (2000), arXiv:astro-ph/0004375 .
  * Negele and Vautherin [1973] J. W. Negele and D. Vautherin, Neutron star matter at sub-nuclear densities, Nucl. Phys. A 207, 298 (1973).
  * Baym  _et al._ [1971] G. Baym, C. Pethick, and P. Sutherland, The Ground State of Matter at High Densities: Equation of State and Stellar Models, ApJ 170, 299 (1971).
  * Alvarez-Castillo and Blaschke [2017] D. E. Alvarez-Castillo and D. B. Blaschke, High-mass twin stars with a multipolytrope equation of state, Phys. Rev. C 96, 045809 (2017), arXiv:1703.02681 [nucl-th] .
  * Demorest  _et al._ [2010] P. Demorest, T. Pennucci, S. Ransom, M. Roberts, and J. Hessels, Shapiro Delay Measurement of A Two Solar Mass Neutron Star, Nature 467, 1081 (2010), arXiv:1010.5788 [astro-ph.HE] .
  * Antoniadis  _et al._ [2013] J. Antoniadis _et al._ , A Massive Pulsar in a Compact Relativistic Binary, Science 340, 6131 (2013), arXiv:1304.6875 [astro-ph.HE] .
  * Fonseca  _et al._ [2016] E. Fonseca _et al._ , The NANOGrav Nine-year Data Set: Mass and Geometric Measurements of Binary Millisecond Pulsars, Astrophys. J. 832, 167 (2016), arXiv:1603.00545 [astro-ph.HE] .
  * Cromartie  _et al._ [2019] H. T. Cromartie _et al._ (NANOGrav), Relativistic Shapiro delay measurements of an extremely massive millisecond pulsar, Nature Astron. 4, 72 (2019), arXiv:1904.06759 [astro-ph.HE] .
  * Fonseca  _et al._ [2021] E. Fonseca _et al._ , Refined Mass and Geometric Measurements of the High-mass PSR J0740+6620, Astrophys. J. Lett. 915, L12 (2021), arXiv:2104.00880 [astro-ph.HE] .
  * Carroll [2019] S. M. Carroll,  _Spacetime and Geometry: An Introduction to General Relativity_ (Cambridge University Press, 2019).
  * Akmal  _et al._ [1998] A. Akmal, V. R. Pandharipande, and D. G. Ravenhall, The Equation of state of nucleon matter and neutron star structure, Phys. Rev. C 58, 1804 (1998), arXiv:nucl-th/9804027 .
  * Bludman and Ruderman [1968] S. A. Bludman and M. A. Ruderman, Possibility of the speed of sound exceeding the speed of light in ultradense matter, Phys. Rev. 170, 1176 (1968).
  * Alford  _et al._ [2013] M. G. Alford, S. Han, and M. Prakash, Generic conditions for stable hybrid stars, Phys. Rev. D 88, 083013 (2013), arXiv:1302.4732 [astro-ph.SR] .
  * Moustakidis  _et al._ [2017] C. C. Moustakidis, T. Gaitanos, C. Margaritis, and G. A. Lalazissis, Bounds on the speed of sound in dense matter, and neutron star structure, Phys. Rev. C 95, 045801 (2017), [Erratum: Phys.Rev.C 95, 059904 (2017)], arXiv:1608.00344 [nucl-th] .
  * Reed and Horowitz [2020] B. Reed and C. J. Horowitz, Large sound speed in dense matter and the deformability of neutron stars, Phys. Rev. C 101, 045803 (2020), arXiv:1910.05463 [astro-ph.HE] .
  * Kurkela  _et al._ [2010] A. Kurkela, P. Romatschke, and A. Vuorinen, Cold Quark Matter, Phys. Rev. D 81, 105021 (2010), arXiv:0912.1856 [hep-ph] .
  * Özel and Freire [2016] F. Özel and P. Freire, Masses, Radii, and the Equation of State of Neutron Stars, Ann. Rev. Astron. Astrophys. 54, 401 (2016), arXiv:1603.02698 [astro-ph.HE] .
  * Haensel  _et al._ [2016] P. Haensel, M. Bejger, M. Fortin, and L. Zdunik, Rotating neutron stars with exotic cores: masses, radii, stability, Eur. Phys. J. A 52, 59 (2016), arXiv:1601.05368 [astro-ph.HE] .
  * Raithel [2019] C. A. Raithel, Constraints on the Neutron Star Equation of State from GW170817, Eur. Phys. J. A 55, 80 (2019), arXiv:1904.10002 [astro-ph.HE] .
  * Radice  _et al._ [2020] D. Radice, S. Bernuzzi, and A. Perego, The Dynamics of Binary Neutron Star Mergers and GW170817, Ann. Rev. Nucl. Part. Sci. 70, 95 (2020), arXiv:2002.03863 [astro-ph.HE] .
  * Paschalidis  _et al._ [2011b] V. Paschalidis, Y. T. Liu, Z. Etienne, and S. L. Shapiro, The merger of binary white dwarf–neutron stars: Simulations in full general relativity, Phys. Rev. D 84, 104032 (2011b), arXiv:1109.5177 [astro-ph.HE] .
  * Loffler  _et al._ [2012] F. Loffler _et al._ , The Einstein Toolkit: A Community Computational Infrastructure for Relativistic Astrophysics, Class. Quant. Grav. 29, 115001 (2012), arXiv:1111.3344 [gr-qc] .
  * Werneck  _et al._ [2023] L. Werneck _et al._ , The einstein toolkit (2023), to find out more, visit http://einsteintoolkit.org.
  * Goodale  _et al._ [2003] T. Goodale, G. Allen, G. Lanfermann, J. Massó, T. Radke, E. Seidel, and J. Shalf, The Cactus framework and toolkit: Design and applications, in  _Vector and Parallel Processing – VECPAR’2002, 5th International Conference, Lecture Notes in Computer Science_ (Springer, Berlin, 2003).
  * [128] The simulation factory, https://simfactory.org.
  * Schnetter  _et al._ [2004] E. Schnetter, S. H. Hawley, and I. Hawke, Evolutions in 3-D numerical relativity using fixed mesh refinement, Class. Quant. Grav. 21, 1465 (2004), arXiv:gr-qc/0310042 .
  * Bozzola [2021] G. Bozzola, kuibit: Analyzing Einstein Toolkit simulations with Python, J. Open Source Softw. 6, 3099 (2021), arXiv:2104.06376 [gr-qc] .
  * Stergioulas and Friedman [1995] N. Stergioulas and J. L. Friedman, Comparing Models of Rapidly Rotating Relativistic Stars Constructed by Two Numerical Methods, ApJ 444, 306 (1995), arXiv:astro-ph/9411032 [astro-ph] .
  * Stergioulas [2003] N. Stergioulas, Rotating Stars in Relativity, Living Rev. Rel. 6, 3 (2003), arXiv:gr-qc/0302034 .
  * Baumgarte and Shapiro [1998] T. W. Baumgarte and S. L. Shapiro, On the numerical integration of Einstein’s field equations, Phys. Rev. D 59, 024007 (1998), arXiv:gr-qc/9810065 .
  * Shibata and Nakamura [1995] M. Shibata and T. Nakamura, Evolution of three-dimensional gravitational waves: Harmonic slicing case, Phys. Rev. D 52, 5428 (1995).
  * Sperhake [2007] U. Sperhake, Binary black-hole evolutions of excision and puncture data, Phys. Rev. D 76, 104015 (2007), arXiv:gr-qc/0606079 .
  * Bona  _et al._ [1995] C. Bona, J. Masso, E. Seidel, and J. Stela, A New formalism for numerical relativity, Phys. Rev. Lett. 75, 600 (1995), arXiv:gr-qc/9412071 .
  * Alcubierre  _et al._ [2003] M. Alcubierre, B. Bruegmann, P. Diener, M. Koppitz, D. Pollney, E. Seidel, and R. Takahashi, Gauge conditions for long term numerical black hole evolutions without excision, Phys. Rev. D 67, 084023 (2003), arXiv:gr-qc/0206072 .
  * Etienne  _et al._ [2015] Z. B. Etienne, V. Paschalidis, R. Haas, P. Mösta, and S. L. Shapiro, IllinoisGRMHD: An Open-Source, User-Friendly GRMHD Code for Dynamical Spacetimes, Class. Quant. Grav. 32, 175009 (2015), arXiv:1501.07276 [astro-ph.HE] .
  * [139] Z. B. Etienne, Cactus code thorn illinoisgrmhd, https://github.com/zachetienne/nrpytutorial/tree/master/IllinoisGRMHD.
  * Paschalidis  _et al._ [2009] V. Paschalidis, M. MacLeod, T. W. Baumgarte, and S. L. Shapiro, Merger of white dwarf-neutron star binaries: Prelude to hydrodynamic simulations in general relativity, Phys. Rev. D 80, 024006 (2009), arXiv:0910.5719 [astro-ph.HE] .
  * Arnowitt  _et al._ [2008] R. L. Arnowitt, S. Deser, and C. W. Misner, The Dynamics of general relativity, Gen. Rel. Grav. 40, 1997 (2008), arXiv:gr-qc/0405109 .
  * Duez  _et al._ [2005] M. D. Duez, Y. T. Liu, S. L. Shapiro, and B. C. Stephens, Relativistic magnetohydrodynamics in dynamical spacetimes: Numerical methods and tests, Phys. Rev. D 72, 024028 (2005), arXiv:astro-ph/0503420 .
  * Stephens  _et al._ [2008] B. C. Stephens, S. L. Shapiro, and Y. T. Liu, Collapse of magnetized hypermassive neutron stars in general relativity: Disk evolution and outflows, Phys. Rev. D 77, 044001 (2008), arXiv:0802.0200 [astro-ph] .
  * Baumgarte and Shapiro [2010] T. W. Baumgarte and S. L. Shapiro,  _Numerical Relativity: Solving Einstein’s Equations on the Computer_ (Cambridge University Press, 2010).
  * Fukushima and Sasaki [2013] K. Fukushima and C. Sasaki, The phase diagram of nuclear and quark matter at high baryon density, Prog. Part. Nucl. Phys. 72, 99 (2013), arXiv:1301.6377 [hep-ph] .
  * Voskresensky  _et al._ [2003] D. N. Voskresensky, M. Yasuhira, and T. Tatsumi, Charge screening at first order phase transitions and hadron quark mixed phase, Nucl. Phys. A 723, 291 (2003), arXiv:nucl-th/0208067 .
  * Maruyama  _et al._ [2007] T. Maruyama, S. Chiba, H.-J. Schulze, and T. Tatsumi, Hadron-quark mixed phase in hyperon stars, Phys. Rev. D 76, 123015 (2007), arXiv:0708.3277 [nucl-th] .
  * Xia  _et al._ [2019] C.-J. Xia, T. Maruyama, N. Yasutake, and T. Tatsumi, Constraining quark-hadron interface tension in the multimessenger era, Phys. Rev. D 99, 103017 (2019), arXiv:1902.08766 [hep-ph] .
  * Xia  _et al._ [2020] C.-J. Xia, T. Maruyama, N. Yasutake, T. Tatsumi, H. Shen, and H. Togashi, Systematic study on the quark-hadron mixed phase in compact stars, Phys. Rev. D 102, 023031 (2020), arXiv:2005.02273 [hep-ph] .
  * Mariani and Lugones [2023] M. Mariani and G. Lugones, Quark-hadron pasta phase in neutron stars: the role of medium-dependent surface and curvature tensions (2023), arXiv:2308.13973 [nucl-th] .
  * Essick  _et al._ [2020b] R. Essick, I. Tews, P. Landry, S. Reddy, and D. E. Holz, Direct Astrophysical Tests of Chiral Effective Field Theory at Supranuclear Densities, Phys. Rev. C 102, 055803 (2020b), arXiv:2004.07744 [astro-ph.HE] .
  * Chodos  _et al._ [1974] A. Chodos, R. L. Jaffe, K. Johnson, C. B. Thorn, and V. F. Weisskopf, A New Extended Model of Hadrons, Phys. Rev. D 9, 3471 (1974).
  * Nambu and Jona-Lasinio [1961a] Y. Nambu and G. Jona-Lasinio, Dynamical Model of Elementary Particles Based on an Analogy with Superconductivity. 1., Phys. Rev. 122, 345 (1961a).
  * Nambu and Jona-Lasinio [1961b] Y. Nambu and G. Jona-Lasinio, Dynamical model of elementary particles based on an analogy with superconductivity. II., Phys. Rev. 124, 246 (1961b).
  * Haensel  _et al._ [2007] P. Haensel, A. Y. Potekhin, and D. G. Yakovlev,  _Neutron stars 1: Equation of state and structure_, Vol. 326 (Springer, New York, USA, 2007).
  * Suwa  _et al._ [2018] Y. Suwa, T. Yoshida, M. Shibata, H. Umeda, and K. Takahashi, On the minimum mass of neutron stars, Mon. Not. Roy. Astron. Soc. 481, 3305 (2018), arXiv:1808.02328 [astro-ph.HE] .
  * Noble  _et al._ [2006] S. C. Noble, C. F. Gammie, J. C. McKinney, and L. D. Zanna, Primitive variable solvers for conservative general relativistic magnetohydrodynamics, The Astrophysical Journal 641, 626 (2006).
  * Tews  _et al._ [2019] I. Tews, J. Margueron, and S. Reddy, Confronting gravitational-wave observations with modern nuclear physics constraints, Eur. Phys. J. A 55, 97 (2019), arXiv:1901.09874 [nucl-th] .
  * Abgaryan  _et al._ [2018] V. Abgaryan, D. Alvarez-Castillo, A. Ayriyan, D. Blaschke, and H. Grigorian, Two Novel Approaches to the Hadron-Quark Mixed Phase in Compact Stars, Universe 4, 94 (2018), arXiv:1807.08034 [astro-ph.HE] .
  * Zhang and Li [2023] N.-B. Zhang and B.-A. Li, Properties of first-order hadron-quark phase transition from inverting neutron star observables, Phys. Rev. C 108, 025803 (2023), arXiv:2304.07381 [nucl-th] .
  * Harada  _et al._ [2023] R. Harada, K. Cannon, K. Hotokezaka, and K. Kyutoku, On the Testability of the Quark-Hadron Transition Using Gravitational Waves From Merging Binary Neutron Stars (2023), arXiv:2310.13603 [astro-ph.HE] .
  * Raaijmakers  _et al._ [2021] G. Raaijmakers, S. K. Greif, K. Hebeler, T. Hinderer, S. Nissanke, A. Schwenk, T. E. Riley, A. L. Watts, J. M. Lattimer, and W. C. G. Ho, Constraints on the Dense Matter Equation of State and Neutron Star Properties from NICER’s Mass–Radius Estimate of PSR J0740+6620 and Multimessenger Observations, Astrophys. J. Lett. 918, L29 (2021), arXiv:2105.06981 [astro-ph.HE] .
  * Christian and Schaffner-Bielich [2022] J.-E. Christian and J. Schaffner-Bielich, Confirming the Existence of Twin Stars in a NICER Way, Astrophys. J. 935, 122 (2022), arXiv:2109.04191 [astro-ph.HE] .
  * Jakobus  _et al._ [2021] P. Jakobus, A. Motornenko, R. O. Gomes, J. Steinheimer, and H. Stoecker, The possibility of twin star solutions in a model based on lattice QCD thermodynamics, Eur. Phys. J. C 81, 41 (2021), arXiv:2004.07026 [nucl-th] .
  * Sorkin [1981] R. Sorkin, A Criterion for the Onset of Instability at a Turning Point, ApJ 249, 254 (1981).
  * Friedman  _et al._ [1988] J. L. Friedman, J. R. Ipser, and R. D. Sorkin, Turning Point Method for Axisymmetric Stability of Rotating Relativistic Stars, ApJ 325, 722 (1988).
  * Schiffrin and Wald [2014] J. S. Schiffrin and R. M. Wald, Turning Point Instabilities for Relativistic Stars and Black Holes, Class. Quant. Grav. 31, 035024 (2014), arXiv:1310.5117 [gr-qc] .
  * Demircik  _et al._ [2022] T. Demircik, C. Ecker, and M. Järvinen, Dense and Hot QCD at Strong Coupling, Phys. Rev. X 12, 041012 (2022), arXiv:2112.12157 [hep-ph] .
  * Contrera  _et al._ [2022] G. A. Contrera, D. Blaschke, J. P. Carlomagno, A. G. Grunfeld, and S. Liebing, Quark-nuclear hybrid equation of state for neutron stars under modern observational constraints, Phys. Rev. C 105, 045808 (2022), arXiv:2201.00477 [nucl-th] .
  * Sammarruca and Millerson [2022] F. Sammarruca and R. Millerson, The Equation of State of Neutron-Rich Matter at Fourth Order of Chiral Effective Field Theory and the Radius of a Medium-Mass Neutron Star, Universe 8, 133 (2022), arXiv:2202.13475 [nucl-th] .
  * Minamikawa  _et al._ [2021] T. Minamikawa, T. Kojo, and M. Harada, Quark-hadron crossover equations of state for neutron stars: constraining the chiral invariant mass in a parity doublet model, Phys. Rev. C 103, 045205 (2021), arXiv:2011.13684 [nucl-th] .
  * Tan  _et al._ [2020] H. Tan, J. Noronha-Hostler, and N. Yunes, Neutron Star Equation of State in light of GW190814, Phys. Rev. Lett. 125, 261104 (2020), arXiv:2006.16296 [astro-ph.HE] .
  * Komoltsev and Kurkela [2022] O. Komoltsev and A. Kurkela, How Perturbative QCD Constrains the Equation of State at Neutron-Star Densities, Phys. Rev. Lett. 128, 202701 (2022), arXiv:2111.05350 [nucl-th] .
  * Gorda  _et al._ [2023b] T. Gorda, O. Komoltsev, and A. Kurkela, Ab-initio QCD Calculations Impact the Inference of the Neutron-star-matter Equation of State, Astrophys. J. 950, 107 (2023b), arXiv:2204.11877 [nucl-th] .
  * Pirner [1992] H.-J. Pirner, The Color dielectric model of QCD, Prog. Part. Nucl. Phys. 29, 33 (1992).
  * Drago  _et al._ [1995] A. Drago, M. Fiolhais, and U. Tambini, Quark matter in the chiral color dielectric model, Nucl. Phys. A 588, 801 (1995), arXiv:hep-ph/9503462 .
  * Klevansky [1992] S. P. Klevansky, The Nambu-Jona-Lasinio model of quantum chromodynamics, Rev. Mod. Phys. 64, 649 (1992).
  * Johnson [1975] K. Johnson, The M.I.T. Bag Model, Acta Phys. Polon. B 6, 865 (1975).
  * Oppenheimer and Volkoff [1939] J. R. Oppenheimer and G. M. Volkoff, On massive neutron cores, Phys. Rev. 55, 374 (1939).
  * Tolman [1939] R. C. Tolman, Static solutions of einstein’s field equations for spheres of fluid, Phys. Rev. 55, 364 (1939).
  * Ozel  _et al._ [2012] F. Ozel, D. Psaltis, R. Narayan, and A. S. Villarreal, On the Mass Distribution and Birth Masses of Neutron Stars, Astrophys. J. 757, 55 (2012), arXiv:1201.1006 [astro-ph.HE] .
  * Cierniak and Blaschke [2021] M. Cierniak and D. Blaschke, Hybrid neutron stars in the mass-radius diagram, Astron. Nachr. 342, 819 (2021), arXiv:2106.06986 [nucl-th] .
  * Sotani  _et al._ [2022] H. Sotani, N. Nishimura, and T. Naito, New constraints on the neutron-star mass and radius relation from terrestrial nuclear experiments, PTEP 2022, 041D01 (2022), arXiv:2203.05410 [nucl-th] .
  * Lattimer [2019] J. M. Lattimer, Neutron Star Mass and Radius Measurements, Universe 5, 159 (2019).
  * Steiner  _et al._ [2013] A. W. Steiner, J. M. Lattimer, and E. F. Brown, The Neutron Star Mass-Radius Relation and the Equation of State of Dense Matter, Astrophys. J. Lett. 765, L5 (2013), arXiv:1205.6871 [nucl-th] .
  * Glendenning [1992] N. K. Glendenning, First order phase transitions with more than one conserved charge: Consequences for neutron stars, Phys. Rev. D 46, 1274 (1992).
  * Li  _et al._ [2024] J. J. Li, A. Sedrakian, and M. Alford, Hybrid star models in the light of new multi-messenger data (2024), arXiv:2401.02198 [astro-ph.HE] .