# Stability of the protoneutron stars toward black hole formation

Hajime Sotani1,2 and Kohsuke Sumiyoshi3
1Astrophysical Big Bang Laboratory, RIKEN, Saitama 351-0198, Japan
2Interdisciplinary Theoretical & Mathematical Science Program (iTHEMS), RIKEN,
Saitama 351-0198, Japan
3National Institute of Technology, Numazu College, Ooka 3600, Numazu, Shizuoka
410-8501, Japan<EMAIL_ADDRESS>
###### Abstract

We examine the protoneutron star (PNS) stability in this study by solving the
radial oscillation equations. For this purpose, we adopt the numerical results
of massive PNS toward the black hole formation obtained by spherically
symmetric numerical simulations for core-collapse supernova with general
relativistic neutrino-radiation hydrodynamics. We find that the PNSs are
basically stable in their evolution against the radial perturbations, while
the PNS finally becomes unstable before the apparent horizon appears inside
the PNS. We also examine the gravitational wave frequencies from PNS with the
relativistic Cowling approximation. Then, we derive the empirical formula for
the $f$-mode frequency, which weakly depends on the PNS models. This kind of
universality tells us the PNS property, which is a combination of the PNS mass
and radius in this study, once one would observe the $f$-mode gravitational
waves.

###### keywords:

stars: neutron – gravitational waves – equation of state – stars: oscillations

##  1 Introduction

Success of the direct detection of gravitational waves opens a new era, i.e.,
the gravitational wave astronomy, where gravitational waves can be used for
extracting a celestial information together with the electromagnetic waves and
neutrino signals. In fact, the gravitational waves from the binary black hole
mergers and the binary neutron star merger(s) have already detected by the
Laser Interferometer Gravitational wave Observatory (LIGO) (Aasi et al., 2015)
and Virgo (Acernese et al., 2014), which tells us the mass of the compact
objects and constraint on the equation of state (EOS) for dense matter. In
particular, when the gravitational wave has been detected from the binary
neutron star merger, GW170817, the electromagnetic counterpart has also been
observed (Abbott et al., 2017a, b). Owing to this observation, it becomes
obvious that the short gamma-ray burst comes from the binary neutron star
merger. In addition, the Japanese gravitational wave detector, KAGRA (Akutsu
et al., 2020), has finally joined the LIGO and Virgo collaboration network at
the end of the third observing run, while the discussion for the third-
generation detectors, such as Einstein Telescope and Cosmic Explorer (Punturo,
Lück & Beker, 2014; Abbott et al., 2017c), has already been started. Thus, the
number of gravitational wave events will increase more and more in the near
future.

Supernova explosions must be another promising source of gravitational waves.
The gravitational waves from core-collapse supernova have been mainly
discussed with the numerical simulations of gravitational collapse of massive
star, the bounce of central core, and subsequent accreting matter to the
central compact object, i.e., protoneutron star (PNS), produced via supernova
(e.g., Murphy, Ott & Burrows (2009); Müller, Janka & Marek (2013); Ott et al.
(2013); Cerdá-Durán et al. (2013); Yakunin et al. (2015); Kuroda, Kotake &
Takiwaki (2016); Andresen et al. (2017); O’Connor & Couch (2018)). For
instance, the three-dimensional (3D) numerical simulations for the core-
collapse supernova tell us the existence of gravitational wave signals related
to the oscillations of PNS (or core region) (Kuroda, Kotake & Takiwaki, 2016;
Radice et al., 2019; Andresen et al., 2019; Mezzacappa et al., 2020), where
the frequency increases with time from a few hertz up to kilohertz after
corebounce. This signal is originally considered to be the evidence of the
Brunt-Väisälä frequency at the PNS surface, i.e., the so-called surface
gravity ($g$-) mode (Müller, Janka & Marek, 2013; Cerdá-Durán et al., 2013),
but it may be natural to consider that the signals come from the global
oscillations of the PNSs because the Brunt-Väisälä frequency is a local value
and depends strongly on the definition of PNS radius. In fact, the
gravitational wave signals appearing in the numerical simulations can be
identified with the global oscillations, such as fundamental ($f$-mode)
oscillations of PNS (Morozova et al., 2018; Sotani et al., 2019; Sotani &
Takiwaki, 2020b, c) or the $g$-mode like oscillations on the region inside the
shock radius (Torres-Forné et al., 2018, 2019a, 2019b). In addition to the
ramp up signals, the 3D simulations tell us the existence of another
gravitational wave signal with $\sim 100$ Hz, which seems to be associated
with the standing accretion-shock instability (Kuroda, Kotake & Takiwaki,
2016; Andresen et al., 2017; O’Connor & Couch, 2018; Radice et al., 2019;
Andresen et al., 2019; Mezzacappa et al., 2020).

The numerical simulations must be a suitable technique for examining the
dynamics of astronomical phenomena, while the perturbative approach is also
powerful technique for extracting the physics behind the numerical results
obtained by simulations. In particular, since a frequency from the
astronomical objects strongly depend on their properties, one would inversely
know the properties of objects through the observation of the frequency. This
technique is known as asteroseismology, which is similar to seismology on the
Earth and helioseismology on the Sun. As for (cold) neutron stars, the crust
properties are constrained by identified the quasiperiodic oscillations
observed in the afterglow following magnetar giant flares with the crustal
torsional oscillations (e.g., Gearheart et al. (2011); Sotani et al. (2012,
2013a, 2013b); Sotani, Iida & Oyamatsu (2016, 2017, 2018, 2019)).
Additionally, once the gravitational waves from compact objects will be
detected, one could constrain the stellar mass, radius, and EOS of source
objects (e.g., Andersson & Kokkotas (1996, 1998); Sotani, Tominaga & Maeda
(2001); Sotani, Kohri & Harada (2004); Tsui & Leung (2005); Sotani et al.
(2011); Passamonti & Andersson (2012); Doneva et al. (2013); Sotani (2020a, b,
2021)).

On the other hand, the study on the PNS asteroseismology is relatively limited
(e.g., Ferrari, Miniutti & Pons (2003); Burgio et al. (2011); Fuller et al.
(2015); Sotani & Takiwaki (2016); Camelio et al. (2017); Sotani et al. (2017);
Morozova et al. (2018); Torres-Forné et al. (2018, 2019a, 2019b); Sotani et
al. (2019); Sotani & Takiwaki (2020a, b, c); Bizouard et al. (2021)), compared
to that for cold neutron stars. Nevertheless, owing to the development of the
numerical simulations of core-collapse supernova, the number of studies about
the PNS asteroseismology gradually increases. Up to now, considering the PNS
asteroseismology, two different approaches have been adopted. With one of two
approaches, the PNS surface is defined with a specific surface density, which
is typically $\sim 10^{11}$ g/cm3, and the linear analysis is done inside the
PNS (Morozova et al., 2018; Sotani et al., 2019; Sotani & Takiwaki, 2020a, b,
c). With this approach, there is uncertainty how to select the surface
density, but at least the $f$\- and $g_{1}$-mode frequencies for late phase
after corebounce ($\sim 0.3$ sec.) are almost independent of the surface
density (Morozova et al., 2018; Sotani & Takiwaki, 2020b). As an advantage of
this approach, since the boundary condition is the same as the standard
asteroseismology, the classification of eigenmode is simple, i.e., the
eigenmode is classified by counting the radial nodal number in the
eigenfunction. Meanwhile, with the other approach, the linear analysis is done
in the region inside the shock radius, i.e., there is no uncertainty for
selection of the outer boundary (Torres-Forné et al., 2018, 2019a, 2019b;
Bizouard et al., 2021). But then, the zero radial displacement condition is
imposed at the outer boundary, which is different boundary condition from the
standard asteroseismology. That is, the problem to solve with the second
approach is completely different mathematical problem, where one has to
reclassify the eigenmode with another definition. So, in this study we simply
adopt the first approach and will discuss the stability of PNS and
gravitational wave frequencies.

