11institutetext:  Institute for Theoretical Physics, University of Wrocław,
Plac Maksa Borna 9, 50-204 Wrocław, Poland 22institutetext: Departament
d’Astronomia i Astrofísica, Universitat de Valéncia, C/ Dr Moliner 50, 46100,
Burjassot (Valéncia), Spain 33institutetext: Observatori Astronòmic,
Universitat de València, E-46980, Paterna (València), Spain

# Neutron stars in accreting systems – signatures of the QCD phase transition

Noshad Khosravi Largani email<EMAIL_ADDRESS>Tobias
Fischer email<EMAIL_ADDRESS>Shota Shibagaki 11 Pablo Cerdá-
Durán 2233 Alejandro Torres-Forné 2233

(Received date / Accepted date )

Neutron stars (NS) that are born in binary systems, with a main sequence star
companion, can experience mass transfer, resulting in the accumulation of
material at the NS’s surface. This, in turn, leads to the continuous growth of
the NS mass and the associated steepening of the gravitational potential. If
the central density surpasses the onset for the phase transition from nuclear,
generally hadronic matter to deconfined quark-gluon plasma—a quantity
currently constrained solely from an upper limit by asymptotic freedom in
QCD—the system may experience a dynamic response due to the appearance of
additional degrees of freedom in the equation of state (EOS). This might give
rise to a rapid softening of the EOS during the transition in the hadron-quark
matter co-existence region. Whilst this phenomenon has long been studied in
the context of hydrostatic configurations, the dynamical implications of this
problem are yet incompletely understood. It is the essence of the present
paper to simulate the dynamics of NS, with previously accreted envelopes,
caused by the presence of a first-order QCD phase transition. Therefore, the
neutrino radiation hydrodynamics treatment is employed based on the fully
general relativistic approach in spherical symmetry, implementing three-flavor
Boltzmann neutrino transport and a microscopic model EOS that contains a
first-order hadron-quark phase transition. The associated neutrino signal
shows a sudden rise of the neutrino fluxes and average energies, becoming
observable for the present generation of neutrino detectors for a galactic
event, and a gravitational wave mode analysis reveals the behaviors of the
dominant $f$ and first gravity $g$ modes that are being excited during the NS
evolution across the QCD phase transition.

###### Key Words.:

Stars: neutron – Equation of state – Dense matter – Neutrinos – Gravitational
waves

##  1 Introduction

Accreting neutron stars (NS) in binary systems have long been considered
subjects of gravitational wave (GW) emission (c.f. Wagoner, 1984), e.g.,
through ellipticity and internal oscillations, which has been studied
extensively in the context of low-mass X-ray binary systems by Watts et al.
(2008) with potential future GW detection prospects (see also Andersson et
al., 2011). Unstable oscillation modes have been studied too in the context of
accreting NS (c.f. Lasky, 2015, and references therein) as well as GW emission
due to star quakes (see Ruderman, 1969; Giliberti & Cambiotti, 2022, and
references therein) and in addition due to high magnetic fields (c.f.
Bonazzola & Gourgoulhon, 1996, and references therein); for recent works about
magnetically deformed NS, see Haskell et al. (2008) and Singh et al. (2020).
One of the largest uncertainties in all these considerations is the yet-
incompletely understood equation of state (EOS), concerning the NS crust,
especially the stability of the crust (Chamel & Haensel, 2008), as well as the
high-density phase in excess of nuclear saturation density. In particular, the
latter concerns the presence of hadronic resonances as well as a transition to
the quark-gluon plasma. From numerical solutions of quantum chromodynamics
(QCD), the theory of strong interactions with quarks and gluons as fundamental
degrees of freedom, a cross-over transition occurs at a pseudocritical
temperature, in the range of 150–160 MeV, and at vanishing baryon density
(c.f. Bazavov et al., 2019). Even though perturbative QCD has been probed
useful to potentially constrain the high-density EOS for compact stars (c.f.
Kurkela et al., 2014; Annala et al., 2020, and references therein), the
conditions encountered at the NS interior correspond to the non-perturbative
regime of QCD and hence phenomenological quark matter models have commonly
been employed in astrophysical studies. These assume a first-order transition
with a phase transition construction from a given hadronic model EOS
(Glendenning, 2012). Observations of NS glitches, c.f. Espinoza et al. (2011),
is one example that has been speculated to be related to possible phase
transitions at the NS interior, e.g., due to the appearance of deconfined QCD
degrees of freedom as well as a transition to a state of superfluidity have
been investigated (for classical works, see Anderson & Itoh, 1975; Alpar et
al., 1984, and references therein); for a comprehensive review about the
possible mechanisms for glitches, see Haskell & Melatos (2015).

Complementary, methods of NS asteroseismology have long been explored in the
context of compact stars, implementing linear perturbation theory (see
Andersson & Kokkotas, 1998, and references therein) in order to predict
possible GW emission characteristics. Such an approach has been demonstrated
to provide reliable frequency estimates for the dominant modes, known as $f$
mode, as well as for gravity $g_{l}$ modes and pressure $p_{l}$ modes. While
the former has been referred to the damping of oscillations (c.f. Benhar et
al., 2004), the presence of $g$ (transverse) and $p$ (longitudinal) modes are
related to buoyancy and sound waves, respectively. Furthermore,
astroseismology has been proposed as a tool to probe the stellar properties,
which implies the EOS. This has been done by Sotani et al. (2012) for NS and
for proto-NS by Sotani & Takiwaki (2016) Sotani et al. (2017), Torres-Forné et
al. (2018), and Torres-Forné et al. (2019b) as well as in Sotani & Sumiyoshi
(2019) toward black-hole formation in failed core-collapse supernovae (see
also Cerdá-Durán et al., 2013). In addition, various universal relations of NS
bulk properties have been studied to remove the EOS dependence, e.g., relating
the maximum mass and corresponding radii of rapidly and non-rotating
configurations as well as the moment of inertia and NS tidal deformability
(for a recent review, see Yagi & Yunes, 2017, and references therein),
extending also to finite temperature studies (Raduta et al., 2020; Khosravi
Largani et al., 2022). Universal relations have also been suggested to exist
for several GW modes deduced from detailed astroseismology analyses in the
content of NS in Sotani & Kumar (2021) and for proto-NS by Torres-Forné et al.
(2019a), Warren et al. (2020), Sotani et al. (2021) and Mori et al. (2023).
These studies, despite some tension about the exact dependencies, e.g.,
surface gravity and mean density, will eventually reveal details about the EOS
from future GW detections.

A detailed analysis within the context of QCD phase transitions at the NS
interior has been considered within general relativistic hydrodynamics
simulations of isolated NS in axial symmetry by Dimmelmeier et al. (2009),
where the phase transition occurs due to the loss of angular momentum, and in
hydrostatic models by Sotani & Kojo (2023) probing the high-density EOS. The
latter study found that the configurations with quark matter at their interior
have generally lower dominant $f$-mode frequencies, by several 100 Hz,
compared to the hadronic reference case, depending on the compact star’s mass.
On the other hand, hydrodynamical simulations, featuring a first-order QCD
phase transition, have been performed extensively in the context of core-
collapse supernovae (Sagert et al., 2009; Fischer et al., 2011, 2018) and
binary neutron star mergers (Bauswein et al., 2019; Most et al., 2019).
Special emphasis has been in these studies on predictions for potentially
observable multi-messenger signals in neutrino and GW emissions (see also Zha
et al., 2020; Kuroda et al., 2022) as well as to constrain bulk properties
from a potential future observation of the core-collapse supernova neutrino
signal in Khosravi Largani et al. (2023) and the binary NS merger GW signals
in Blacker et al. (2020). Kuroda et al. (2022) found loud core-collapse
supernova GW signals from the first-order QCD phase transition with
frequencies in the range of several kHz with amplitudes in excess of the
ordinary neutrino-driven supernova explosions and even of black hole
formation.

In this article, we extend the study of Sotani & Kojo (2023) by performing
neutrino radiation hydrodynamics simulations of NS in accreting systems as an
alternative astrophysical scenario to probe a possible QCD phase transition. A
well calibrated nuclear matter EOS from Hempel & Schaffner-Bielich (2010)
based on the DD2F relativistic mean field parametrization of Typel et al.
(2010), with density-dependent nucleon-meson couplings, is implemented
together with a first-order phase transition to the quark-gluon plasma based
on the relativistic density function EOS of Bastian (2021). The resulting
hybrid EOS is consistent with present nuclear physics constraints such as
nuclear saturation properties and the supersaturation density constrain
derived from the elliptic flow analysis of heavy-ion collision experiments by
Danielewicz et al. (2002), as well as astrophysical constraints derived from
observations of pulsar timing.

From population synthesis it is presently understood that about 30–50% of
massive stars are in binary systems (c.f. Eldridge, 2017, and references
therein), some of which result in high-mass X-ray binaries (c.f. Tauris & van
den Heuvel, 2006, and references therein). Even though binary stellar
evolution predicts potentially very different tracks for single-and binary-
star systems (c.f. Sana et al., 2012), it is natural to assume that some of
the binary systems will end up having a NS, i.e. after the first core-collapse
supernova explosion took place and a main sequence companion star. The latter
can well be a low- or intermediate-mass star. The lifetime of such systems
depends on the properties of the system, e.g., age, metallicity, geometry such
as the spatial separation and angular velocities as well as that of the
companion star such as its mass. If the binary system’s configuration enables
mass transfer from the companion star onto the NS, also in particular if a
common envelope phase occurs (Keegans et al., 2019), which is enabled when
Roche-Lobe overflow occurs (c.f. Tauris et al., 2013, and references therein).

The details leading to accreting NS systems are left aside in the present
investigation. Instead, it is considered that a NS has undergone mass transfer
for a certain period and further, that mass accretion has ceased. This is
similar to what is considered in accreting white dwarf systems leading to
Novae (Chomiuk et al., 2021). Furthermore, it is expected that mass transfer
from the companion star heats the neutron star’s surface and, to a lesser
extent, the NS interior. In excess of temperatures of several $10^{9}$ K
($T=6\times 10^{9}$ K corresponding to $0.517$ MeV), neutrinos are being
produced from a variety of weak reactions. Taking those into account in the
simulations, featuring three-flavor Boltzmann neutrino transport allows for
the prediction not only of the neutrino fluxes and their spectra but to study
also the potential impact on mass ejection associated with the NS evolution.

The manuscript is organized as follows. In sec. 2 the simulation setup is
introduced, followed by the discussion of the neutrino radiation hydrodynamics
simulations in sec. 3, including the neutrino emission. Sec. 4 discusses the
GW mode analysis. The paper closes with a summary in sec. 5.

##  2 Equation of state and initial conditions

For the construction of the initial conditions for the numerical study of
systems that are composed of neutron stars with previously accreted envelopes,
the hadronic EOS catalog of Hempel et al. (2012) is employed together with the
relativistic density functional (RDF) quark-matter hadron EOS catalog of
Bastian (2021), both of which are available as multi-purpose tabulations for
astrophysical application at the CompOSE EOS database111The EOS data and
routines provided by CompOSE can be downloaded from the web site
https://compose.obspm.fr. For the hadronic EOS, the DD2F relativistic mean
field (RMF) model of Typel et al. (2010) is selected, and for the quark-hadron
EOS the DD2F-RDF-1.1 parametrization is employed representative for the entire
class of RDF models. At densities below normal nuclear matter density222The
saturation density for DD2F is 0.145 fm-3 or equivalently in cgs units,
$2.44\times 10^{14}$ g cm-3 and low temperatures, the RMF EOS is extended into
the modified nuclear statistical equilibrium (NSE) model of Hempel &
Schaffner-Bielich (2010) with several thousand nuclear species (for a review
of this class of hadronic EOS in supenrova simulations, see Fischer et al.,
2017, and references therein). Figure 1 (left panel) compares these two EOS at
the condition of $\beta$-equilibrium and for two constant values of the
entropy per baryon of $s=0.1\leavevmode\nobreak\ k_{\rm B}$ and
$s=0.5\leavevmode\nobreak\ k_{\rm B}$.

The DD2F-RDF-1.1 hybrid EOS of Bastian (2021) implements a first-order hadron-
quark phase transition from the DD2F hadronic EOS. At the present conditions
of $s=0.5\leavevmode\nobreak\ k_{\rm B}$ an onset density for the phase
transition is obtained of $\rho=8.8\times 10^{14}$ g cm-3 and an onset mass of
1.582 M⊙, as is illustrated via the horizontal grey line in the right panel of
Fig. 1. Note that the EOS bulk property differences for different entropies in
the range of $s=0.1$–$0.5\leavevmode\nobreak\ k_{\rm B}$, and $T=0$
configurations are on the order of one percent (see Table I in Bastian, 2021).
Larger differences arise on the phase transition region, as is illustrated in
the left panel in Fig. 1, comparing the $T=0$ and $s=0.5\leavevmode\nobreak\
k_{\rm B}$ cases where the onset density for the phase transition decreases as
well as the change of the pressure slope is less abrupt for the finite entropy
situation. This feature of the RDF class of hadron-quark matter hybrid EOS has
been discussed in detail in the example of core-collapse supernova simulations
in Fischer et al. (2018), Fischer (2021) and recently in Khosravi Largani et
al. (2023), as well as for (non)rotating hot and cold NS in Khosravi Largani
et al. (2022), however, for somewhat larger entropies of
$s=3\leavevmode\nobreak\ k_{\rm B}$ leaving a stronger impact on the hybrid
(proto)NS maximum mass and radius.

Figure 1:  Comparison of the EOS at $\beta$-equilibrium for zero temperature
$(T=0)$ and for a constant entropy per particle of $s=0.5\leavevmode\nobreak\
k_{\rm B}$, for DD2F-RDF-1.1 (solid light and dark blue lines) and DD2F
(dashed red and magenta lines), in comparison to the corresponding DD2F
hadronic EOS (red and magenta dashed lines respectively, which lay on top of
each other for both configurations), showing pressure $P$ versus restmass
density $\rho$ (left panel) and the mass $M$ radius $R$ relation (right
panel), including as example the high-precision maximum mass constraint of
$2.01\pm 0.04$ M⊙ (light grey band) deduced from the radio observations of PSR
J0348+0432 of Antoniadis et al. (2013). The horizontal grey line in the mass-
radius diagram marks the onset mass for the quark-hadron phase transition at
$M=1.582$ M⊙ for the EOS with $s=0.5\leavevmode\nobreak\ k_{\rm B}$. Figure 2:
Initial conditions for the NS simulations, based on the DD2F-RDF-1.1 hadron-
quark hybrid EOS, for a NS with total baryon mass of $M=1.781$ M⊙
(corresponding to a gravitational mass of $M_{\rm G}=1.588$ M⊙), showing the
temperature $T$ in the upper left panel, restmass density $\rho$ in the upper
right panel, entropy per baryon $s$ in the lower left panel, and electron
abundance $Y_{e}$ in the lower right panel, as a function of the enclosed
baryon mass. We compare two configurations, with constant central entropy of
$s=0.1\leavevmode\nobreak\ k_{\rm B}$ (black solid lines) and
$s=0.5\leavevmode\nobreak\ k_{\rm B}$ (grey dashed lines).

The quark matter phase of the RDF EOS is based on the concept of the string-
flip model, originally developed by Goddard et al. (1973) and Johnson & Thorn
(1976), and later further extended by Röpke et al. (1986) and Horowitz et al.
(1985). It mimics confinement through divergent medium-dependent quark masses,
which has been extended recently for applications to hydrodynamical
simulations of heavy-ion collisions in Bastian et al. (2023). Of particular
relevance for astrophysics are repulsive interactions, implemented here
through leading and higher order vector interaction channels, following Benić
et al. (2015) and Kaltenborn et al. (2017), similar as is the case in
simplistic quark-matter model EOS of the Nambu–Jona-Lasinio type (Nambu &
Jona-Lasinio, 1961) as well as vector-interaction enhanced bag models of Klähn
& Fischer (2015). These are essential in order to obtain maximum neutron star
masses—hybrid stars if quark matter is present at their interior—consistent
with the current constraints of the maximum mass of about 2 M⊙ deduced from
high-precision radio observations (c.f. Antoniadis et al., 2013; Fonseca et
al., 2021). The DD2F-RDF-1.1 hybrid EOS features a maximum mass of 2.150 M⊙
for $s=0.5\leavevmode\nobreak\ k_{\rm B}$ and 2.155 M⊙ for $T=0$ (for a
detailed discussion about the differences of finite entropy per particle and
$T=0$ NS and hybrid stars, see Khosravi Largani et al., 2022).