Most of the PNS asteroseismology among previous studies are considered with
the PNSs produced by the successful supernova, while in this study we focus on
the accreting PNSs toward black hole formation due to the massive progenitor
(Sumiyoshi et al., 2006; O’Connor & Ott, 2011; Pan et al., 2018; Kuroda et
al., 2018), i.e., the case for the failed supernova, as in Sotani & Sumiyoshi
(2019). In the previous study (Sotani & Sumiyoshi, 2019), we examine the
gravitational wave frequency from the accreting PNSs and propose the
possibility how to identify the PNS model with the increasing mass by using
the simultaneous observation of the neutrino signal and gravitational waves
from PNS (See Yokozawa et al. (2015); Kuroda et al. (2017); Takiwaki & Kotake
(2018); Westernacher-Schneider et al. (2019) for the simultaneous observations
of the neutrinos and gravitational waves, for example). However, the details
of the gravitational waves from PNSs in the final phase just before the black
hole formation are not discussed well in our previous study. In addition, the
stability analysis on PNSs has been limited, although the presence of unstable
quasi-radial modes has actually been discussed in the context of supernova
simulations in Cerdá-Durán et al. (2013); Torres-Forné et al. (2019a), and
Torres-Forné et al. (2019a) computed the frequencies of the radial modes
(among others) and estimated their gravitational wave emission depending on
the deformation of the rotating PNS. So, in this study we mainly examine the
stability of PNS by considering the radial oscillation, as in cold neutron
stars (Chandrasekhar, 1964). Then, we also discuss the gravitational waves
from the PNSs, especially focusing on the final phase just before the black
hole formation, where we study the dependence of the gravitational wave
frequencies on the surface density and derive a new fitting formula by adding
new time periods from the previous analysis in Sotani & Sumiyoshi (2019).

The paper is organized as follows. In section 2, we describe the setup of PNS
models, which become background models for linear analysis, and their
properties, especially focusing on the final stage in their evolution just
before the black hole formation. In section 3, we examine the stability of
protoneutron stars by considering the radial perturbations. In section 4, we
report the eigenfrequencies of gravitational waves radiating from the PNS
models. Finally, we summarize this study in section 5. Unless otherwise
mentioned, we adopt geometric units in the following, $c=G=1$, where $c$
denotes the speed of light, and the metric signature is $(-,+,+,+)$.

##  2 PNS Models

In order to make a linear analysis, first we have to prepare the PNS models as
a background. The PNS properties depend on not only the density and pressure
profiles but also the distributions of temperature (or entropy per baryon) and
electron fraction inside the PNS, while such profiles can be determined only
via the numerical simulation of the core-collapse supernova explosion. In this
study, as in Sotani & Sumiyoshi (2019), we particularly adopt the profiles
obtained via the numerical simulations performed by solving the general
relativistic neutrino-radiation hydrodynamics under the spherical symmetry. In
the simulations, hydrodynamics and neutrino transfer in general relativity are
solved simultaneously (Yamada, 1997; Yamada, Janka & Suzuki, 1999; Sumiyoshi
et al., 2005). To describe the neutrino transfer, the Boltzmann equation is
directly solved with the multiangle and multienegy neutrino distributions for
four species, $\nu_{e}$, $\bar{\nu}_{e}$, $\nu_{\mu/\tau}$, and
$\bar{\nu}_{\mu/\tau}$, i.e., we implement 6 species of neutrinos by assuming
$\mu$-type and $\tau$-type (anti-)neutrinos have identical distributions. For
the collision term associated with neutrino emission, absorption, and
scattering with leptons, nucleons and nuclei handles, the basic neutrino
reactions are adopted (Sumiyoshi et al., 2005; Bruenn, 1985). The metric
adopted in the numerical code is given by

$ds^{2}=-e^{2\Phi(t,m_{b})}dt^{2}+e^{2\Lambda(t,m_{b})}dm_{b}^{2}+r^{2}(t,m_{b})(d\theta^{2}+\sin^{2}\theta
d\phi^{2}),$ (1)

where $t$ and $m_{b}$ denote the coordinate time and the baryon mass
coordinate, respectively (Misner & Sharp, 1964). In addition, $m_{b}$ is
related to the circumference radius ($r$) via the baryon mass conservation,
while the metric functions, $\Phi(t,m_{b})$ and $\Lambda(t,m_{b})$, are
evolved together with hydrodynamical variables in the numerical simulations
(Yamada, 1997). The numerical simulations for core-collapse supernovae have
been done with 255 grid points in the radial mass coordinate, 6 grid points in
the neutrino angle, and 14 grid points in the neutrino energy. The rezoning of
radial mesh is made during the simulations to resolve the accreting matter. We
remark that the radial grids of mass coordinate are nonuniformly arranged to
cover not only the dense region inside the central object but also the region
for accreting matter.

As in Sotani & Sumiyoshi (2019), in order to examine the massive PNS and see
the dependence on the progenitor models, we focus on two different progenitor
models, i.e., a $40M_{\odot}$ star based on Woosley & Weaver (1995) and a
$50M_{\odot}$ star based on Tominaga, Umeda & Nomoto (2007), which are
hereafter refereed to as W40 and T50, respectively (See Nakazato et al.
(2013); Sukhbold et al. (2016); Horiuchi et al. (2018) for the various models
of progenitors and their consequences). In addition, we adopt three different
EOSs in this study, i.e., the Shen EOS (Shen et al., 1998), LS180, and LS220
(Lattimer & Swesty, 1991). The Shen EOS is based on the relativistic mean
field theory with the TM1 nuclear interaction, while LS180 and LS220 are
constructed with the compressible liquid drop model. In Table 1 we show the
EOS parameters adopted in this study, i.e., the incompressibility, $K_{0}$,
the symmetry energy, $S_{0}$, and the slope parameter of the nuclear symmetry
energy, $L$, together with the maximum (gravitational) mass for the cold
neutron stars constructed with the corresponding EOS (Oertel et al., 2017).
Even though the maximum mass with LS180 can not reach the $2M_{\odot}$
observations (Demorest et al., 2010; Antoniadis et al., 2013; Cromartie et
al., 2020), we also adopt LS180 to examine the dependence of the gravitational
wave frequencies on the EOS stiffness (Sumiyoshi et al., 2004). On the other
hand, according to the terrestrial experimental constraints on the nuclear
saturation parameters, $K_{0}$ is to be $230\pm 40$ MeV (Khan & Margueron,
2003), while $L$ still seems to have a large uncertainty, such as $L=58.9\pm
16$ MeV (Li et al., 2019) and $L=106\pm 37$ MeV (Reed et al., 2021). The three
EOSs have been popularly used to explore the EOS dependence, although some of
the EOS parameters are extreme in the experimental ranges. Anyway, one can see
the details about the evolution of PNS toward the black hole formation and
properties of neutrino emissions with these progenitor models and EOSs in
Sumiyoshi, Yamama & Suzuki (2007, 2008); Sumiyoshi et al. (2009); Nakazato et
al. (2010). With using the numerical data obtained in these simulations, we
will prepare the PNS models and consider their oscillations in this study. As
in Sotani & Sumiyoshi (2019), we consider the PNS models provided by five sets
of numerical simulations, i.e., W40 with Shen, LS180, and LS220 (three sets)
and T50 with Shen and LS180 (two sets). To identify the PNS model, we
hereafter label the model by combining the progenitor model with EOS, e.g.,
W40-Shen for the PNS model provided with W40 and Shen.