Precision NS radius constraints are available since the observation analysis
of quiescent low-mass X-ray binaries by Steiner et al. (2010), with
$R=11$–$12$ km for intermediate-mass neutron stars of $M\simeq 1.4$ M⊙. These
constraints have been further strengthened due to the analysis of the GW data
from the binary neutron star merger event GW170817, provided by Abbott et al.
(2018) and reanalyzed by De et al. (2018) taking into account the binary NS
mass ratio. Further, high-precision NS radii constraints are available from
the NICER NASA mission, consistent with the present DD2F-RDF-1.1 EOS. The
radii of massive neutron stars of around 2.0 M⊙ is $R_{2.0}=11.1$ km, in
agreement with the values derived by Riley et al. (2021) and Miller et al.
(2021), and that of 1.5 M⊙ intermediate-mass NS is $R_{1.5}=12.4$ km, in
agreement with the analysis of Bilous et al. (2019) and Miller et al. (2019).

Initial conditions are constructed for the hydrodynamical simulations of NS
with an accreted envelope. Therefore, hydrostatic NS solutions, based on the
Tolmann-Oppenheimer-Volkov equations, are computed for a NS of selected baryon
mass of $M=1.781$ M⊙ based on the hadronic DD2F EOS. It corresponds to a
gravitational mass of $M_{\rm G}=1.588$ M⊙. A constant central entropy of
$s=0.1\leavevmode\nobreak\ k_{\rm B}$ is assumed and compared with a
configuration of $s=0.5\leavevmode\nobreak\ k_{\rm B}$ in Fig. 2, showing the
corresponding radial profiles of selected quantities as a function of the
enclosed baryon mass $M$. It becomes evident that the structure of the NS is
independent of all thermodynamic quantities, such as the restmass density and
$Y_{e}$, for this entropy range, except central temperatures, which vary
between about 3 and 16 MeV, respectively (lower left panel in in Fig. 2).
However, the thermal contributions to the high-density EOS are negligible,
such that both configurations have the same central densities of
$\rho=8.5\times 10^{14}$ g cm-3 (upper right panel in Fig. 1). Towards the
star’s surface, the temperature decreases towards around $T=0.5$ MeV, and the
electron abundance (lower right panel) increases from the minimum value of
$Y_{e}=0.018$ towards values of $Y_{e}=0.4642$ at the lowest density
considered in this study. For simplicity, in the hydrodynamical simulations,
matter is assumed to be in nuclear statistical equilibrium everywhere in the
NS envelope and crust, which, according to modified NSE EOS of Hempel &
Schaffner-Bielich (2010), gives 56Ni as the nuclear composition. The NS crust-
core transition is therefore implemented via the excluded volume approach,
constructing therefore a first-order phase transition, which ensures the
disappearance of all nuclei at nuclear saturation density. A comparison of
this approach with other, microscopic nuclear calculations has been discussed
in Fischer et al. (2020b).

An accreted envelope is assumed to contain $M_{\rm envelope}=0.018$ M⊙,
constructed such that the resulting NS, with $M_{\rm G}=1.588$ M⊙, belongs to
the hybrid branch in the mass-radius diagram. Note that the onset mass is 1.57
M⊙ for the DD2F-RDF-1.1 EOS (see Table I in Bastian, 2021). The continues mass
transfer from the secondary companion star is expected to gradually heat the
stellar interior as the gravitational potential steepens continuously during
this process, for which we select the aforementioned values of central
entropy. However, the entropy of the outer parts of the envelope rises sharply
at low density, reaching values up to $s=1$–$3\leavevmode\nobreak\ k_{\rm B}$,
see Fig. 1, which is dominated by radiation.

##  3 Neutron star simulations with phase transition

The initial conditions constructed in Sec. 2, corresponding to the NS with
accreted envelope, are at the boundary between the hadronic phase and the
onset conditions for the hadron-quark phase transition. In fact only the
envelope mass of 0.0181 M⊙ is above the mass onset. These solutions of the
hydrostatic equations for $s=0.1\leavevmode\nobreak\ k_{\rm B}$ are
implemented into the neutrino radiation hydrodynamics model as initial
conditions, considering the impact of the hadron-quark phase transition. In
appendix A we confirm the hydrostatic equilibrium solution obtained within the
hydrodynamics approach based on the DD2F EOS, implementing a high radial
resolution of 300 mass mesh points.

Figure 3:  Evolution of 230 selected baryonic mass shells corresponding to the
material being ejected, showing the radius $R$, velocity $u$, electron
fraction $Y_{e}$ and restmass density $\rho$. Figure 4:  Evolution of selected
quantities, from top to bottom: central quark matter volume fraction
$\chi_{\rm central}^{\rm quark}$, central restmass density $\rho_{\rm
central}$, central temperature $T_{\rm central}$, central lapse function
$\alpha_{\rm central}$, and explosion energy estimate $E_{\rm diagnostic}$ as
well as the ejected baryonic mass $M_{\rm ejecta}$. Figure 5:  Evolution of
the neutrino luminosities (left panels) and average energies (right panels)
for all flavors, sampled in the co-moving frame of reference at a radius of
200 km, for the DD2F-RDF-1.1 run.

Therefore, the spherically symmetric general relativistic neutrino radiation
hydrodynamics model AGILE-BOLTZTRAN is employed. It was developed by
Liebendörfer et al. (2004), see also Mezzacappa & Bruenn (1993c, a, b). It
solves the neutrino radiation hydrodynamics equations in co-moving coordinates
(Liebendörfer et al., 2001), featuring an adaptive baryon mass mesh (see
Liebendörfer et al., 2002; Fischer et al., 2009). The general relativistic
three-flavor Boltzmann equation of Lindquist (1966) is solved for the
transport of neutrinos. The collision integral includes the list of weak
reactions considered here (see Table 1 of Fischer et al., 2020a), which
includes the Urca processes including the (inverse) neutron decay within the
full kinematics approach including weak magnetism contributions, following Guo
et al. (2020), elastic neutrino nucleon/nucleus scattering following Bruenn
(1985) with the approximate inclusion of inelastic contributions and weak
magnetism following Horowitz (2002), inelastic neutrino electron/positron
scattering of Mezzacappa & Bruenn (1993b) and neutrino pair processes. The
latter consider electron/positron annihilation, nucleon-nucleon bremsstrahlung
of Friman & Maxwell (1979) and Hannestad & Raffelt (1998), including $\pi NN$
vertex corrections (see also Bartl et al., 2016; Guo & Martínez-Pinedo, 2019)
implemented by Fischer (2016), and the annihilation of electron neutrino-
antineutrino pairs to $\mu/\tau$ neutrino-antineutrino pairs of Buras et al.
(2006) implemented in Fischer et al. (2009).

The initial neutrino phase space distributions are set to zero at all neutrino
energies and propagation angles, for all flavors. The time resolution is then
adjusted to a very short scale, with timesteps on the order of sub-
nanoseconds, to encompass all reaction rates, including the scattering and
pair reaction kernels. This leads to the launch of an artificial startup
neutrino burst, however, does not cause any impact on the evolution of the NS
structure. It is also present in the reference hadronic run, discussed in
Appendix A. The later neutrino evolution and impact of neutrino transport will
be discussed below.

Figure 3 shows the evolution of 230 selected baryonic mass shells, between
mass coordinates $M=1.776$ M⊙ and $M=1.780676$ M⊙, which corresponds to the
material being ejected in the simulations launched based on the DD2F-RDF-1.1
hybrid EOS. After an initial period during which the accreted envelope expands
slightly, as a direct response to the thermal pressure in the radiation
dominated low-density domain, the central quark matter fraction, denoted as
$\chi^{\rm quark}$, rises slowly, as illustrated in the top panel of Fig. 4
after about 3 ms, however, never reaching and remaining in the pure quark
matter phase. The latter is defined as $\chi^{\rm quark}=1$, where $\chi^{\rm
quark}=0$ corresponds to matter being in the hadronic phase and in the hadron-
quark mixed phase $0<\chi^{\rm qurk}<1$ (further details about the definitions
of $\chi^{\rm quarj}$ can be found in Bastian, 2021). The matter velocities,
shown in Fig. 3, oscillate accordingly on the order of a few times $10^{4}$ km
s-1 reaching in the peaks $1.5\times 10^{5}$ km s-1.