Table 1: EOS parameters adopted in this study, and the maximum (gravitational) mass for a cold neutron star constructed with the corresponding EOS. EOS |  $K_{0}$ (MeV) |  $S_{0}$ (MeV) |  $L$ (MeV) | $M_{\rm max}/M_{\odot}$
---|---|---|---|---
Shen | 281 | 36.9 | 110.8 | 2.2
LS180 | 180 | 29.3 | 73.8 | 1.8
LS220 | 220 | 29.3 | 73.8 | 2.0
Figure 1: Time evolution of the PNS radius (top) and gravitational mass
(bottom) is shown as a function of $T_{\rm pb}$ for the model of W40-LS180
with $\rho_{\rm s}=10^{10}$ (open-circles), $10^{11}$ (filled-circles), and
$2\times 10^{11}$ g/cm3 (filled-squars). In the top panel, an enlarged figure
focusing on the final phase before the gravitational collapse is also added.

Unlike the cold neutron stars, the stellar surface of the PNS is ambiguous
because the accreting matter exists even outside the PNS “surface”. So,
usually the PNS surface is identified with the position where the density
becomes a specific surface density, $\rho_{s}$, e.g., typically
$\rho_{s}\simeq 10^{11}$ g/cm3. In Fig. 1, we show the $\rho_{s}$-dependence
of the PNS radius (top panel) and gravitational mass (bottom panel) for the
model of W40-LS180 with $\rho_{s}=10^{10}$, $10^{11}$, and $2\times 10^{11}$
g/cm3, where $T_{\rm pb}$ denotes a time after the core bounce in the unit of
seconds, i.e., the moment at the core bounce corresponds to $T_{\rm pb}=0$.
From this figure, one can observe a wiggling behavior in the PNS radius with
lower surface density. This is because the numerical simulations have been
done with the Lagrangian coordinate in baryon mass instead of Eulerian radial
coordinate. Resolving the gradually accreting matter in the moving baryon mass
grids is well done by the rezoning method, however, determining the radial
position is not as easy as in the fixed radial grids. As a result, the number
of grid points in the low density region is limited, i.e., the resolution is
not enough to get the smooth radius evolution with lower surface density.
Nevertheless, since the radius evolution for the PNS model with
$\rho_{s}=2\times 10^{11}$ g/cm3 is relatively smooth even in the final phase,
we adopt $\rho_{s}=2\times 10^{11}$ g/cm3 in this study and discuss the
stability and gravitational wave frequencies on such PNS models. We remark
that the PNS mass seems to be almost independent of the selection of surface
density.

In Fig. 2, we show the time evolution of the PNS radius (top panel),
gravitational mass (middle panel), and average density defined by $M_{\rm
PNS}/R_{\rm PNS}^{3}$ (bottom panel) as a function of $T_{\rm pb}-T_{\rm BH}$
for various models with $\rho_{s}=2\times 10^{11}$ g/cm3 111The label on the
$y$ axis in the right panel of Fig. 3 in Sotani & Sumiyoshi (2019) is typo,
which should be the stellar average density, $M_{\rm PNS}/R_{\rm PNS}^{3}$,
instead of the square root of it, $(M_{\rm PNS}/R_{\rm PNS}^{3})^{1/2}$, where
$T_{\rm BH}$ denotes the black hole formation time with $T_{\rm pb}$ and the
concrete values are listed in Table 2. In the current study, we define the
black hole formation time by finding the apparent horizon, i.e., the PNS
models adopted in this study cover the evolution up to the final collapse to
the black hole formation222The simulations have been performed even a little
after the formation of apparent horizon in most of simulations except for
T50-Shen.. The apparent horizon is formed in the snapshot at the final time
step except for the case of T50-Shen. We note that the values of $T_{\rm BH}$
shown in Table 2 are a little larger than those shown in Table I in Sotani &
Sumiyoshi (2019), where the final snapshot is just before the collapse to the
black hole. Toward the moment of the black hole formation, $T_{\rm BH}$, the
PNS radius suddenly decreases as shown in Figs. 1 and 2. From Table 2, one can
see that the black hole formation time is shorter for the PNS models with
softer EOS. In Fig. 2, the time variation of the PNS radius and the average
density change faster with softer EOS. In addition, we remark that the maximum
gravitational mass for PNSs is $\sim 10\%$ larger than that for cold neutron
stars, which comes from the thermal effect and the contribution of lepton,
such as an electron and neutrino.

Figure 2: Time evolution of the PNS radius (top), gravitational mass (middle), and average density (bottom) for the final phase before the gravitational collapse is shown as a function of $T_{\rm pb}-T_{\rm BH}$ for various models with $\rho_{\rm s}=2\times 10^{11}$ g/cm3.  Table 2: $T_{\rm BH}$ and $T_{f_{\xi}=0}$ for various models, where $T_{\rm BH}$ and $T_{f_{\xi}=0}$ are respectively the black hole formation time, i.e., the time when the apparent horizon appears inside the PNS, and the time when the fundamental frequency of the radial oscillation becomes zero. EOS | Model |  $T_{\rm BH}$ (sec) |  $T_{f_{\xi}=0}$ (sec)
---|---|---|---
Shen | W40 | 1.345 | 1.345
| T50 | 1.511 | 1.509
LS180 | W40 | 0.566 | 0.561
| T50 | 0.507 | 0.506
LS220 | W40 | 0.784 | 0.782

Finally, as in the previous studies (e.g., Sotani & Takiwaki (2016); Morozova
et al. (2018); Sotani et al. (2019); Torres-Forné et al. (2018, 2019a); Sotani
& Sumiyoshi (2019)), we assume that the PNS model are in a static equilibrium
at each time step and prepare the PNS model as a background model for
considering the linear analysis in this study. In this case, the metric is
given with the spherical coordinate as

$ds^{2}=-e^{2\Phi}dt^{2}+e^{2\Lambda}dr^{2}+r^{2}\left(d\theta^{2}+\sin^{2}d\phi^{2}\right).$
(2)

Now, $\Phi$ and $\Lambda$ are functions only of $r$ and $e^{2\Lambda}$ is
directly associated with the mass function $m(r)$ via $e^{-2\Lambda}=1-2m/r$.

##  3 Stability of Protoneutron stars

As seen in the previous section, the mass of the protoneutron star increases
with the accretion, which eventually approaches the maximum mass allowed with
the adopted EOSs. Then, the protoneutron star would gravitationally collapse
to a black hole. The moment when the protoneutron star approaches its maximum
mass corresponds to the onset of the instability. In order to determine the
onset of instability in the evolution of the protoneutron stars, we make a
linear analysis with the radial perturbation on the protoneutron star models
at each time step after core bounce. For this purpose, one can derive the
perturbation equations as