After the softening of the central NS EOS, due to the onset of the phase
transition—note that the hybrid EOS softens significantly in the hadron-quark
mixed phase where the pressure slope decreases (see Fig. 1)—the central
density rises rapidly. It oscillates in accordance with the central quark
matter fraction around a mean value of about $10^{15}$ g cm-3, as is
illustrated in the top panel of Fig. 4. The onset of the phase transition at
the NS interior is accompanied by a decreasing central temperature, as shown
in Fig. 4. This phenomenon is known for this class of RDF hybrid EOS, i.e.
during the adiabatic compression of the central NS fluid elements, the
temperature decreases. It has been discussed at the example of core-collapse
supernovae in Fischer et al. (2018) and Fischer (2021). In the hadron-quark
mixed phase, after around 3 ms, the temperature of the central fluid elements
oscillate accordingly, though on a larger magnitude than before reaching the
mixed phase, similar to the central temperatures, due to the substantially
larger EOS gradients encountered in the mixed phase.

As a consequence of the evolution after reaching quark matter at the neutron
star interior, the steeper gravitational potential results in increased
temperatures at the neutron star surface where the accreted envelope expands,
reaching supersonic matter velocities. This rapid initial expansion is
illustrated in Fig. 3 (top panel). The entire envelope expands quickly to
radii on the order of several 100 km. At about 15 ms, the outer envelope
reaches a radius of about 400 km. This evolution is accompanied by the
constant change of the gravitational potential, illustrated via the central
lapse function $\alpha_{\rm central}$ in Fig. 4 dropping from 0.62 to slightly
below 0.5, and hence oscillations of the neutron star surface, which in turn,
affects the central evolution as shown in Fig. 4. The associated millisecond
timescale of these oscillations is related to the high central densities on
the order of $10^{15}$ g cm-3 and the CFL conditions for the numerical scheme
solving the system of a partial differential equation, i.e. the time
resolution decreases by a factor of 2 when the density increases from
$\rho=10^{14}$ g cm-3 to $\rho=10^{15}$ g cm-3.

At the moment of the hadron-quark phase transition, when the NS envelope
expands, positive diagnostic explosion energies are obtained. These are
evaluated by integrating the total energy, composed of gravitational, kinetic,
and internal energies, the latter given by the EOS, from the surface towards
the center, following thereby the standard procedure of core-collapse
supernova phenomenology (c.f. Fischer et al., 2009, and references theren).
Values of about $E_{\rm diagnostic}\simeq 5\times 10^{49}$ erg are obtained at
about 30 ms of simulation time, as illustrated in Fig. 4 (bottom panel),
however, still the asymptotic value has not been reached yet. The total mass
ejected is about $M_{\rm ejecta}\simeq 1\times 10^{-3}$ M⊙, for the simulation
time. These ejecta feature high escape velocities, on the order of 10–60% of
the speed of light (see Fig. 3).

It is interesting to note that these ejecta are subject to neutrino
interactions as they expand, which give rise to a broad $Y_{e}$ distribution
(see Fig. 3). There is an early component featuring slightly neutron-rich
conditions with $Y_{e}\simeq 0.3-0.45$ and high entropies per baryon on the
order of $s\simeq 100-175\leavevmode\nobreak\ k_{\rm B}$ as illustrated in the
bottom panel of Fig. 3. It is followed by an isospin symmetric component with
$Y_{e}\simeq 0.5$ and moderately low entropies per baryon of $s\simeq
50\leavevmode\nobreak\ k_{\rm B}$. The neutrino luminosities and average
energies are shown in Fig. 5 for all flavors, sampled in the co-moving
reference frame at a radius of 200 km. This corresponds to the free-streaming
regime. Luminosities rise suddenly after about 6 ms, coinciding when quark
matter appears at the NS interior, reaching values of a few times $10^{52}$
erg s-1 with average energies on the order of several tens of MeV. Initially,
i.e. before the appearance of quark matter, the neutrino luminosities are on
the order of $10^{49}$ erg s-1 and the average energies are on the order of
few MeV, as is discussed and illustrated in Appendix A. Note further, the
reason for the early ejecta to remain low $Y_{e}$ is due to the high expansion
velocities on the order of 60% of the speed of light, while for the later
ejecta the neutrino fluxes and their average energies are responsible for the
rising $Y_{e}$, similar as has been found semi-analytically for the late-time
supernova ejecta by Qian & Woosley (1996). Since the spectral and flux
differences remain generally small, in particular, the difference $\langle
E_{\bar{\nu}_{e}}\rangle-\langle E_{\nu_{e}}\rangle$ remains small, on the
order of few MeV, despite the generally high average energies on the order of
15–25 MeV for both $\nu_{e}$ and $\bar{\nu}_{e}$, the spectral difference is
insufficient to keep the expanding material neutron-rich.

(a) Eigenmodes including $f$ and $g_{1}$ mode frequency fits (b) Mass, radius
and mass-radius relations

Figure 6:  Results of the mode analysis (left panel), showing the solutions of
the GREAT analysis (black dots) and the fit expression (2) for the dominant
mode (solid purple line), see text for details, as well as the corresponding
evolution of the NS (right panel), enclosed mass $M_{\rm NS}$ and radius
$R_{\rm NS}$, sampled at a restmass density of $\rho=10^{11}$ g cm-3. Figure
7:  Radial profiles of selected quantities, from top to bottom, showing
restmass density $\rho$, radial solutions of the eigenfrequencies $\eta_{r}$
for the highest three modes found, indicating the corresponding frequencies,
further the Lamb $L$ and Brunt-Väisälä $N$ frequencies as well as $Y_{e}$ and
speed of sound $cs$, at three different times (see Fig. 6), corresponding to
the hadronic evolution in graph (a), the transition when quark matter appears
in graph (b) and at late times in graph (c).

##  4 Mode analysis

In this section results of an in-depth mode analysis are conducted following
the approach of Torres-Forné et al. (2019b) employing the numerical GREAT
code. It solves the eigenvalue problem of hydrodynamic and metric
perturbations of a spherically symmetric self-gravitating system equilibrium
model in General Relativity. Assuming spherical symmetry, these eigenmodes are
obtained through linear adiabatic perturbation analysis of the hydrodynamics
and metric equations with a background in equilibrium. This linear Eulerian
perturbation is expanded with a harmonic time dependence of frequency and in
spherical harmonics for the angular dependence (see Torres-Forné et al.,
2019b). The corresponding Lagrangian fluid element displacement vector can
then be expressed as follows,

$\displaystyle\boldsymbol{\xi}$ $\displaystyle=$
$\displaystyle\bigg{[}\eta_{r}(r)Y_{lm}(\theta,\varphi)\boldsymbol{\hat{r}}$
(1)
$\displaystyle+\,\eta_{\perp}(r)\left(\frac{\partial_{\theta}Y_{lm}(\theta,\varphi)}{r^{2}}\boldsymbol{\hat{\theta}}+\frac{\partial_{\varphi}Y_{lm}(\theta,\varphi)}{r^{2}\sin^{2}\theta}\boldsymbol{\hat{\varphi}}\right)\bigg{]}e^{i\sigma
t}$

where we use spherical coordinates $(r,\theta,\varphi)$ with its corresponding
orthonormal basis
$(\boldsymbol{\hat{r}},\boldsymbol{\hat{\theta}},\boldsymbol{\hat{\varphi}})$
and the spherical harmonics $Y_{lm}$. Here, $\sigma$ is the oscillation
frequency in the harmonic time-dependent part. $\eta_{r}$ and $\eta_{\perp}$
are two functions of the radius that characterize the radial and angular
components of the displacement, respectively.

Employing the relativistic Brunt-Väisälä frequency and the relativistic Lamb
frequency, denoted as $N^{2}$ and $L^{2}$, respectively (see Eqs. (33) and
(34) in Torres-Forné et al., 2018), is a valuable approach to identify regions
of possible convection. The former is related to the criterion of Ledoux
stability and quantified by the presence of lepton number and entropy
gradients (c.f. Mirizzi et al., 2016, and references therein). In general,
$N^{2}$ governs gravity ($g$) modes and $L^{2}$ dictates pressure ($p$) modes
and the fundamental ($f$) mode.

Distinct modes can emerge contingent upon the temporal evolution of diverse
physical parameters. Gravity modes manifest in regions where buoyancy serves
as a restoring force, characterized by $N^{2}>0$, indicating stability to
convective motion. In the absence of buoyancy, $p$ modes are excited by the
propagation of sound waves. Our focus lies not in the modes supported by sound
waves, but rather in the interior of the NS. Although the NS is not fully
isolated, and the choice of boundary conditions influences $p$-mode
eigenfrequency calculations, our attention is directed solely toward the NS
interior. We impose boundary conditions from the innermost fluid elements to
the surface of the NS. It is crucial to note that the definition of the NS
surface is based on the density profile reaching $\rho=10^{11}$ g cm-3. The
inner boundary is set to be the central fluid element. Altering the inner
boundary to somewhat lower densities, e.g., $\rho=10^{14}$ g cm-3, has a
negligible impact on the solutions, given the lack of significant density
fluctuations within the NS inner core.

The frequency outcomes derived from mode analysis via the GREAT tool are
depicted in the left panel of Fig. 6. These results encompass solutions
obtained for all nodes with $l=2$, spanning the frequency range of 0 to 4 kHz.
Due to the chosen boundary conditions, the majority of modes associated with
$p$ modes fall beyond this range. Following the mode classification of Cowling
(1941), the foremost dominant mode, characterized by zero nodes and identified
as the $f$ mode, initiates at approximately 1.9 kHz. Subsequently, reaching
the onset density for the first-order QCD phase transition, it experiences a
slight reduction and transiently follows a descending trajectory. At lower
frequencies, another prominent mode with $n=1$ node emerges as the first $g$
mode, with a broader frequency spectrum compared to the $f$ mode and
distinctive behavior. It commences at approximately 0.5 kHz and exhibits a
sudden, drastic rise to frequencies around 1 kHz after reaching the onset
density for the QCD phase transition. Despite the neutron star never fully
reaching a quark matter core and undergoing oscillations in mass and radius,
the g-mode, on the whole, evolves toward higher frequencies.

In accordance with the node definition proposed by Torres-Forné et al. (2018),
the quantity $n$ is herein defined as the count of sign changes within the
radial function $\eta_{r}$. These are shown in Fig 7 at three different times,
corresponding to the hadronic evolution before the appearance of quark matter
(left panel), during the transition when quark matter appeared at the NS
interior (middle panel) and during the long-term evolution after (right
panel). The associated frequencies are listed in the legend, coinciding with
the definition of Torres-Forné et al. (2018), confirming the highest frequency
mode as $f$ mode and the other two lower frequencies as $g$ modes. Radial
profiles of both, $N$ and $L$ frequencies, are shown in Fig. 7 in the panels
below. Note that the plane-wave approximation imposes the condition,
$L^{2}>N^{2}$, for acceptable solutions. In other words, oscillatory solutions
can occur when $\sigma^{2}>N^{2},L^{2}$ (representing acoustic waves), or when
$\sigma^{2}<N^{2},L^{2}$ (indicating gravity waves). Intermediate values lead
to non-oscillatory evanescent solutions. It is essential to acknowledge that
this simplifying assumption may not be applicable in the present context. The
NS system features finite entropy, potentially deviating from the conditions
under which this approximation is valid. In particular, it becomes evident
from Fig. 7, that the rise of $N^{2}$ is dominated by the presence of a large
lepton number gradients, indicated by the rising $Y_{e}$ profile at radii
above around 12 km, since lepton number and entropy gradients at the NS
interior are negligible as, in particular, the initial conditions assume
constant entropy per particle at the NS interior. Entropy fluctuations remain
overall small during the entire simulation and the speed of sound, $c_{s}$,
shown in the bottom panel of Fig. 7 is small, $(c_{s}/{\rm c})^{2}\simeq
0.1$–$0.3$, and nearly constant towards larger radii. At the interior,
however, the value of the speed of sound exceeds the perturbative QCD limit of
$(c_{s}/{\rm c})^{2}=1/3$, since the conditions correspond to the highly non-
perturbative regime of QCD. Note further the evolution of the $f$ and $g_{1}$
frequencies via the associated radial solutions $\eta_{r}$ in Fig. 7, which
quantifies the results shown in Fig. 6 at the three selected times.

Table 1: Parameters for the frequency fits expression (2). mode | $x$ | $a_{0}$ | $a_{1}$ | $a_{2}$
---|---|---|---|---
$f$ | $\sqrt{M_{\rm NS}/R_{\rm NS}^{3}}$ | $1470$ | $1.5\times 10^{4}$ | $0$
$g_{1}$ | $M_{\rm NS}/R_{\rm NS}^{2}$ | 1075 | $8\times 10^{3}$ | $-5.4\times 10^{6}$

In order to fit the frequency results to the evolution of NS bulk properties,
such as the enclosed mass and its radius, we closely follow the discussion in
Torres-Forné et al. (2019a), and express linear and quadratic frequency fits
for the dominant and the first sub-dominant modes, respectively, in terms of
the mean density inside the NS, $x\equiv\sqrt{M_{\rm NS}/R_{\rm NS}^{3}}$, and
in terms of surface gravity, $x\equiv M_{\rm NS}/R_{\rm NS}^{2}$, as follows:

$\omega_{\rm fit}(x)=a_{0}+a_{1}\,x+a_{2}\,x^{2}\leavevmode\nobreak\ ,$ (2)

with the fit parameters $a_{0}$, $a_{1}$ and $a_{2}$ are given in Table 1 for
the dominant and the first sub-dominant modes identified. This indicates that
the frequency fits found, the dominant one initially at around 2 kHz and
descending towards 1.6 kHz after the phase transition is an $f$ mode, and the
one starting at around 500 Hz and the rapidly rising frequency to about 1000
Hz, is the $g_{1}$ mode, the magenta and turquoise lines, respectively. The
corresponding evolution of the enclosed NS mass and radius are shown in the
right panel of Fig. 6, sampled at a restmass density of $\rho=10^{11}$ g cm-3,
to avoid $p$ modes. The general tendency, i.e. the ratio of $M/R$ decreases
after the appearance of quark matter at the NS interior, is due to the
ejection of material from the NS surface as a consequence of the dynamical
response of the system to the central NS oscillations discussed in sec. 3. All
other sub-dominant modes are higher $g_{l}$-modes.

##  5 Summary and conclusions

The present paper discusses NS in accreting systems as a novel astrophysical
scenario for the appearance of QCD degrees of freedom. To this end, general
relativistic neutrino radiation hydrodynamics simulations are performed based
on three flavor Boltzmann neutrino transport in spherical symmetry. Initial
conditions are constructed based on the well selected DD2F nuclear EOS for a
NS with a total baryonic mass of 1.781 M⊙, featuring an accreted envelope of
$M_{\rm envelope}=0.018$ M⊙, corresponding to the hybrid branch of the mass-
radius relation of the DD2F-RDF-1.1 hybrid EOS employed in this study. It
features a first-order phase transition based on the Maxwell construction from
the DD2F RMF nuclear matter EOS to the RDF-1.1. model. The latter includes
vector repulsion which is responsible for the stiffening of the EOS with
increasing density and hence yielding maximum hybrid star masses, which are
consistent with the current constraint of about 2 M⊙. The appearance of quark
matter at the NS interior causes the NS’s central density to rapidly rise and
oscillate on a timescale on the order of milliseconds. In response to the
central contraction, matter is expelled from the hybrid star’s surface,
resulting in the ejection of approximately $10^{-3}$ M⊙ of baryonic mass
during the considered simulation times. Explosion energy estimates are
obtained of $5\times 10^{49}$ erg, computed following the procedure commonly
used in the core-collapse supernova context.

Due to the sudden compression and later, long-term oscillations of the hybrid
star interior, that lasts at least for several tens of milliseconds, the
hybrid star experiences a temperature increase, in particular at low densities
towards the surface. This, in turn, results in the production of neutrinos
with luminosities on the order of $10^{52}$ erg s-1 and average energies of
several tens of MeV. This represents an observable signal for the current
generation of kiloton water-Cherenkov radiation detectors for a galactic
event.