$\displaystyle\frac{d\xi}{dr}$ $\displaystyle=$
$\displaystyle-\left[\frac{3}{r}+\frac{p^{\prime}}{p+\varepsilon}\right]\xi-\frac{\eta}{r\Gamma},$
(3) $\displaystyle\frac{d\eta}{dr}$ $\displaystyle=$
$\displaystyle\left[r(p+\varepsilon)e^{2\Lambda}\left(\frac{\omega^{2}}{p}e^{-2\Phi}-8\pi\right)-\frac{4p^{\prime}}{p}+\frac{r(p^{\prime})^{2}}{p(p+\varepsilon)}\right]\xi-\left[\frac{\varepsilon
p^{\prime}}{p(p+\varepsilon)}+4\pi r(p+\varepsilon)e^{2\Lambda}\right]\eta,$
(4)

where $p$, $\varepsilon$, and $\Gamma$ denote the pressure, energy density,
and adiabatic index for the background protoneutron star models, while $\xi$
and $\eta$ are perturbative variables given by $\xi\equiv\Delta r/r$ and
$\eta\equiv\Delta p/p$ with the radial displacement, $\Delta r$, and the
Lagrangian perturbation of pressure, $\Delta p$ (Chandrasekhar, 1964;
Chanmugan, 1977; Gondek, Haensel & Zdunik, 1997) 333The perturbation equations
for radial oscillations considered in this study are derived without the
Cowling approximation, while it is assumed only in the calculation of the non-
radial oscillations discussed in the section 4. In fact, if one assumes the
Cowling approximation for the radial oscillations, one can not discuss the
onset of a black hole formation (Torres-Forné et al., 2018).. The prime in the
equations denotes the radial derivative and the adiabatic index is given by

$\Gamma\equiv\left(\frac{\partial\ln p}{\partial\ln n_{\rm
b}}\right)_{s}=\frac{p+\varepsilon}{p}c_{s}^{2},$ (5)

where $n_{\rm b}$, $s$, and $c_{s}$ denote the baryon number density, entropy
per baryon, and sound velocity, respectively. We remark that one can derive
the Sturm-Liouville type second order differential equation with respect to
$\xi$ from Eqs. (3) and (4). To solve the eigenvalue problem with respect to
the eigenvalue $\omega^{2}$, one should impose the appropriate boundary
conditions. The boundary condition at the stellar surface comes from the
condition to remove the singularity in Eq. (4) (Chanmugan, 1977), i.e.,

$\eta=-\left[\left(\frac{\omega^{2}R_{\rm PNS}^{3}}{M_{\rm PNS}}+\frac{M_{\rm
PNS}}{R_{\rm PNS}}\right)\left(1-\frac{2M_{\rm PNS}}{R_{\rm
PNS}}\right)^{-1}+4\right]\xi,$ (6)

while the boundary condition at the center is the regularity condition, i.e.,

$3\Gamma\xi+\eta=0.$ (7)

In addition, as normalization of the eigenfunction, we set $\xi=1$ at the
stellar center. Then, with the resultant eigenvalue $\omega^{2}$, the
frequency of radial oscillations are given by

$f_{\xi}={\rm sgn}(\omega^{2})\sqrt{|\omega^{2}|}/2\pi,$ (8)

where the system is unstable when $\omega^{2}$ becomes negative.

For the case of cold neutron stars, it is well known that the stellar models
constructed with the central density, which is more than that for the neutron
star model with the maximum mass, are unstable. As an example (and as a test
of our code for the eigenvalue problem with respect to the radial
oscillations), we show the results for the neutron star models constructed
with Togashi EOS (Togashi et al., 2017) in Fig. 3, where top panels correspond
to the equilibrium models, i.e., the stellar mass is shown as a function of
the central density, $\rho_{c}$, normalized by the nuclear saturation density,
$\rho_{0}$, in the left-top panel and as a function of the stellar radius in
the right-top panel. The open-circle denotes the stellar model with the
maximum mass. In the bottom-left panel, we show the frequencies of the lowest
three radial oscillations as a function of $\rho_{c}/\rho_{0}$. From this
figure, it is obvious that the frequency of the lowest radial oscillation
becomes negative for the stellar models with the central density, which is
larger than that for the stellar model with the maximum mass, i.e., the
corresponding stellar models are unstable.

Figure 3: For the neutron stars constructed with Togashi EOS (Togashi et al.,
2017), the mass of stars is shown as a function of the central density,
$\rho_{c}$, normalized by the nuclear saturation density, $\rho_{0}$, in the
left-top panel and as a function of the stellar radius in the right-top panel,
while the corresponding frequencies of the lowest three radial oscillations
are shown as a function of $\rho_{c}/\rho_{0}$ in left-bottom panel. The open-
circle denotes the stellar model with the maximum mass. It is obvious that the
frequency of the lowest (1st) radial oscillation becomes negative for the
stellar model constructed with the central density, which is more than that
for the neutron star with the maximum mass. We remark that the central density
in this figure is not a rest-mass density, but an energy density.

Now, we consider the case of PNSs. In Fig. 4, we show the frequencies of the
lowest three radial oscillations on the various PNS models for the final phase
just before the black hole formation. From this figure, one can obviously see
that, for any PNS models considered in this study, the lowest frequency is
basically positive, i.e., PNS is stable, but it becomes negative at the last
moment. That is, the PNS model for which the lowest frequency is zero is the
PNS model with the maximum mass allowed with the adopted EOS, where the
corresponding time is shown in Table 2. Since one can see that $T_{\rm
BH}>T_{f_{\xi}=0}$ from Table 2 and Fig. 4, we find that the PNS becomes
unstable before the apparent horizon appears inside the PNS. In addition, we
find that the lowest frequency monotonically decreases with time at least in
the phase considered in this study, while the second and third lowest
frequencies for the models with soft EOSs increase. Incidentally, the second
and third lowest frequencies become negative at the end as for the model of
T50-Shen.

Furthermore, we focus on the lowest frequency in Fig. 5, where the lowest
frequency is shown as a function of $T_{\rm pb}-T_{\rm BH}$ (left panel) and
the PNS compactness (right panel). From the left panel, we see that the time
evolution of the lowest frequency (especially in the last $\sim 40$ msec)
weakly depend on the PNS model, although the second and third lowest
frequencies depend on the PNS model, as mentioned above. We found that, from
the right panel, the gradient of the lowest frequency with respect to the PNS
compactness seems to be almost independent of the PNS models at the last
moment, although as shown in Fig. 2 the PNS properties in the final phase
depend on PNS models. In practice, the radial oscillation itself is not
associated with the gravitational radiations, because radial oscillations of a
spherical background do not induce a quadrupolar deformation and hence no
gravitational waves, but one can observe the aspect of the radial oscillations
in the gravitational waves through the nonlinear coupling between the radial
and nonradial oscillations (Passamonti, Stergioulas & Nagar, 2007). On the
other hand, apart from non-linear coupling, deformations due to stellar
rotation can also make the radial oscillations visible (Cerdá-Durán et al.,
2013; Torres-Forné et al., 2019b). In presence of rotation, the radial
oscillations become quasi-radial, however, for sufficiently slow rotation
their frequencies should be close to the radial modes computed here.

Figure 4: The frequencies of the lowest three radial oscillations are shown as
a function of $T_{\rm pb}-T_{\rm BH}$, where the top and bottom panels
correspond to the PNS models with W40 and T50, while the left, middle and
right panels correspond to the PNS models with Shen, LS180, and LS220,
respectively.  Figure 5: For various PNS models, the frequencies of the 1st
radial oscillations for the final phase are shown as a function of $T_{\rm
pb}-T_{\rm BH}$ (left panel) and $M_{\rm PNS}/R_{\rm PNS}$ (right panel).

##  4 Protoneutron star asteroseismology

In this study, we also examine the gravitational wave frequency from the PNSs
in the final phase just before the black hole formation. For this purpose, as
in Sotani & Sumiyoshi (2019) we simply adopt the relativistic Cowling
approximation, i.e, the metric perturbations are neglected during the fluid
oscillations, where one can derive the perturbation equations by linearizing
the energy-momentum conservation law. Then, by adopting the appropriate
boundary conditions at the stellar center and the radius, the problem to solve
becomes an eigenvalue problem with respect to the eigenvalue $\omega$, with
which the gravitational wave frequency, $f$, is determined via
$f=\omega/(2\pi)$. The perturbation equations and boundary conditions are the
same as in Sotani et al. (2019). That is, we make a linear analysis in a
similar way to the previous study (Sotani & Sumiyoshi, 2019), but we
especially focus on the final phase toward the black hole formation in this
study. In this study, we only consider the $\ell=2$ mode oscillations, which
are considered to be energetically dominant in gravitational wave radiation.

First, we examine how the gravitational wave frequencies from the PNSs depend
on the PNS surface density. This dependence has already been discussed in
Morozova et al. (2018); Sotani & Takiwaki (2020b) for the PNSs produced by the
successful supernova explosions, where the $f$\- and $g_{1}$-mode frequencies
are independent of the selection of the surface density after the avoided
crossing between the both modes, while the other modes, especially the
$p_{i}$-modes with $i\geq 2$, depend on the surface density. This behavior may
be understood by considering the radial profile of pulsation energy density
(or eigenfunction) (Sotani & Takiwaki, 2020b). That is, the $f$\- and
$g_{1}$-modes significantly oscillate inside the star, while the other modes
oscillate not only inside the star but also in the region closed to the PNS
surface.

In the case for the PNSs in the failed supernovae considered in this study, we
show the gravitational wave frequencies of the $f$-, $p_{i}$-, and
$g_{i}$-modes with $i=1-3$ in Fig. 6 for the PNS models of W40-Shen (left
panel) and W40-LS220 (right panel) with $\rho_{s}=10^{11}$ (open-marks) and
$2\times 10^{11}$ g/cm3 (filled mark). From this figure, we confirm that the
$f$\- and $g_{1}$-mode frequencies are basically independent of the selection
of the surface density even for the failed supernovae after the avoided
crossing between the $f$\- and $g_{1}$-modes. We also find that the dependence
of the surface density can be seen in the $g_{1}$-mode frequency in the final
phase just before the black hole formation, which is the period after the
avoided crossing between the $g_{1}$\- and $g_{2}$-modes. In addition, as in
the previous studies for the PNSs produced via the successful supernova
explosions, the other modes seem to depend on the selection of surface
density, where we find that the $p_{1}$-mode frequency significantly depends
on the surface density in the late phase after the core bounce. We remark that
a wiggling in the $p_{1}$-mode frequency with $\rho_{s}=10^{11}$ g/cm3 may
come from the fact that the numerical resolution is not so high in a lower
density region.

Figure 6: Comparison between the eigenfrequencies of the $f$-, $p_{i}$-, and
$g_{i}$-modes with $i=1-3$ for the PNSs models with $\rho_{s}=2.0\times
10^{11}$ g/cm3 (filled-mark with dashed line) and those with
$\rho_{s}=10^{11}$ g/cm3 (open-mark with dashed line) for W40-Shen (left
panel) and W40-LS220 (right panel). Although we do not show, the additional
eigenfrequencies exist above the $p_{3}$\- and below the $g_{3}$-modes, which
are respectively the $p_{i}$\- and the $g_{i}$-modes with $i\geq 4$.  Figure
7: Time evolution of gravitational wave frequencies from the PNS model of
W40-Shen with $\rho_{s}=2\times 10^{11}$ g/cm3, where the $f$-, $p_{i}$-, and
$g_{i}$-modes with $i=1-5$ are shown up to the black hole formation. The right
panel is just an enlarged view of the left panel. In order to focus on the PNS
models in the final phase just before the black hole formation. Note that we
adopt a small time interval for the snapshots for the final phase.

Next, we focus on the gravitational wave frequencies from the PNS in the final
phase just before the black hole formation. In Fig. 7 we show the
gravitational wave frequencies from the PNS model of W40-Shen with
$\rho_{s}=2\times 10^{11}$ g/cm3. Note that we adopt a small time interval for
the PNS models in the final phase to study in detail the black hole formation.
From this figure, one can observe that the gravitational wave frequencies
steeply increase at the final moment when the PNS collapses to a black hole.
This behavior can be seen not only in the PNS model of W40-Shen but also in
the other PNS models as shown in appendix A. Thus, this steep increase in the
gravitational wave frequencies would tell us the moment of the black hole
formation. This feature is revealed in the current study, more clearly than
Sotani & Sumiyoshi (2019) where the analysis is made only before the final
collapse.

We explore the empirical formulae for the evolution of frequency of
gravitational waves. In Sotani & Sumiyoshi (2019), we just derive a fitting
formula for the $f$-mode frequency as a function of the square root of the PNS
average density according to Andersson & Kokkotas (1996, 1998), because the
$f$-mode oscillation is a kind of an acoustic oscillation, which could be
associated with the density of the object. On the other hand, Tsui & Leung
(2005); Sotani (2021) shows another possibility for fitting of the $f$-mode
frequency, which seems to be better than that in Andersson & Kokkotas (1996,
1998). In fact, as in Tsui & Leung (2005); Sotani (2021), we find that the
$f$-mode frequency multiplied by the PNS mass, $f_{f}(M/1.4M_{\odot})$, can be
expressed as a function of the PNS compactness, ${\cal C}\equiv M_{\rm
PNS}/R_{\rm PNS}$, almost independently of the PNS models even for the last
PNS model in the evolution considered in this study, such as

$f_{f}\left(\frac{M_{\rm
PNS}}{1.4M_{\odot}}\right){\rm(kHz)}=0.08087+3.262\left(\frac{\cal
C}{0.2}\right)-1.740\left(\frac{\cal C}{0.2}\right)^{2}+2.204\left(\frac{\cal
C}{0.2}\right)^{3}.$ (9)

In Fig. 8 we show $f_{f}(M/1.4M_{\odot})$ as a function of ${\cal C}$ for
various PNS models, where the thick-solid line denotes the fitting formula
given by Eq. (9). In this figure, the PNSs evolve from the bottom left to the
top right. With this empirical formula, one can obtain a constraint on the
relation between the PNS mass and radius at each time step by observing the
$f$-mode gravitational wave from the PNS 444The similar universal relation has
been also derived in Torres-Forné et al. (2019b), where the mode
classification is different from ours because the position of outer boundary
and the boundary condition are different from us. Additionally, in Bizouard et
al. (2021) it is also shown that PNS asteroseismology using gravitational
waves is possible in practice.. We remark that what is called here $f$-mode is
also interpreted as a surface $g$-mode (Müller, Janka & Marek, 2013; Cerdá-
Durán et al., 2013) or a $g$-mode like oscillations by adopting the outer
boundary at the shock radius (Torres-Forné et al., 2019a, b), where the
corresponding frequencies can be expressed well with the PNS surface gravity,
$M_{\rm PNS}/R_{\rm PNS}^{2}$.

Figure 8: The $f$-mode frequency multiplied by the normalized PNS mass is show
as a function of the PNS compactness for various PNS models up to the black
hole formation. The thick-solid line denote the fitting formula given by Eq.
(9).