The Urca processes, electron captures on protons and positron capture on
neutrons, as well as their reverse $\nu_{e}$ and $\bar{\nu}_{e}$ captures, and
the (inverse)neutron decay, cause the rise of the $Y_{e}$ of the ejected
material, from $Y_{e}\leq 0.46$ to $Y_{e}\simeq 0.5$. However, the release of
a millisecond burst-like neutrino signature, as was reported for the QCD phase
transition in the context of core-collapse supernovae (c.f. Fischer et al.,
2018, and references therein), has not been found here. On the other hand, a
detailed GW mode analysis shows the evolution of the decreasing dominant
$f$-mode frequency from around 1.9 kHz, in the hadronic phase, towards around
1.7 kHz, after the appearance of quark matter. The magnitude of the decrease
reported here for the particular NS mass is consistent with the findings of
Sotani & Kojo (2023), who further studied a range of compact stars between
1.0–2.3 M⊙ comparing models with QCD degrees of freedom to a reference
hadronic case, finding the increase of the decreasing $f$-mode frequencies for
more massive stars. The present study reports in addition about the $g_{1}$
mode evolution, for which we find a sudden rise from 0.5 kHz, in the hadronic
phase, to about 1 kHz after the appearance of quark matter. The lower
$g_{l}$-mode frequencies, i.e. $l\geq 2$, follow qualitatively the behaviour
of $g_{1}$. We provide linear and quadratic fits of the $f$\- and $g_{1}$-mode
frequencies, following the argumentation of Torres-Forné et al. (2019a), in
terms of the mean density and surface gravity, respectively.

It is interesting to note that a complete phase transition, i.e. the central
NS reaches and stays in the pure quark matter, was not found, even though the
initial conditions belong to the hybrid branch of the mass-radius relation, by
an amount of exactly $M_{\rm envelope}$. Mass transfer in binary systems—from
a main sequence star to a secondary NS—results in the slow accumulation of
material at the NS surface and hence the growth of the NS total mass. A low-
mass NS born in a core-collapse supernova explosion without QCD phase
transition will slowly move along the mass-radius relation towards increasing
mass. The present investigation demonstrates quantitatively, for a particular
hadron-quark hybrid EOS, the hydrodynamical evolution when such system reaches
the conditions for the onset of the QCD phase transition. It remains to be
explored what is the role of the EOS, in terms of varying onset densities for
quark matter as well as bulk properties of the phase transition such as the
pressure slope and other dependencies of the hadron-quark mixed phase, such as
the possible presence of shear (Sotani et al., 2013) and associated presence
of geometric structures collectively known as ’Pasta’ phases (c.f. Yasutake et
al., 2014, and references therein), which we leave for future exploration.

Furthermore, the impact of potentially high magnetic fields and rapid
rotation, associated with the very existence of magnetars, require the multi-
dimensional magneto neutrino radiation hydrodynamics framework. The present
work is limited to NS with low surface magnetic fields, referring to old NS.
However, the response of the system reported here due to the QCD phase
transition is expected to remain qualitatively, independent of the magnetic
field, featuring the ejection of material from the NS surface and the rise of
the neutrino fluxes to observable values as well as the emission of GW with
kHz frequencies. The magnitude of the frequencies might alter as well as the
amount of mass ejected as the oscillations reported here will be damped by the
emission of GW in multi-dimensional simulations, as was demonstrated in the
multi-dimensional hydrodynamics simulations of Dimmelmeier et al. (2009),
based on polytropic parametrised and MIT bag model EOS, settling the NS back
into the initial stable state on a longer timescale.

## Acknowledgements

The authors (N.K.L., T.F., and S.S.) acknowledge support from the Polish
National Science Centre (NCN) under grant number 2020/37/B/ST9/00691. P.C.D
and A.T.F acknowledge support from the Spanish Agencia Estatal de
Investigación (Grants No. PGC2018-095984-B-I00 and PID2021-125485NB-C21)
funded by MCIN/AEI/10.13039/501100011033 and ERDF A way of making Europe, by
the Generalitat Valenciana (PROMETEO/2019/071), and by COST Actions CA16104
and CA16214.

## ORCID iDs

Noshad Khosravi Largani https://orcid.org/0000-0003-1551-0508
Tobias Fischer https://orcid.org/0000-0003-2479-344X Shota Shibagaki
https://orcid.org/0000-0002-8196-0745
Pablo Cerdá-Durán https://orcid.org/0000-0001-7354-6897
Alejandro Torres-Forné https://orcid.org/0000-0001-8709-5118

##  Appendix A Hydrostatic equilibrium solution

Figure 8:  Evolution of mass shells for the reference DD2F hadronic model,
showing the evolution of (from top to bottom) radius, velocity, restmass
density, temperature and electron fraction.

To quantify hydrostatic equilibrium, neutrino radiation hydrodynamics
simulations are launched for the reference DD2F hadronic EOS, featuring a high
radial resolution of 300 radial mass shells implemented in AGILE-BOLTZTRAN.
Otherwise, the same input physics is used as discussed in Sec. 3. After an
initial relaxation phase lasts for about 5–10 ms, the central compact NS
settles into complete hydrostatic equilibrium. The results are illustrated in
Fig. 8 for all 300 mass shells, showing the evolution of selected quantities
(from top to bottom), radius, velocity, restmass density, temperature and
$Y_{e}$. The evolution of all central mass shells remains constant values for
all quantities, corresponding to the neutron star interior for densities as
low as $\rho\simeq 10^{6}$ g cm-3. The accreted envelope, on the other hand,
does not obey complete hydrostatic equilibrium, illustrated via expanding and
contracting mass shells, up to radii of about 40 km and the velocities
correspondingly reach values of a few times $10^{4}$ km/s. Densities and
temperatures follow accordingly, oscillating around $\rho\simeq 10^{4}$ g cm-3
and $T\simeq 0.2$ MeV. The oscillations are damped towards later times on the
order of several tens of milliseconds.

Figure 9:  Radial profiles of selected quantities, velocity $u$ in units of
the speed of light ${\rm c}$, restmass density $\rho$, temperature $T$,
electron abundance $Y_{e}$, and neutrino luminosities $L_{\nu}$ and root mean
square energies $\langle E\rangle_{\rm rms}$ of all flavors, comparing the
reference DD2F hadronic simulation at $t=0$ and $t\simeq 38$ ms of simulation
time.

All matter remains gravitationally bound at all times during the simulation
and hence no positive explosion energy estimate is obtained. The situation is
illustrated via radial profiles of selected quantities in Fig. 9, comparing
the initial conditions (thick grey lines) and at about 38 ms of simulation
time (thin red lines), showing the expansion of the envelope to radii of about
100 km and the subsequent density, temperature and $Y_{e}$ profiles.

Figure 10:  Evolution of neutrino luminosities (left panels) and average
energies (right panels) for all $\nu_{e}$ and $\bar{\nu}_{e}$ (top panels) and
heavy lepton flavors (bottom panels), for the DD2F hadronic reference
simulations (grey lines), in comparison to the diffusion approximation (red
lines) following Expression (3). Figure 11:  Evolution of the neutrinospheres
$R_{\nu}$ of all flavors.

The evolution of the neutrino luminosities and average energies are shown in
Fig. 10 for all flavors. The values of the luminosities remain generally low,
on the order of a few $10^{49}$ erg s-1, with average energies of $\langle
E_{\nu}\rangle\simeq 4$ MeV, after the initial relaxation phase of about 10
ms. The neutrinos stem from the neutron star’s surface, where the neutrino
luminosities rise. Values of the neutrinoshere radii, $R_{\nu}$, are given in
Fig. 11, following the standard definition via the optical depth $\tau_{\nu}$,
where $R_{\nu}:=r(\tau_{\nu}\equiv 2/3)$, which settle towards constant values
between $R_{\nu_{\mu/\tau}}\simeq 10.75$ km, $R_{\nu_{e}}\simeq 10.81$ and
$R_{\bar{\nu}_{e}}=10.85$ after the initial relaxation period. The radial
profile of the luminosities and root-mean-square energies $\langle
E\rangle_{\rm rms}$ are illustrated in the bottom panels of Fig. 9, which is
the transition between neutron star crust and previously accreted envelope,
corresponding to the sharp rise of the envelope temperature.