In addition, it is pointed out that the ratio of the $p_{1}$-mode frequency to
the $f$-mode frequency, $f_{p_{1}}/f_{f}$, can be characterized by a new
variable, $Q_{0}$, defined in Camelio et al. (2017) as

$Q_{0}=\frac{1}{f_{f}}\left(\frac{M_{\rm
PNS}}{1.4M_{\odot}}\right)^{1/2}\left(\frac{R_{\rm PNS}}{10\ {\rm
km}}\right)^{-3/2}.$ (10)

In Sotani & Sumiyoshi (2019), we could not find any good relation between
$f_{p_{1}}/f_{f}$ and $Q_{0}$ 555Explicitly speaking, we could not find any
relation between $f_{f}/f_{p_{1}}$ (instead of $f_{p_{1}}/f_{f}$) and $Q_{0}$
in Sotani & Sumiyoshi (2019)., where the PNS surface density is chosen to be
$\rho_{s}=10^{11}$ g/cm3. On the other hand, in this study with
$\rho_{s}=2\times 10^{11}$ g/cm3, we can get a good relation between
$f_{p_{1}}/f_{f}$ and $Q_{0}$. In Fig. 9, $f_{p_{1}}/f_{f}$ is plotted as a
function of $Q_{0}$ for various PNS models. In this figure, the PNSs evolve
from left to right for the early phase, but their evolutionary direction
changes from right to left for the later phase. From this figure, one can
observe that $f_{p_{1}}/f_{f}$ is almost independent of the PNS models for
$Q_{0}\,\,\raisebox{-3.44444pt}{$\stackrel{{\scriptstyle\textstyle<}}{{\sim}}$
}0.2$ kHz-1 and weakly depends on the PNS models for
$Q_{0}\,\raisebox{-3.44444pt}{$\stackrel{{\scriptstyle\textstyle>}}{{\sim}}$
}0.2$ kHz-1. With using the data for
$Q_{0}\,\raisebox{-3.44444pt}{$\stackrel{{\scriptstyle\textstyle>}}{{\sim}}$
}0.13$ kHz-1, we can obtain the fitting formula such as

$f_{p_{1}}/f_{f}=2.845-22.608\left(\frac{Q_{0}}{1{\rm/kHz}}\right)+76.226\left(\frac{Q_{0}}{1{\rm/kHz}}\right)^{2},$
(11)

with which the expected value of $f_{p_{1}}/f_{f}$ is plotted with the thick-
solid line in Fig. 9. As mentioned above, since the $p_{1}$-mode frequency
depends on the selection of the PNS surface density, it is still uncertain how
universal this fitting formula is. But, if this type of fitting formula truly
exists, the simultaneous observation of the $f$\- and $p_{1}$-mode
gravitational waves tells us another constraint on the relation between the
PNS mass and radius at each time step. The dependence of the $p_{1}$-mode
frequency on the PNS surface density comes from the fact that the pulsation
energy density for the $p_{1}$-mode oscillates not only inside the PNS but
also in the vicinity of the PNS surface (Sotani & Takiwaki, 2020b). To
determine the $p_{1}$-mode frequency with a realistic situation, one may have
to impose the outer boundary condition at the shock radius, instead of the PNS
surface, as in Torres-Forné et al. (2018, 2019a, 2019b).

Figure 9: For various PNS models up to the black hole formation, the ratio of
the $p_{1}$-mode frequency to the $f$-mode frequency, $f_{p_{1}}/f_{f}$, is
shown as a function of $Q_{0}$ defined by Eq. (10). The thick-solid line
denotes the fitting formula given by Eq. (11), which is derived with using the
data for
$Q_{0}\,\raisebox{-3.44444pt}{$\stackrel{{\scriptstyle\textstyle>}}{{\sim}}$
}0.13$ kHz-1.

Finally, we also check the dependence of the ratio of the $g_{1}$-mode
frequency to the $f$-mode frequency, $f_{g_{1}}/f_{f}$, on the PNS
compactness. In Sotani & Sumiyoshi (2019), it is already shown that
$f_{g_{1}}/f_{f}$ can be characterized by the PNS compactness almost
independently of the PNS models. We found that this universality with the
frequencies calculated in this study holds up to the final moment of black
hole formation, as shown in Fig. 10. Again, the simultaneous observation of
the $f$\- and $g_{1}$-mode gravitational waves tells us a constraint on the
PNS compactness at each time step.

Figure 10: For various PNS models up to the black hole formation, the ratio of
the $g_{1}$-mode frequency to the $f$-mode frequency, $f_{g_{1}}/f_{f}$, is
shown as a function of the PNS compactness.

##  5 Conclusion

We examine the PNS stability and gravitational wave frequencies from the
massive PNS produced by failed supernovae. In order to prepare the PNS model,
we adopt the numerical data obtained by the spherically symmetric numerical
simulation performed by solving the general relativistic neutrino-radiation
hydrodynamics. Then, as in the previous studies, we assume that the PNS models
are static equilibrium at each time step. On such a PNS model, first we
examine the PNS stability by solving the radial oscillations. We show that
most of PNS models in their evolution are stable against the radial
perturbations, i.e., all the frequencies of the radial oscillations are
positive, while the PNS model finally becomes unstable before the apparent
horizon appears inside the PNS. We also show that the time evolution of the
lowest frequency of the radial oscillation for the last $\sim 40$ msec weakly
depend on the PNS model. In addition, we examine the gravitational wave
frequencies from the PNS models by adopting the relativistic Cowling
approximation. In this study we especially focus on the final phase in this
study, although how to analysis is the same as in the previous study. We find
that the frequencies steeply increase at the final moment when a black hole
forms. So, such a signal in gravitational wave frequencies would tell us the
moment of black hole formation. We also show that the $f$-mode frequency
multiplied by the PNS mass can be characterized by the PNS compactness almost
independently of the PNS models and derive the empirical formula for the
$f$-mode frequency. Moreover, we find that the ratio of the $p_{1}$-mode to
the $f$-mode weakly can be characterized by the variable, which is defined as
the square root of the normalized PNS average density divided by the $f$-mode
frequency, while the ratio of the $g_{1}$-mode to the $f$-mode can be
characterized by the PNS compactness. So, via the simultaneous observation of
the $f$\- and $p_{1}$-mode (or the $f$\- and $g_{1}$-mode) gravitational wave,
one would constrain the PNS mass and radius at each time step, which helps us
to understand the dense matter physics. In this study, as a first step, we
simply assume that the PNS models are static equilibrium at each time step
even for the final phase, which could be highly dynamical phase. Since this
static assumption may not be valid on such a dynamical phase, we will study
the PNS oscillations without the static assumption somewhere in the future. In
fact, by comparing the dynamical timescale of the system defined as $\tau_{\rm
dyn}=(R_{\rm PNS}^{3}/M_{\rm PNS})^{1/2}$ with the typical timescale of the
$f$-mode oscillation, i.e., $1/f_{f}$, we find that the value of
$f_{f}\tau_{\rm dyn}$ becomes $\sim 0.3$ in the final dynamical phase.

## Acknowledgements

This work is supported in part by Japan Society for the Promotion of Science
(JSPS) KAKENHI Grant Numbers JP17H06357, JP17H06365, JP18H05236, JP19KK0354,
JP19K03837, JP20H01905, JP20H04753, and JP21H01088 and by Pioneering Program
of RIKEN for Evolution of Matter in the Universe (r-EMU). For providing high
performance computing resources, Computing Research Center, KEK, JLDG on
SINET4 of NII, Research Center for Nuclear Physics, Osaka University, Yukawa
Institute of Theoretical Physics, Kyoto University, Nagoya University,
Hokkaido University, and Information Technology Center, University of Tokyo
are acknowledged. This work was partly supported by MEXT as “Program for
Promoting Researches on the Supercomputer Fugaku” (Toward a unified view of
the universe: from large scale structures to planets) and the Particle,
Nuclear and Astro Physics Simulation Program (Nos. 2020-004, 2021-004) of
Institute of Particle and Nuclear Studies, High Energy Accelerator Research
Organization (KEK).

## Data availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

##  Appendix A Time evolution of gravitational wave frequencies from various
PNS models

We discuss the evolution of gravitational wave frequencies with the PNS model
of W40-Shen in main text, as shown in Fig. 7. In this appendix, we show the
evolution of gravitational wave frequencies from the other PNS models in Fig.
11. We remark that the top-left panel is the same as the left panel in Fig. 7.
We cover the time periods up to the black hole formation (the formation of
apparent horizon) with a small time interval for the final phase.

Figure 11: Time evolution of gravitational wave frequencies from the various
PNS models. The left and right panels correspond to the results with W40 and
T50, while the top, middle, and bottom panels correspond to the results from
the PNS constructed with Shen, LS180, and LS220, respectively. In all panels,
the $f$-, $g_{i}$-, and $p_{i}$-modes with $i=1-5$ are shown. Note that we
adopt a small time interval for the snapshots for the final phase.