The Urca processes drive the electron abundance of the expanding envelope from
initially $Y_{e}\simeq 0.46$ towards $Y_{e}\simeq 0.5$, and even above at
later times, as shown in Fig. 9. The increasing $Y_{e}$ of the crust-envelope
interface is also shown in the evolution of the mass shells in the bottom
panel of Fig. 8. This is the reason for the hierarchy of the neutrino
luminosities, $L_{\nu_{e}}\gtrsim L_{\bar{\nu}_{e}}$, and average energies,
$\langle E_{\nu_{e}}\rangle\gtrsim\langle E_{\bar{\nu}_{e}}\rangle$, however,
comparable to other weak processes, which is the cause for the heavy lepton
flavor neutrino luminosities and average energies being of similar order as
those of the electron flavors. Note that the neutrino luminosities approach
asymptotically the diffusion limit,

$L_{\rm T}\simeq 4\pi\sigma_{\rm
SB}\,\kappa_{\nu}\,R_{\nu}^{2}\,T_{\nu}^{4}\leavevmode\nobreak\ ,$ (3)

where $\kappa_{\nu}$ is an empirical factor that is of order $20\%$; $1/6$ for
$\nu_{e}$, $10/52$ for $\bar{\nu}_{e}$ and $10/46$ for the heavy lepton
neutirno flavors. The values of the neutrino luminosities (3) are shown in
Fig. 10 as red lines, for comparison. This approximation has been explored in
the core-collapse supernova case for both the fail branch in Fischer et al.
(2009) as well as for the proto-neutron star deleptnoization and later cooling
phases in Fischer et al. (2012). Here, particularly during the early
evolution, (3) particularly reproduced the behavior of the heavy lepton
neutrino flavor. The evolution of the neutrinospheres of the last scattering
is shown in Fig. 11 for which the method based on the integration of the
transport mean free path of Fischer et al. (2012) has been employed.