## References

  * Aasi et al. (2015) Aasi J. et al., 2015, Classical and Quantum Gravity, 32, 074001
  * Acernese et al. (2014) Acernese F. et al., 2014, Classical and Quantum Gravity, 32, 024001
  * Abbott et al. (2017a) Abbott B. P. et al. (LIGO Scientific and Virgo Collaboration), 2017a, Phys. Rev. Lett., 119, 161101
  * Abbott et al. (2017b) Abbott B. P. et al. (LIGO Scientific and Virgo Collaboration), 2017b, ApJ, 848, L12
  * Abbott et al. (2017c) Abbott B. P. et al. (LIGO Scientific and Virgo Collaboration), 2017c, Class. Quantum Grav., 34, 044001
  * Andresen et al. (2017) Andresen H., Müller B., Müller E., Janka H.-T., 2017, MNRAS, 468, 2032
  * Andresen et al. (2019) Andresen H., Müller E., Janka H. -Th., Summa A., Gill K., Zanolin M., 2019, MNRAS, 486, 2238
  * Andersson & Kokkotas (1996) Andersson N., Kokkotas K. D., 1996, Phys. Rev. Lett., 77, 4134
  * Andersson & Kokkotas (1998) Andersson N., Kokkotas K. D., 1998, MNRAS, 299, 1059
  * Akutsu et al. (2020) Akutsu T. et al. 2020, Prog. Theor. Exp. Phys., 2020, 05A103
  * Antoniadis et al. (2013) Antoniadis J. et al., 2013, Science, 340, 1233232
  * Bizouard et al. (2021) Bizouard M. A. et al. 2021, Phys. Rev. D, 103, 063006
  * Burgio et al. (2011) Burgio G. F., Ferrari V., Gualtieri L., Schulze H.-J., 2011, Phys. Rev. D, 84, 044017
  * Bruenn (1985) Bruenn S. W., 1985, ApJS, 58, 771
  * Camelio et al. (2017) Camelio G., Lovato A., Gualtieri L., Benhar O., Pons J. A., Ferrari V., 2017, Phys. Rev. D, 96, 043015
  * Cerdá-Durán et al. (2013) Cerdá-Durán P., DeBrye N., Aloy M. A., Font J. A., Obergaulinger M., 2013, ApJ, 779, L18
  * Chandrasekhar (1964) Chandrasekhar S., 1964, ApJ, 140, 417
  * Chanmugan (1977) Chanmugan G., 1977, ApJ, 217, 799
  * Cromartie et al. (2020) Cromartie H. T. et al., 2020, Nature Astronomy, 4, 72
  * Demorest et al. (2010) Demorest P., Pennucci T., Ransom S., Roberts M., Hessels J., 2010, Nature, 467, 1081
  * Doneva et al. (2013) Doneva D. D., Gaertig E., Kokkotas K. D., Krüger C., 2013, Phys. Rev. D, 88, 044052
  * Ferrari, Miniutti & Pons (2003) Ferrari V., Miniutti G., Pons J. A., 2003, MNRAS, 342, 629
  * Fuller et al. (2015) Fuller J., Klion H., Abdikamalov E., Ott C. D., 2015, MNRAS, 450, 414
  * Gearheart et al. (2011) Gearheart M., Newton W. G., Hooker J., Li B. A., 2011, MNRAS, 418, 2343
  * Gondek, Haensel & Zdunik (1997) Gondek D., Haensel P., Zdunik J. L., 1997, Astron. Astrophys., 325, 217
  * Horiuchi et al. (2018) Horiuchi S., Sumiyoshi K., Nakamura K., Fischer T., Summa A., Takiwaki T., Janka H.-T., Kotake K., 2018, MNRAS, 475, 1363
  * Khan & Margueron (2003) Khan E., Margueron J., 2013, Phys. Rev. C, 88, 034319
  * Kokkotas & Ruoff (2001) Kokkotas K. D., Ruoff J., 2001, Astron. Astrophys., 366, 565
  * Kuroda, Kotake & Takiwaki (2016) Kuroda T., Kotake K., Takiwaki T., 2016, ApJ, 829, L14
  * Kuroda et al. (2017) Kuroda T., Kotake K., Hayama K., Takiwaki T., 2017, ApJ, 851, 62
  * Kuroda et al. (2018) Kuroda T., Kotake K., Takiwaki T., Thielemann F.-K., 2018, MNRAS, 477, L80
  * Lattimer & Swesty (1991) Lattimer J. M., Swesty F. D., 1991, Nucl. Phys. A 535, 331
  * Li et al. (2019) Li B. A., Krastev P. G., Wen D. H., Zhang N. B., 2019, Eur. Phys. J. A 55, 117
  * Mezzacappa et al. (2020) Mezzacappa A., Marronetti P., Landfield R. E., Lentz E. J., Yakunin K. N., Bruenn S. W., Hix W. R., Messer O. E. B., Endeve E., Blondin J. M., Harris J. A., 2020, Phys. Rev. D, 102, 023027
  * Misner & Sharp (1964) Misner C. W., Sharp D. H., 1964, Phys. Rev. B, 136, 571
  * Morozova et al. (2018) Morozova V., Radice D., Burrows A., Vartanyan D., 2018, ApJ, 861, 10
  * Müller, Janka & Marek (2013) Müller B., Janka H. -T., Marek A., 2013, ApJ, 766, 43
  * Murphy, Ott & Burrows (2009) Murphy J. W., Ott C. D., Burrows A., 2009, ApJ, 707, 1173
  * Nakazato et al. (2010) Nakazato K., Sumiyoshi K., Suzuki H., Yamada S., 2010, Phys. Rev. D, 81, 083009
  * Nakazato et al. (2013) Nakazato K., Sumiyoshi K., Suzuki H., Totani T., Umeda H., Yamada S., 2013, ApJS, 205, 2
  * O’Connor & Ott (2011) O’Connor E., Ott C. D., 2011, ApJ, 730, 70
  * O’Connor & Couch (2018) O’Connor E. P., Couch S. M., 2018, ApJ, 865, 81
  * Oertel et al. (2017) Oertel M., Hempel M., Klähn T., Typel S., 2017, Rev. Mod. Phys., 89, 015007
  * Ott et al. (2013) Ott, C. D. Abdikamalov E., Mösta P., Haas R., Drasco S., O’Connor E. P., Reisswig C., Meakin C. A., Schnetter E., 2013, ApJ 768, 115
  * Pan et al. (2018) Pan K.-C., Liebendörfer M., Couch S. M., Thielemann F.-K., 2018, ApJ, 857, 13
  * Passamonti, Stergioulas & Nagar (2007) Passamonti A., Stergioulas N., Nagar A., 2007, Phys. Rev. D, 75, 084038
  * Passamonti & Andersson (2012) Passamonti A., Andersson N., 2012, MNRAS, 419, 638
  * Punturo, Lück & Beker (2014) Punturo M., Lück H., Beker M., A Third Generation Gravitational Wave Observatory: The Einstein Telescope, edited by M. Bassan, Advanced Interferometers and the Search for Gravitational Waves. Astrophysics and Space Science Library Vol. 404 (Springer, Cham, 2014)
  * Radice et al. (2019) Radice D., Morozova V., Burrows A., Vartanyan D., Nagakura H., 2019, ApJ, 876, L9
  * Reed et al. (2021) Reed B. T., Fattoyev F. J., Horowitz C. J., Piekarewicz J., 2021, Phys. Rev. Lett., 126, 172503
  * Shen et al. (1998) Shen H., Toki H., Oyamatsu K., Sumiyoshi K., 1998, Nucl. Phys. A, 637, 435
  * Sotani, Tominaga & Maeda (2001) Sotani H., Tominaga K., Maeda K. I., 2001, Phys. Rev. D, 65, 024010
  * Sotani, Kohri & Harada (2004) Sotani H., Kohri K., Harada T., 2004, Phys. Rev. D, 69, 084008
  * Sotani et al. (2011) Sotani H., Yasutake N., Maruyama T., Tatsumi T., 2011, Phys. Rev. D, 83, 024014
  * Sotani et al. (2012) Sotani H., Nakazato K., Iida K., Oyamatsu K., 2012, Phys. Rev. Lett., 108, 201101
  * Sotani et al. (2013a) Sotani H., Nakazato K., Iida K., Oyamatsu K., 2013a, MNRAS, 428, L21
  * Sotani et al. (2013b) Sotani H., Nakazato K., Iida K., Oyamatsu K., 2013b, MNRAS, 434, 2060
  * Sotani, Iida & Oyamatsu (2016) Sotani H., Iida K., Oyamatsu K., 2016, New Astron., 43, 80
  * Sotani & Takiwaki (2016) Sotani H., Takiwaki T., 2016, Phys. Rev. D, 94, 044043
  * Sotani, Iida & Oyamatsu (2017) Sotani H., Iida K., Oyamatsu K., 2017, MNRAS, 464, 3101
  * Sotani et al. (2017) Sotani H., Kuroda T., Takiwaki T., Kotake K., 2017, Phys. Rev. D, 96, 063005
  * Sotani, Iida & Oyamatsu (2018) Sotani H., Iida K., Oyamatsu K., 2018, MNRAS, 479, 4735
  * Sotani, Iida & Oyamatsu (2019) Sotani H., Iida K., Oyamatsu K., 2019, MNRAS, 489, 3022
  * Sotani et al. (2019) Sotani H., Kuroda T., Takiwaki T., Kotake K., 2019, Phys. Rev. D, 99, 123024
  * Sotani & Sumiyoshi (2019) Sotani H., Sumiyoshi K., 2019, Phys. Rev. D, 100, 083008
  * Sotani & Takiwaki (2020a) Sotani H., Takiwaki T., 2020a, Phys. Rev. D, 102, 023028
  * Sotani & Takiwaki (2020b) Sotani H., Takiwaki T., 2020b, MNRAS, 498, 3503
  * Sotani & Takiwaki (2020c) Sotani H., Takiwaki T., 2020c, Phys. Rev. D, 102, 063025
  * Sotani (2020a) Sotani H., 2020a, Phys. Rev. D, 102, 063023
  * Sotani (2020b) Sotani H., 2020b, Phys. Rev. D, 102, 103021
  * Sotani (2021) Sotani H., 2021, Phys. Rev. D, 103, 123015
  * Sumiyoshi et al. (2004) Sumiyoshi K., Suzuki H., Yamada S., Toki H., 2004, Nucl. Phys. A, 730, 227
  * Sumiyoshi et al. (2005) Sumiyoshi K., Yamada S., Suzuki H., Shen H., Chiba S., Toki H., 2005, ApJ, 629, 922
  * Sumiyoshi et al. (2006) Sumiyoshi K., Yamada S., Suzuki H., Chiba S., 2006, Phys. Rev. Lett., 97, 091101
  * Sumiyoshi, Yamama & Suzuki (2007) Sumiyoshi K., Yamada S., Suzuki H., 2007, ApJ, 667, 382
  * Sumiyoshi, Yamama & Suzuki (2008) Sumiyoshi K., Yamada S., Suzuki H., 2008, ApJ, 688, 1176
  * Sumiyoshi et al. (2009) Sumiyoshi K., Ishizuka C., Ohnishi A., Yamada S., Suzuki H., 2009, ApJ, 690, L43
  * Sukhbold et al. (2016) Sukhbold T., Ertl T., Woosley S. E., Brown J. M., Janka H.-T., 2016, ApJ, 821, 38
  * Takiwaki & Kotake (2018) Takiwaki T., Kotake K., 2018, MNRAS, 475, L91
  * Togashi et al. (2017) Togashi H., Nakazato K., Takehara Y., Yamamuro S., Suzuki H., Takano M., 2017, Nucl. Phys. A, 961, 78
  * Tominaga, Umeda & Nomoto (2007) Tominaga N., Umeda H., Nomoto K., 2007, ApJ, 660, 516
  * 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. (2019a) Torres-Forné A., Cerdá-Durán P., Passamonti A., Obergaulinger M., Font J. A., 2019a, MNRAS, 482, 3967
  * Torres-Forné et al. (2019b) Torres-Forné A., Cerdá-Durán P., Obergaulinger M., Müller B., Font J. A., 2019, Phys. Rev. Letts., 123, 051102
  * Tsui & Leung (2005) Tsui L. K., Leung P. T., 2005, MNRAS, 357, 1029
  * Westernacher-Schneider et al. (2019) Westernacher-Schneider J. R., O’Connor E., O’Sullivan E., Tamborra I., Wu M.-R., Couch S. M., Malmenbeck F., 2019, Phys. Rev. D, 100, 123009
  * Woosley & Weaver (1995) Woosley S. E., Weaver T. A., 1995, ApJ. Suppl., 101, 181
  * Yakunin et al. (2015) Yakunin K. N., Mezzacappa A., Marronetti P., Yoshida S., Bruenn S. W., Hix W. R., Lentz E. J., Bronson Messer O. E., Harris J. A., Endeve E., Blondin J. M., Lingerfelt E. J., 2015, Phys. Rev. D, 92, 084040
  * Yamada (1997) Yamada S., 1997, ApJ, 475, 720
  * Yamada, Janka & Suzuki (1999) Yamada S., Janka H.-T., Suzuki H., 1999, Astron. Astro., 344, 533
  * Yokozawa et al. (2015) Yokozawa T., Asano M., Kayano T., Suwa Y., Kanda N., Koshio Y., Vagins M. R., 2015, ApJ, 811, 86