## References

  * Abbott et al. (2018) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2018, Phys. Rev. Lett., 121, 161101
  * Alpar et al. (1984) Alpar, M. A., Langer, S. A., & Sauls, J. A. 1984, ApJ, 282, 533
  * Anderson & Itoh (1975) Anderson, P. W. & Itoh, N. 1975, Nature, 256, 25
  * Andersson et al. (2011) Andersson, N., Ferrari, V., Jones, D. I., et al. 2011, General Relativity and Gravitation, 43, 409
  * Andersson & Kokkotas (1998) Andersson, N. & Kokkotas, K. D. 1998, MNRAS, 299, 1059
  * Annala et al. (2020) Annala, E., Gorda, T., Kurkela, A., Nättilä, J., & Vuorinen, A. 2020, Nature Physics, 16, 907
  * Antoniadis et al. (2013) Antoniadis, J., Freire, P. C. C., Wex, N., et al. 2013, Science, 340, 448
  * Bartl et al. (2016) Bartl, A., Bollig, R., Janka, H. T., & Schwenk, A. 2016, Phys. Rev. D, 94, 083009
  * Bastian (2021) Bastian, N.-U. F. 2021, Phys. Rev. D, 103, 023001
  * Bastian et al. (2023) Bastian, N.-U. F., Huovinen, P., & Nazarova, E. 2023, Nucl. Phys. A, 1038, 122718
  * Bauswein et al. (2019) Bauswein, A., Bastian, N.-U. F., Blaschke, D. B., et al. 2019, Phys. Rev. Lett., 122, 061102
  * Bazavov et al. (2019) Bazavov, A., Ding, H. T., Hegde, P., et al. 2019, Physics Letters B, 795, 15
  * Benhar et al. (2004) Benhar, O., Ferrari, V., & Gualtieri, L. 2004, Phys. Rev. D, 70, 124015
  * Benić et al. (2015) Benić, S., Blaschke, D., Alvarez-Castillo, D. E., Fischer, T., & Typel, S. 2015, A&A, 577, A40
  * Bilous et al. (2019) Bilous, A. V., Watts, A. L., Harding, A. K., et al. 2019, ApJ, 887, L23
  * Blacker et al. (2020) Blacker, S., Bastian, N.-U. F., Bauswein, A., et al. 2020, Phys. Rev. D, 102, 123023
  * Bonazzola & Gourgoulhon (1996) Bonazzola, S. & Gourgoulhon, E. 1996, A&A, 312, 675
  * Bruenn (1985) Bruenn, S. W. 1985, ApJS, 58, 771
  * Buras et al. (2006) Buras, R., Janka, H.-T., Rampp, M., & Kifonidis, K. 2006, A&A, 457, 281
  * Cerdá-Durán et al. (2013) Cerdá-Durán, P., DeBrye, N., Aloy, M. A., Font, J. A., & Obergaulinger, M. 2013, ApJ, 779, L18
  * Chamel & Haensel (2008) Chamel, N. & Haensel, P. 2008, Living Reviews in Relativity, 11, 10
  * Chomiuk et al. (2021) Chomiuk, L., Metzger, B. D., & Shen, K. J. 2021, ARA&A, 59, 391
  * Cowling (1941) Cowling, T. G. 1941, MNRAS, 101, 367
  * Danielewicz et al. (2002) Danielewicz, P., Lacey, R., & Lynch, W. G. 2002, Science, 298, 1592
  * De et al. (2018) De, S., Finstad, D., Lattimer, J. M., et al. 2018, Phys. Rev. Lett., 121, 091102
  * Dimmelmeier et al. (2009) Dimmelmeier, H., Bejger, M., Haensel, P., & Zdunik, J. L. 2009, MNRAS, 396, 2269
  * Eldridge (2017) Eldridge, J. J. 2017, in Handbook of Supernovae, ed. A. W. Alsabti & P. Murdin, 671
  * Espinoza et al. (2011) Espinoza, C. M., Lyne, A. G., Stappers, B. W., & Kramer, M. 2011, MNRAS, 414, 1679
  * Fischer (2016) Fischer, T. 2016, A&A, 593, A103
  * Fischer (2021) Fischer, T. 2021, EuPhJA, 57, 270
  * Fischer et al. (2017) Fischer, T., Bastian, N.-U., Blaschke, D., et al. 2017, PASA, 34, e067
  * Fischer et al. (2018) Fischer, T., Bastian, N.-U. F., Wu, M.-R., et al. 2018, Nature Astronomy, 2, 980
  * Fischer et al. (2020a) Fischer, T., Guo, G., Dzhioev, A. A., et al. 2020a, Phys. Rev. C, 101, 025804
  * Fischer et al. (2012) Fischer, T., Martínez-Pinedo, G., Hempel, M., & Liebendörfer, M. 2012, Phys. Rev. D, 85, 083003
  * Fischer et al. (2011) Fischer, T., Sagert, I., Pagliara, G., et al. 2011, ApJS, 194, 39
  * Fischer et al. (2020b) Fischer, T., Typel, S., Röpke, G., Bastian, N.-U. F., & Martínez-Pinedo, G. 2020b, Phys. Rev. C, 102, 055807
  * Fischer et al. (2009) Fischer, T., Whitehouse, S. C., Mezzacappa, A., Thielemann, F.-K., & Liebendörfer, M. 2009, A&A, 499, 1
  * Fonseca et al. (2021) Fonseca, E., Cromartie, H. T., Pennucci, T. T., et al. 2021, ApJ, 915, L12
  * Friman & Maxwell (1979) Friman, B. L. & Maxwell, O. V. 1979, ApJ, 232, 541
  * Giliberti & Cambiotti (2022) Giliberti, E. & Cambiotti, G. 2022, MNRAS, 511, 3365
  * Glendenning (2012) Glendenning, N. 2012, Compact Stars: Nuclear Physics, Particle Physics and General Relativity, Astronomy and Astrophysics Library (Springer New York)
  * Goddard et al. (1973) Goddard, P., Goldstone, J., Rebbi, C., & Thorn, C. B. 1973, Nuclear Physics B, 56, 109
  * Guo & Martínez-Pinedo (2019) Guo, G. & Martínez-Pinedo, G. 2019, ApJ, 887, 58
  * Guo et al. (2020) Guo, G., Martínez-Pinedo, G., Lohs, A., & Fischer, T. 2020, Phys. Rev. D, 102, 023037
  * Hannestad & Raffelt (1998) Hannestad, S. & Raffelt, G. 1998, ApJ, 507, 339
  * Haskell & Melatos (2015) Haskell, B. & Melatos, A. 2015, International Journal of Modern Physics D, 24, 1530008
  * Haskell et al. (2008) Haskell, B., Samuelsson, L., Glampedakis, K., & Andersson, N. 2008, MNRAS, 385, 531
  * Hempel et al. (2012) Hempel, M., Fischer, T., Schaffner-Bielich, J., & Liebendörfer, M. 2012, ApJ, 748, 70
  * Hempel & Schaffner-Bielich (2010) Hempel, M. & Schaffner-Bielich, J. 2010, Nucl. Phys. A, 837, 210
  * Horowitz (2002) Horowitz, C. J. 2002, Phys. Rev. D, 65, 043001
  * Horowitz et al. (1985) Horowitz, C. J., Moniz, E. J., & Negele, J. W. 1985, Phys. Rev. D, 31, 1689
  * Johnson & Thorn (1976) Johnson, K. & Thorn, C. B. 1976, Phys. Rev. D, 13, 1934
  * Kaltenborn et al. (2017) Kaltenborn, M. A. R., Bastian, N.-U. F., & Blaschke, D. B. 2017, Phys. Rev. D, 96, 056024
  * Keegans et al. (2019) Keegans, J., Fryer, C. L., Jones, S. W., et al. 2019, MNRAS, 485, 620
  * Khosravi Largani et al. (2023) Khosravi Largani, N., Fischer, T., & Bastian, N. U. F. 2023, arXiv e-prints, arXiv:2304.12316
  * Khosravi Largani et al. (2022) Khosravi Largani, N., Fischer, T., Sedrakian, A., et al. 2022, MNRAS, 515, 3539
  * Klähn & Fischer (2015) Klähn, T. & Fischer, T. 2015, ApJ, 810, 134
  * Kurkela et al. (2014) Kurkela, A., Fraga, E. S., Schaffner-Bielich, J., & Vuorinen, A. 2014, ApJ, 789, 127
  * Kuroda et al. (2022) Kuroda, T., Fischer, T., Takiwaki, T., & Kotake, K. 2022, ApJ, 924, 38
  * Lasky (2015) Lasky, P. D. 2015, PASA, 32, e034
  * Liebendörfer et al. (2004) Liebendörfer, M., Messer, O. E. B., Mezzacappa, A., et al. 2004, ApJS, 150, 263
  * Liebendörfer et al. (2001) Liebendörfer, M., Mezzacappa, A., & Thielemann, F.-K. 2001, Phys. Rev. D, 63, 104003
  * Liebendörfer et al. (2002) Liebendörfer, M., Rosswog, S., & Thielemann, F.-K. 2002, ApJS, 141, 229
  * Lindquist (1966) Lindquist, R. W. 1966, Annals of Physics, 37, 487
  * Mezzacappa & Bruenn (1993a) Mezzacappa, A. & Bruenn, S. W. 1993a, ApJ, 405, 669
  * Mezzacappa & Bruenn (1993b) Mezzacappa, A. & Bruenn, S. W. 1993b, ApJ, 410, 740
  * Mezzacappa & Bruenn (1993c) Mezzacappa, A. & Bruenn, S. W. 1993c, ApJ, 405, 637
  * Miller et al. (2019) Miller, M. C., Lamb, F. K., Dittmann, A. J., et al. 2019, ApJ, 887, L24
  * Miller et al. (2021) Miller, M. C., Lamb, F. K., Dittmann, A. J., et al. 2021, ApJ, 918, L28
  * Mirizzi et al. (2016) Mirizzi, A., Tamborra, I., Janka, H.-T., et al. 2016, Nuovo Cimento Rivista Serie, 39, 1
  * Mori et al. (2023) Mori, M., Suwa, Y., & Takiwaki, T. 2023, Phys. Rev. D, 107, 083015
  * Most et al. (2019) Most, E. R., Papenfort, L. J., Dexheimer, V., et al. 2019, Phys. Rev. Lett., 122, 061101
  * Nambu & Jona-Lasinio (1961) Nambu, Y. & Jona-Lasinio, G. 1961, PhRv, 122, 345
  * Qian & Woosley (1996) Qian, Y. Z. & Woosley, S. E. 1996, ApJ, 471, 331
  * Raduta et al. (2020) Raduta, A. R., Oertel, M., & Sedrakian, A. 2020, MNRAS, 499, 914
  * Riley et al. (2021) Riley, T. E., Watts, A. L., Ray, P. S., et al. 2021, ApJ, 918, L27
  * Röpke et al. (1986) Röpke, G., Blaschke, D., & Schulz, H. 1986, Phys. Rev. D, 34, 3499
  * Ruderman (1969) Ruderman, M. 1969, Nature, 223, 597
  * Sagert et al. (2009) Sagert, I., Fischer, T., Hempel, M., et al. 2009, Phys. Rev. Lett., 102, 081101
  * Sana et al. (2012) Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Science, 337, 444
  * Singh et al. (2020) Singh, N., Haskell, B., Mukherjee, D., & Bulik, T. 2020, MNRAS, 493, 3866
  * Sotani & Kojo (2023) Sotani, H. & Kojo, T. 2023, Phys. Rev. D, 108, 063004
  * Sotani & Kumar (2021) Sotani, H. & Kumar, B. 2021, Phys. Rev. D, 104, 123002
  * Sotani et al. (2017) Sotani, H., Kuroda, T., Takiwaki, T., & Kotake, K. 2017, Phys. Rev. D, 96, 063005
  * Sotani et al. (2013) Sotani, H., Maruyama, T., & Tatsumi, T. 2013, Nucl. Phys. A, 906, 37
  * Sotani et al. (2012) Sotani, H., Nakazato, K., Iida, K., & Oyamatsu, K. 2012, Phys. Rev. Lett., 108, 201101
  * Sotani & Sumiyoshi (2019) Sotani, H. & Sumiyoshi, K. 2019, Phys. Rev. D, 100, 083008
  * Sotani & Takiwaki (2016) Sotani, H. & Takiwaki, T. 2016, Phys. Rev. D, 94, 044043
  * Sotani et al. (2021) Sotani, H., Takiwaki, T., & Togashi, H. 2021, Phys. Rev. D, 104, 123009
  * Steiner et al. (2010) Steiner, A. W., Lattimer, J. M., & Brown, E. F. 2010, ApJ, 722, 33
  * Tauris et al. (2013) Tauris, T. M., Langer, N., Moriya, T. J., et al. 2013, ApJ, 778, L23
  * Tauris & van den Heuvel (2006) Tauris, T. M. & van den Heuvel, E. P. J. 2006, in Compact stellar X-ray sources, Vol. 39, 623–665
  * Torres-Forné et al. (2019a) Torres-Forné, A., Cerdá-Durán, P., Obergaulinger, M., Müller, B., & Font, J. A. 2019a, Phys. Rev. Lett., 123, 051102
  * Torres-Forné et al. (2018) Torres-Forné, A., Cerdá-Durán, P., Passamonti, A., & Font, J. A. 2018, MNRAS, 474, 5272
  * Torres-Forné et al. (2019b) Torres-Forné, A., Cerdá-Durán, P., Passamonti, A., Obergaulinger, M., & Font, J. A. 2019b, MNRAS, 482, 3967
  * Typel et al. (2010) Typel, S., Röpke, G., Klähn, T., Blaschke, D., & Wolter, H. H. 2010, Phys. Rev. C, 81, 015803
  * Wagoner (1984) Wagoner, R. V. 1984, ApJ, 278, 345
  * Warren et al. (2020) Warren, M. L., Couch, S. M., O’Connor, E. P., & Morozova, V. 2020, ApJ, 898, 139
  * Watts et al. (2008) Watts, A. L., Krishnan, B., Bildsten, L., & Schutz, B. F. 2008, MNRAS, 389, 839
  * Yagi & Yunes (2017) Yagi, K. & Yunes, N. 2017, Phys. Rep, 681, 1
  * Yasutake et al. (2014) Yasutake, N., Łastowiecki, R., Benić, S., et al. 2014, Phys. Rev. C, 89, 065803
  * Zha et al. (2020) Zha, S., O’Connor, E. P., Chu, M.-c., Lin, L.-M., & Couch, S. M. 2020, Phys. Rev. Lett., 125, 051102