# A New Kilohertz Gravitational-Wave Feature from Rapidly Rotating Core-
Collapse Supernovae

He-Feng Hsieh Physics Division, National Center for Theoretical Sciences,
National Taiwan University, Taipei 10617, Taiwan  Department of Physics,
National Taiwan University, Taipei 10617, Taiwan  Institute of Astronomy,
National Tsing Hua University, Hsinchu 30013, Taiwan  Center for Informatics
and Computation in Astronomy, National Tsing Hua University, Hsinchu 30013,
Taiwan  Rubén Cabezón Center for Scientific Computing - sciCORE, Universität
Basel, Klingelbergstrasse 82, CH-4056 Basel, Switzerland  Li-Ting Ma (馬麗婷)
Institute of Astronomy, National Tsing Hua University, Hsinchu 30013, Taiwan
Center for Informatics and Computation in Astronomy, National Tsing Hua
University, Hsinchu 30013, Taiwan  Department of Physics, National Tsing Hua
University, Hsinchu 30013, Taiwan  Kuo-Chuan Pan (潘國全) Physics Division,
National Center for Theoretical Sciences, National Taiwan University, Taipei
10617, Taiwan  Institute of Astronomy, National Tsing Hua University, Hsinchu
30013, Taiwan  Center for Informatics and Computation in Astronomy, National
Tsing Hua University, Hsinchu 30013, Taiwan  Department of Physics, National
Tsing Hua University, Hsinchu 30013, Taiwan  Center for Theory and
Computation, National Tsing Hua University, Hsinchu 30013, Taiwan

###### Abstract

We present self-consistent three-dimensional core-collapse supernova
simulations of a rotating $20M_{\odot}$ progenitor model with various initial
angular velocities from $0.0$ to $4.0$ rad s-1 using a smoothed particle
hydrodynamics code, SPHYNX, and a grid-based hydrodynamics code, FLASH. We
identify two strong gravitational-wave features, with peak frequencies of
$\sim 300$ Hz and $\sim 1.3$ kHz in the first $100$ ms postbounce. We
demonstrate that these two features are associated with the $m=1$ deformation
from the proto-neutron star (PNS) modulation induced by the low-$T/|W|$
instability, regardless of the simulation code. The $300$ Hz feature is
present in models with an initial angular velocity between $1.0$ and $4.0$ rad
s-1, while the $1.3$ kHz feature is present only in a narrower range, from
$1.5$ to $3.5$ rad s-1. We show that the $1.3$ kHz signal originates from the
high-density inner core of the PNS, and the $m=1$ deformation triggers a
strong asymmetric distribution of electron anti-neutrinos. In addition to the
$300$ Hz and $1.3$ kHz features, we also observe one weaker but noticeable
gravitational-wave feature from higher-order modes in the range between $1.5$
and $3.5$ rad s-1. Its peak frequency is around $800$ Hz initially and
gradually increases to $900-1000$ Hz. Therefore, in addition to the
gravitational bounce signal, the detection of the $300$ Hz, $1.3$ kHz, the
higher-order mode, and even the related asymmetric emission of neutrinos,
could provide additional diagnostics to estimate the initial angular velocity
of a collapsing core.

Core-collapse supernovae (304); Gravitational wave astronomy (675);
Hydrodynamical simulations (767); Neutron stars (1108)

††software:  SPHYNX (Cabezón et al., 2017; García-Senz et al., 2022), FLASH
(Fryxell et al., 2000; Dubey et al., 2008), Matplotlib (Hunter, 2007), NumPy
(van der Walt et al., 2011), PyCWT (Torrence & Compo, 1998), SciPy (Virtanen
et al., 2019), yt (Turk et al., 2011)

##  1 INTRODUCTION

Detection of gravitational waves (GWs) together with neutrino emissions from
nearby core-collapse supernovae (CCSNe) will place meaningful constraints on
the supernova engine(s) and the nuclear equation of state. In the past decade,
tremendous efforts have been dedicated to producing GW waveform predictions
from CCSNe (Burrows & Vartanyan, 2021), including non-rotating progenitors
(O’Connor & Couch, 2018; Vartanyan & Burrows, 2020; Mezzacappa et al., 2020,
2023), light massive stars (Melson et al., 2015; Burrows et al., 2019),
rotating progenitors (Takiwaki & Kotake, 2018; Andresen et al., 2019; Powell &
Müller, 2020; Shibagaki et al., 2021; Takiwaki et al., 2021), failed
supernovae (Pan et al., 2021; Ivanov & Fernández, 2021), magnetized CCSNe
(Kuroda et al., 2020; Obergaulinger & Aloy, 2021; Jardine et al., 2022; Powell
et al., 2023), and pulsational pair-instability supernovae (Powell et al.,
2021; Rahman et al., 2022). However, these studies suggest that only galactic
CCSNe or fast-rotating progenitors have a reasonable chance of being detected
by the ongoing O4 (or the next O5) observational run of the Advanced LIGO,
Advanced Virgo, and KAGRA network (Abbott et al., 2020; Szczepańczyk et al.,
2021).

To prepare ourselves for future GW detection from a nearby CCSN, we have to
understand the GW features from CCSNe through (magneto-)hydrodynamical
simulations as those described above, and/or GW asteroseismology (Torres-Forné
et al., 2018; Morozova et al., 2018; Sotani et al., 2021). Richers et al.
(2017) conducted $1824$ axisymmetric general-relativistic (GR) hydrodynamical
simulations that include a wide range of angular velocities of a $12M_{\odot}$
progenitor. However, due to limits on computational resources, Richers et al.
(2017) calculated only up to the first $\sim 50$ ms after the core bounce,
focused on the bounce signals, and conducted the calculations exclusively with
the s12 progenitor. Pajkos et al. (2019, 2021) conducted a similar analysis
for different progenitor masses and performed longer simulations during the
accretion phase. However, in their latest work, the analysis is done with
approximately fifty 2D simulations and only four 3D simulations. Three-
dimensional fluid instabilities, such as standing acceleration shock
instability (SASI; Blondin et al. 2003) could play an important role in the GW
features (Kuroda et al., 2016), and the spiral modes of SASI can only be
developed in 3D simulations. Therefore, extending the previous works to full
3D calculations is crucial and necessary to investigate the realistic GW
features from CCSNe.

Unlike simulations of binary neutron star mergers, it is important to
highlight that gravitational waveforms produced from CCSN simulations rely on
grid-based hydrodynamics codes. Additionally, the majority of these grid-based
hydrodynamics codes adopt either Cartesian coordinates (O’Connor & Couch,
2018; Kuroda et al., 2020; Pan et al., 2021; Shibagaki et al., 2021) or
spherical coordinates (Burrows et al., 2019; Powell et al., 2021; Takiwaki et
al., 2021). Spherical grids require an inner boundary condition and/or a fixed
proto-neutron star (PNS) center. Such limitations might restrict the GW
emissions from regions close to the coordinate center or the inner boundary.
On the other hand, Cartesian grids have no inner boundary conditions at the
coordinate center but can generate $m=4$ perturbations (e.g., Ott et al.,
2007) and therefore, artificial GW emissions from these grid effects may
pollute the GW features. Furthermore, sound waves might oscillate and bounce
back between grid-refinement boundaries, which might also induce artificial GW
emissions. An alternative to avoid these artifacts is to use a pure meshless
Lagrangian method, such as smoothed particle hydrodynamics (SPH), which does
not suffer from these grid effects. Furthermore, GW emissions can be recovered
directly from the particle distribution and physical magnitudes carried by
each particle, without time derivatives involved, producing clean and reliable
waveforms (Centrella & McMillan, 1993).

An additional reason for using an SPH code in this work is its intrinsic
conservation properties and, in particular, angular momentum conservation.
Given that we are simulating from slow to fast rotators, angular momentum
conservation and transfer are key aspects to be confident in our results. SPH
codes are constructed so that their momentum and energy equations are pairwise
equivalent between particles and their neighbors (Monaghan, 2005). This,
theoretically, ensures conservation to machine precision. Note, nevertheless,
that in real simulations, this conservation is somewhat degraded by the
inclusion of self-gravity, but it is still maintained at consistently high
accuracy.

Therefore, in this work, we investigate the GW features of a wide range of
initial angular velocities of a $20M_{\odot}$ progenitor, using the SPH code
SPHYNX (Cabezón et al., 2017; García-Senz et al., 2022). A subset of these
simulations was also repeated with the grid-based adaptive mesh refinement
(AMR) code, FLASH (Fryxell et al., 2000; Dubey et al., 2008), to ensure that
our conclusions are code-independent, and when differences appear, we could
understand them within the framework of the usage of different hydrodynamics
solvers.

In the following section, we introduce the numerical methods that we employ.
Namely, the hydrodynamics codes, the physics included in them, the initial
setup including a list of all calculated models, and the extraction methods of
the GW emissions and the spherical harmonic modes. Section 3 presents our
results, focusing on the features in the GW emissions, their dependence on the
initial angular velocity, and their potential observability. Finally, Section
4 offers a summary and conclusion of the results.

##  2 NUMERICAL METHODS and MODELS

We describe the numerical codes and the corresponding setup of our simulations
in Section 2.1. In Section 2.2, we present the initial conditions of the
investigated supernova progenitor and the rotation setup. Finally, we describe
the analysis methods for obtaining the GW emissions and spherical harmonic
coefficients in Sections 2.3 and 2.4, respectively.

###  2.1 Hydrodynamics Codes

SPHYNX111https://astro.physik.unibas.ch/sphynx is a state-of-the-art SPH code,
with accurate gradient evaluation (García-Senz et al., 2012; Rosswog, 2015;
García-Senz et al., 2022), pairing-resistant interpolating kernels (Cabezón et
al., 2008), generalized volume elements (Hopkins, 2013; Saitoh & Makino,
2013), adaptive artificial viscosity via switches (Read et al., 2010), and
adaptive spatial and temporal resolution. To our knowledge, this is currently
the only SPH code that can simulate CCSNe with spectral neutrino treatment.

In the context of CCSN, SPHYNX was used for the development and comparison of
a spectral neutrino leakage scheme (Perego et al., 2014). It was also used in
a code-comparison work (Cabezón et al., 2018), where it showed particularly
good agreement with other Eulerian codes, such as FLASH, at least for the
collapse phase and the first $50$ ms postbounce, and for a series of CCSN
simulations with several physics implementations and different progenitors.

The implementation of neutrino treatment in SPHYNX is based on the Isotropic
Diffusion Source Approximation (IDSA; Liebendörfer et al., 2009) and is
coupled with a parametrized deleptonization (Liebendörfer, 2005) for the
collapse phase. Self-gravity is calculated using a Barnes-Hut algorithm on an
octree, and it also includes an effective GR potential correction (Case A in
Marek et al., 2006), adopted to replace the monopole term of the Newtonian
gravitational potential. Finally, we use the LS220 equation of state from
stellarcollapse.org (Lattimer & Swesty, 1991; O’Connor & Ott, 2010) in all of
our simulations. For further details on the coupling of all these physics
modules, we refer the reader to Section 2.4 in Cabezón et al. (2018).

For comparisons, we additionally perform two simulations with the grid-based
Eulerian hydrodynamics code FLASH222https://flash.rochester.edu that includes
the IDSA for neutrino transport (Pan et al., 2016, 2018, 2019, 2021). We also
use the same effective GR potential and LS220 equation of state.

###  2.2 Initial Setup

In this work, we employ the same $20M_{\odot}$ spherically symmetric
progenitor star with solar metallicity as the initial conditions for all
simulations. The initial radial profiles of density, temperature, electron
fraction, and radial velocity are adopted from the 1D, s20 model provided by
Woosley & Heger (2007). The initial 3D particle distribution for the SPH
calculations was achieved following the original 1D density profile but with a
random angular arrangement. Next, the resulting 3D distribution was relaxed,
allowing the particles to move only tangentially (i.e., at a fixed radius),
which smoothed out spurious clumps of particles and provided clean, smooth
initial profiles. A small external pressure ($P_{\mathrm{ext}}=1.8\times
10^{22}$ dyn cm-2) was added to all particles, similar to the pressure exerted
by the outer layers of the star, which are not included in the simulation.
This prevented particles in the very low-density region from experiencing an
artificially high gradient of pressure and escaping from the system. Note that
$P_{\mathrm{ext}}$ is only relevant to a very thin layer of particles in the
outer region of the simulated section of the star, as the pressure profile
ranges between $10^{24}-10^{33}$ dyn cm-2 for the vast majority of the domain.
Regarding neutrino transport, we use $20$ energy bins, logarithmically spaced
from $3$ until $300$ MeV for the electron-flavor neutrinos. The heavier $\mu$
and $\tau$ neutrinos and anti-neutrinos are considered as a single species
that cools the system via a leakage scheme.

Finally, an artificial tangential velocity was imposed to induce rotation. We
assigned the initial rotational profile as follows,

$\Omega(r)=\dfrac{\Omega_{0}}{1+\left(\frac{r}{A}\right)^{2}},$ (1)

where $\Omega$ is the angular velocity as a function of the spherical radius
$r$, $A=1.03376\times 10^{8}$ cm is a characteristic length that is taken from
the fitting formula provided in Pajkos et al. (2019) (see Figure 2 therein),
and $\Omega_{0}$ is the angular velocity at the center of the star, of which
we investigate several values ranging from $0.0$ to $4.0$ rad s-1.

We perform the SPH simulations with three different resolutions, which contain
$200$k, $675$k, and $1600$k particles that cover the central $10^{4}$ km of
the progenitor star. The corresponding resolutions within $10$ km in the
radius are about $600$, $400$, and $300$ m after the core bounce. The primary
analysis and results presented in this paper are based on the high-resolution
simulations with $1600$k particles. Furthermore, we used the simulations with
$675$k particles to investigate a larger parameter space of $\Omega_{0}$.

Table 1: Summary of the models. Columns from left to right show the model name, numerical code, initial angular velocity ($\Omega_{0}$), number of particles, peak frequency of the kHz signal, the value of $T/|W|$ at the initial time, and the highest spatial resolution at $20$ ms postbounce ($\Delta x_{\mathrm{min,20ms}}$). Model | Code | $\Omega_{0}$ | # of particles | $f_{\mathrm{peak}}$ |  $T/|W|$init | $\Delta x_{\mathrm{min,20ms}}$
---|---|---|---|---|---|---
|  | [rad s-1] | [$10^{3}$] | [kHz] | $[10^{-3}]$ | [m]
S00 | SPHYNX | $0.0$ | $675$ | - | $4.84\times 10^{-10}$ | $372.3$
S05 | SPHYNX | $0.5$ | $675$ | - | $0.12$ | $372.4$
S10 | SPHYNX | $1.0$ | $675$ | - | $0.48$ | $374.4$
S15 | SPHYNX | $1.5$ | $675$ | - | $1.09$ | $378.5$
S20 | SPHYNX | $2.0$ | $675$ | $1.35$ | $1.94$ | $383.7$
S25 | SPHYNX | $2.5$ | $675$ | $1.29$ | $3.02$ | $389.3$
S30 | SPHYNX | $3.0$ | $675$ | $1.29$ | $4.35$ | $395.4$
S35 | SPHYNX | $3.5$ | $675$ | $1.34$ | $5.93$ | $403.0$
S40 | SPHYNX | $4.0$ | $675$ | - | $7.74$ | $410.9$
S30L | SPHYNX | $3.0$ | $200$ | $1.28$ | $4.35$ | $591.9$
S20H | SPHYNX | $2.0$ | $1600$ | $1.38$ | $1.94$ | $288.1$
S30H | SPHYNX | $3.0$ | $1600$ | $1.12-1.32$ | $4.36$ | $296.9$
F20 | FLASH | $2.0$ | - | $1.01-1.48$ | $1.89$ | $488.3$
F30 | FLASH | $3.0$ | - | $1.21-1.39$ | $4.25$ | $488.3$

The Cartesian grid setup in our FLASH simulations closely follows the setup
described in Pan et al. (2021). Therefore, we just give here a brief review of
our setup. We use a three-dimensional simulation box that covers the inner
$10^{4}$ km of the CCSN progenitor, and employ nine levels of AMR, which yield
a maximum spatial resolution of $488$ m at the highest AMR level. The central
$r<120$ km sphere has the highest spatial resolution, while we reduce the AMR
level as we move farther away from the center of the progenitor to save
computing time. This results in an effective angular resolution $\sim
0^{\circ}.2-0^{\circ}.4$. For the outer boundary conditions, we use a power-
law profile that depends on the spherical radius. In addition, we adopt the
same $20$ neutrino energy bins as in SPHYNX, spaced logarithmically from $3$
to $300$ MeV for the electron flavor neutrinos and a leakage scheme for the
$\mu$ and $\tau$ neutrinos. We use the same s20 progenitor from Woosley &
Heger (2007) and take the same rotational profile as in Equation (1) with
$\Omega_{0}=2$ and $3$ rad s-1.

Table 1 summarizes the relevant information for all our models. Letters F and
S in the model name denote the code used to perform the simulation, FLASH or
SPHYNX, respectively. The number in the model names shows the adopted value of
the initial central angular velocity, $\Omega_{0}$. Finally, L and H stand for
low and high resolution, respectively, for SPHYNX models. The column
$f_{\mathrm{peak}}$ shows the resulting peak GW frequency around kHz in the
time interval $t_{\mathrm{pb}}=10-100$ ms. We provide a range of peak
frequencies for models that exhibit noticeable variations in their peak
frequencies over time. The sixth column of Table 1 indicates the ratio between
rotational energy and gravitational energy, $T/|W|$, calculated for the
regions with $\rho\geq 10^{6}$ g cm-3 at the initial time. The last column
$\Delta x_{\mathrm{min,20ms}}$ shows the highest resolution for each model,
which is defined as the smallest smoothing length for SPHYNX models or the
smallest cell size for FLASH models, at $20$ ms postbounce.

###  2.3 Extracting the Gravitational Waves

In order to extract the GW emissions in SPHYNX, we adopt the transverse-
traceless gauge to cover the far zone of the source. Hence, we have only two
polarizations of the amplitude of the GWs:

$\displaystyle h_{+}$
$\displaystyle=\frac{1}{D}\frac{G}{c^{4}}\left(\hbox{\"{I}}\kern-7.5pt\kern-5.0pt\raise
2.45004pt\hbox{ \vrule
width=7.5pt,height=0.5pt}_{\theta\theta}-\hbox{\"{I}}\kern-7.5pt\kern-5.0pt\raise
2.45004pt\hbox{ \vrule width=7.5pt,height=0.5pt}_{\phi\phi}\right),$ (2)
$\displaystyle h_{\times}$
$\displaystyle=\frac{1}{D}\frac{G}{c^{4}}\hbox{\"{I}}\kern-7.5pt\kern-5.0pt\raise
2.45004pt\hbox{ \vrule width=7.5pt,height=0.5pt}_{\theta\phi}\,,$ (3)

where $G$ is gravitational constant, $c$ is the speed of light, and $D$ is the
distance to the source. The reduced quadrupole
($\hbox{I}\kern-3.61111pt\kern-5.0pt\raise 3.07497pt\hbox{ \vrule
width=3.61111pt,height=0.5pt}_{lm}$) is defined in Cartesian coordinates as,

$\hbox{I}\kern-3.61111pt\kern-5.0pt\raise 3.07497pt\hbox{ \vrule
width=3.61111pt,height=0.5pt}_{lm}=\int\rho\left(x_{l}x_{m}-\frac{1}{3}\delta_{lm}x_{k}x^{k}\right)d^{3}x\,.$
(4)

where $l$ and $m$ are components of the position vector $\mathbf{x}$,
$\delta_{lm}$ is a $\delta$-Kronecker, and $\rho$ is the local density.

The components of the reduced quadrupole in spherical coordinates are related
to the Cartesian by (Oohara et al., 1997),

$\displaystyle\hbox{\"{I}}\kern-7.5pt\kern-5.0pt\raise 2.45004pt\hbox{ \vrule
width=7.5pt,height=0.5pt}_{\theta\theta}$
$\displaystyle=\left(\hbox{\"{I}}\kern-7.5pt\kern-5.0pt\raise 2.45004pt\hbox{
\vrule
width=7.5pt,height=0.5pt}_{xx}\cos^{2}\phi+\hbox{\"{I}}\kern-7.5pt\kern-5.0pt\raise
2.45004pt\hbox{ \vrule
width=7.5pt,height=0.5pt}_{yy}\sin^{2}\phi+\hbox{\"{I}}\kern-7.5pt\kern-5.0pt\raise
2.45004pt\hbox{ \vrule width=7.5pt,height=0.5pt}_{xy}\sin
2\phi\right)\cos^{2}\theta$
$\displaystyle\qquad+\hbox{\"{I}}\kern-7.5pt\kern-5.0pt\raise 2.45004pt\hbox{
\vrule
width=7.5pt,height=0.5pt}_{zz}\sin^{2}\theta-\left(\hbox{\"{I}}\kern-7.5pt\kern-5.0pt\raise
2.45004pt\hbox{ \vrule
width=7.5pt,height=0.5pt}_{xz}\cos\phi+\hbox{\"{I}}\kern-7.5pt\kern-5.0pt\raise
2.45004pt\hbox{ \vrule width=7.5pt,height=0.5pt}_{yz}\sin\phi\right)\sin
2\theta\,,$ (5) $\displaystyle\hbox{\"{I}}\kern-7.5pt\kern-5.0pt\raise
2.45004pt\hbox{ \vrule width=7.5pt,height=0.5pt}_{\phi\phi}$
$\displaystyle=\hbox{\"{I}}\kern-7.5pt\kern-5.0pt\raise 2.45004pt\hbox{ \vrule
width=7.5pt,height=0.5pt}_{xx}\sin^{2}\phi+\hbox{\"{I}}\kern-7.5pt\kern-5.0pt\raise
2.45004pt\hbox{ \vrule
width=7.5pt,height=0.5pt}_{yy}\cos^{2}\phi-\hbox{\"{I}}\kern-7.5pt\kern-5.0pt\raise
2.45004pt\hbox{ \vrule width=7.5pt,height=0.5pt}_{xy}\sin 2\phi\,,$ (6)
$\displaystyle\hbox{\"{I}}\kern-7.5pt\kern-5.0pt\raise 2.45004pt\hbox{ \vrule
width=7.5pt,height=0.5pt}_{\theta\phi}$
$\displaystyle=-\frac{1}{2}\left(\hbox{\"{I}}\kern-7.5pt\kern-5.0pt\raise
2.45004pt\hbox{ \vrule
width=7.5pt,height=0.5pt}_{xx}-\hbox{\"{I}}\kern-7.5pt\kern-5.0pt\raise
2.45004pt\hbox{ \vrule width=7.5pt,height=0.5pt}_{yy}\right)\cos\theta\sin
2\phi+\hbox{\"{I}}\kern-7.5pt\kern-5.0pt\raise 2.45004pt\hbox{ \vrule
width=7.5pt,height=0.5pt}_{xy}\cos\theta\cos 2\phi$
$\displaystyle\qquad-\left(\hbox{\"{I}}\kern-7.5pt\kern-5.0pt\raise
2.45004pt\hbox{ \vrule
width=7.5pt,height=0.5pt}_{xz}\sin\phi-\hbox{\"{I}}\kern-7.5pt\kern-5.0pt\raise
2.45004pt\hbox{ \vrule
width=7.5pt,height=0.5pt}_{yz}\cos\phi\right)\sin\theta\,.$ (7)

Therefore, to obtain the GW waveforms we compute the six non-zero components
of $\hbox{\"{I}}\kern-7.5pt\kern-5.0pt\raise 2.45004pt\hbox{ \vrule
width=7.5pt,height=0.5pt}_{lm}$ in Cartesian coordinates and then transform
them into spherical coordinates using Equations (5), (6), and (7). Finally, we
substitute these components into Equations (2) and (3) to obtain the final
waveforms.

There are different methods for evaluating the components of
$\hbox{\"{I}}\kern-7.5pt\kern-5.0pt\raise 2.45004pt\hbox{ \vrule
width=7.5pt,height=0.5pt}_{lm}$. Nevertheless, time derivatives cause
numerical difficulties due to two main reasons: the numerical noise introduced
by the discretization and the magnification of the high-frequency components
of the noise. To avoid these problems, we opted to use the method proposed by
Centrella & McMillan (1993), which takes advantage of the Lagrangian nature of
SPH and is similar to the stress formula of Finn & Evans (1990).

A discretized version of Equation (4) is

$\hbox{I}\kern-3.61111pt\kern-5.0pt\raise 3.07497pt\hbox{ \vrule
width=3.61111pt,height=0.5pt}_{lm}=\sum_{i}m^{i}\left[x_{l}^{i}x_{m}^{i}-\frac{1}{3}\delta_{lm}\mathbf{r}^{i}\cdot\mathbf{r}^{i}\right]\,,$
(8)

where the subscripts refer to the three Cartesian components, the superscripts
label the SPH particles, $m^{i}$ is the mass of each particle, and the
summation is over all the particles of the simulation. Taking the second time
derivative of Equation (8) is straightforward, obtaining

$\displaystyle\hbox{\"{I}}\kern-7.5pt\kern-5.0pt\raise 2.45004pt\hbox{ \vrule
width=7.5pt,height=0.5pt}_{lm}$
$\displaystyle=\frac{2}{3}\sum_{i}m^{i}\left[2v_{l}^{i}v_{m}^{i}+a_{l}^{i}x_{m}^{i}+x_{l}^{i}a_{m}^{i}\right.$
$\displaystyle\qquad+\left.\delta_{lm}\left(v_{l}^{i}v_{m}^{i}+x_{l}^{i}a_{m}^{i}-\mathbf{v}^{i}\cdot\mathbf{v}^{i}-\mathbf{r}^{i}\cdot\mathbf{a}^{i}\right)\right]\,.$
(9)

The $x$, $v$, and $a$ terms are the components of the position, velocity, and
acceleration vectors of particle $i$, respectively. Using Equation (9) in
Equations (2)$-$(7) we can find the amplitude of both polarizations of the GW
emissions directly from magnitudes calculated by the SPH code without having
to compute explicit second-time derivatives.

FLASH also uses the quadrupole formula to extract the GW emissions (Equations
2$-$7), but unlike SPHYNX, the GW calculations in FLASH follow the strain
formulation to compute the first time derivative of the quadrupole moment. The
second time derivative of the quadrupole moment is evaluated by a finite
difference method via post-processing (Pan et al., 2018, 2021).

###  2.4 Spherical Harmonic Mode Analysis

It is known that, for a moderately rotating core, the low-$T/|W|$ instability
could develop after the core bounce, which induces an $m=1$ or $m=2$
deformation that results in quasi-sinusoidal oscillations of the GW emissions
(Ott et al., 2005; Shibagaki et al., 2020). To investigate the relationship
between the GW signals and the deformation induced by the low-$T/|W|$
instability, we apply a spherical harmonic decomposition to the resulting
fluid distributions. Following Burrows et al. (2012), we evaluate the
coefficient of each mode using

$a_{lm}=\frac{(-1)^{|m|}}{\sqrt{4\pi(2l+1)}}\frac{\sum
wR(\theta,\phi)Y_{l}^{m}(\theta,\phi)}{\sum w},$ (10)

where $R$ is the distance from the PNS center, and $w$ is a weighting function
taken as the volume. For SPHYNX, we set the weighting function to the
associated volume of the SPH particles, $w=m^{i}/\rho^{i}$, where $m^{i}$ and
$\rho^{i}$ are the mass and density carried by the SPH particles,
respectively. The orthonormal harmonic basis functions, $Y_{l}^{m}$, are
expressed as

$Y_{l}^{m}(\theta,\phi)=\begin{cases}\sqrt{2}N_{l}^{m}P_{l}^{m}(\cos\theta)\cos
m\phi&m>0,\\\ N_{l}^{0}P_{l}^{0}(\cos\theta)&m=0,\\\
\sqrt{2}N_{l}^{|m|}P_{l}^{|m|}(\cos\theta)\sin|m|\phi&m<0,\end{cases}$ (11)

where

$N_{l}^{m}=\sqrt{\frac{2l+1}{4\pi}\frac{(l-m)!}{(l+m)!}}$ (12)

and $P_{l}^{m}(\cos\theta)$ is the associated Legendre polynomial.

In this work, we focus on the dipole mode ($l=1$, $m=1$) and quadrupole mode
($l=2,m=2$), which represent the $m=1$ and $m=2$ deformations, respectively.

##  3 RESULTS

###  3.1 Dynamics Overview

Figure 1: Panels from left to right and from top to bottom describe the time
evolution of the central density ($\rho_{c}$), the averaged shock radius
($R_{\mathrm{sh}}$), the mass accretion rate measured at $r=200$ km and $500$
km ($\dot{M}_{200}$ and $\dot{M}_{500}$), the enclosed mass within a density
contour of $10^{11}$ g cm-3 and $10^{13}$ g cm-3 ($M_{11}$ and $M_{13}$), and
the averaged isodensity radii corresponding to $M_{11}$ and $M_{13}$ ($R_{11}$
and $R_{13}$) in the models S20H, S30H, F20, and F30. Figure 2: Entropy
distribution on the XZ-plane at $t_{\mathrm{pb}}=50$ ms in models S30H (top)
and F30 (bottom). The arrows represent the velocity field. Figure 3:
Spherically-averaged, radial profiles of various fluid and neutrino quantities
at different postbounce times for models S30H (blue) and F30 (orange). The
profiles at different times are shifted cumulatively by the offset labeled in
each panel, where the black dotted lines denote the corresponding zero point.

We first summarize and compare the evolution of the collapsar in our
simulations performed with SPHYNX and FLASH. A broad overview of the time
evolution of several PNS quantities is shown in Figure 1 for models with
$\Omega_{0}=2$ and $3$ rad s-1 (models S20H, S30H, F20, and F30). The panels
in the first row of Figure 1 show the time evolution of the central density
and averaged shock radius, respectively, which are roughly in agreement
between both hydrodynamics codes. As pointed out in Cabezón et al. (2018),
SPHYNX simulations generally produce slightly higher central densities than
FLASH simulations. This is due to the different definitions of the central
density. In FLASH, the central density is taken from the cell-averaged density
of the single densest cell ($\Delta x=488$ m), while for SPHYNX, it is
evaluated as an average of the $50$ densest SPH particles, with an equivalent
cell width of less than $300$ m. From the averaged shock radius evolution, we
find that models with both initial angular velocities exhibit a very rapid
explosion, where the averaged shock radii exceed $500$ km within $70-100$ ms
postbounce. This is much faster than the previous work with slower initial
angular velocities (Takiwaki et al., 2016; Andresen et al., 2019; Pan et al.,
2021; Shibagaki et al., 2021), but comparable with the explosion times in
Kuroda et al. (2014).

In the second row of Figure 1, we show the mass accretion rates measured at
$r=200$ and $500$ km, which are in good agreement between codes in the models
with $\Omega_{0}=2$ rad s-1 (solid lines). The sudden drops in the mass
accretion rates occur when the shock front reaches the measured radii at
$t_{\mathrm{pb}}\sim 40$ and $80$ ms, respectively, transitioning from infall
to outflow. The $\Omega_{0}=3$ cases (dashed lines) show behavior similar to
that of the $\Omega_{0}=2$ cases, though the deviations become larger when the
shock front reaches the measured radii. The more substantial negative mass
accretion rate in the SPHYNX model (S30H) implies a stronger explosion
compared to the FLASH model (F30). This interpretation can also be
comprehended from the entropy distribution as well. Figure 2 shows slice color
plots of the entropy distribution at $50$ ms postbounce along the XZ-plane,
perpendicular to the equatorial plane with the $z$-axis representing the
rotation axis, using models S30H and F30. Both models have similar shock
expansion and convective regions, but S30H has a more spherically symmetric
shock expansion at around $300$ km due to poor shock resolution ($\Delta x\sim
25$ km) at this low-density region. In consequence, the negative accretion
rate contributed from the outflow near the equatorial plane is negated by the
inflow through the pole in model F30, leading to a less negative mass
accretion rate compared to model S30H. Furthermore, the drops in the mass
accretion in model S30H occur later than model F30 due to the slower averaged
shock expansion.

The panels in the third row of Figure 1 show the evolution of the enclosed
mass within a density contour of $10^{11}$ g cm-3 ($M_{11}$) and $10^{13}$ g
cm-3 ($M_{13}$), respectively, and the corresponding isodensity radii
($R_{11}$ and $R_{13}$) in the last row. The radius at $R_{11}$ usually
coincides with the neutrino sphere, and therefore, the enclosed mass $M_{11}$
is used to represent the PNS mass. The enclosed mass $M_{13}$ roughly
describes the PNS inner core. In the models with $\Omega_{0}=2$ rad s-1, model
S20H has a slightly lower enclosed mass $M_{11}$ compared to model F20. This
is mainly due to the lower mass accretion rates in SPHYNX around core bounce,
while the PNS radius $R_{11}$ remains similar. In the cases of $\Omega_{0}=3$
rad s-1, the deviation in the PNS mass $M_{11}$ increases further after
$t_{\mathrm{pb}}=40$ ms when the low-$T/|W|$ instability starts to develop and
generates stronger asymmetric accretion at later stages postbounce. On the
other hand, the PNS inner core mass $M_{13}$ shows the opposite behavior. In
both angular velocities, SPHYNX models have more compact inner cores than
their FLASH counterparts due to better core resolutions in models S20H and
S30H (see Table 1).

Typically, when there is no rotation, the central density and the PNS inner
core mass $M_{13}$ are expected to increase over time due to ongoing mass
accretion and PNS cooling. However, in the case of rotating progenitors, we
observe that model S20H has a slight increase in the inner core mass $M_{13}$
within $100$ ms postbounce, whereas model F20 shows almost no increase within
the first $40$ ms postbounce and is followed by a slight decrease in the inner
core mass. These distinctions may be due to differences in angular momentum
transport and conservation between SPHYNX and FLASH as described in Cabezón et
al. (2018). In the models with $\Omega_{0}=3$ rad s-1, when angular momentum
keeps propagating inward, the centrifugal force eventually overcomes the
gravitational force and, therefore, it induces a decrease in the PNS inner
core mass $M_{13}$. Both models S30H and F30 behave similarly but with
different decreasing rates.

In addition, we also find that FLASH simulations have a notable neutron star
kick and a modulation motion ($v_{\mathrm{pns}}\sim 240$ and $360$ km s-1 in
models F20 and F30, respectively) due to an asymmetric explosion and together
with numerical artifacts on angular momentum non-conservation. These motions
affect the evolution of the PNS in FLASH simulations, especially at late time.
On the other hand, the neutron star kick velocities in S20H and S30H are less
than $1$ km s-1. In the following sections, we show that these differences in
the compactness and motion of the PNS inner core will affect the GW signatures
of the low-$T/|W|$ instability.

A comparison of radial profiles of various fluid and neutrino quantities
between the models S30H and F30 at different postbounce times is shown in
Figure 3. The spherically-averaged profiles of density, electron fraction, and
entropy are consistent in the radius coordinate. The reason for the lower
angular momentum observed in the inner core of the PNS in model F30 is
attributed to the fact that angular momentum conservation has deteriorated,
which is primarily caused by inherent numerical dissipation in grid-based
hydrodynamics codes (as discussed in Cabezón et al. 2018). When considering
the quantities of electron-type neutrinos, it should be noted that model S30H
has a more compact and hotter inner core in its PNS in comparison to the
corresponding PNS in model F30 (see Figure 1). This leads to a higher
production rate and energy level of electron neutrinos within the PNS inner
core in model S30H. Overall, the radial profiles of neutrino quantities at
different times are consistent between the models S30H and F30.

Figure 4: GW strains and spectrograms of the plus mode for the models with
$\Omega_{0}=2$ and $3$ rad s-1, using SPHYNX with $1600$k particles (top
panels) and FLASH (bottom panels), seen along the equatorial plane at a
distance of $10$ kpc. The white line represents the doubled dynamical
frequency at a density of $10^{13}$ g cm-3 (see Section 3.3 for a more
detailed description).

###  3.2 GW Features and Origin in Rapidly Rotating CCSNe

In this section, we discuss the GW features in the models with $\Omega_{0}=2$
and $3$ rad s-1, using SPHYNX with $1600$k particles (models S20H and S30H)
and FLASH (models F20 and F30). The top panels in Figure 4 show the plus mode
of GW emissions from models S20H and S30H, seen along the equatorial plane at
a distance of $10$ kpc. The GW strain is shown at the top, and the
corresponding spectrogram is displayed at the bottom in each panel. The
spectrogram is computed using the wavelet analysis implemented in the
PyCWT333https://github.com/regeirk/pycwt code, a Python package based on
Torrence & Compo (1998), where the power spectrum is divided by the wavelet
scales to rectify the energy bias (Liu et al., 2007). In addition, we divide
the resulting amplitude by the square root of the sampling rate to ensure
consistent strength between different sampling rates.

Figure 5: GW spectrograms of the plus mode emitted from regions with density
$\rho<10^{11}$ g cm-3 (left), $10^{11}\leq\rho<10^{13}$ g cm-3 (middle), and
$\rho\geq 10^{13}$ g cm-3 (right) for models S30H (top) and F30 (bottom), seen
along the pole at a source distance of $10$ kpc. In the middle and right
panels, the white and red contours show the spectrogram of the normalized
spherical harmonic coefficients $a_{11}$ and $a_{22}$, respectively, where the
frequency of $a_{11}$ is doubled to ease comparison.

First, both SPHYNX models show distinctive bounce and ring-down signals in the
first $20$ ms postbounce, which have been extensively studied in previous
works (e.g., Abdikamalov et al., 2014; Richers et al., 2017; Abdikamalov et
al., 2022, and references therein). After the bounce and ring-down signals, we
can observe that the GW emissions exhibit quasi-sinusoidal time oscillations
starting at $t_{\mathrm{pb}}=20-30$ ms in both models. We find that this GW
feature is simultaneous with the so-called low-$T/|W|$ instability (Saijo et
al., 2003; Ott et al., 2005, 2007; Scheidegger et al., 2008, 2010; Kuroda et
al., 2014; Shibagaki et al., 2020). The low-$T/|W|$ instability is a non-
axisymmetric rotational instability that develops in cores with a high degree
of differential rotation around the corotation radius, where the pattern
frequency of the induced oscillation is equal to the local angular frequency
of the background flow (Centrella et al., 2001; Watts et al., 2005; Saijo &
Yoshida, 2006). This low-$T/|W|$ instability can induce $m=1$ and/or $m=2$
deformations that lead to a time-changing quadrupole moment, which is the
ultimate source of the GW emissions.

From the spectrograms of S20H and S30H in Figure 4, we can identify two strong
GW signals in the frequency ranges from $200$ to $400$ Hz and from $1100$ to
$1400$ Hz. Hereafter we refer to them as the $300$ Hz signal and the kHz
signal, respectively. In the bottom panels of Figure 4, we present the GW
strains obtained using the FLASH code with $\Omega_{0}=2$ rad s-1 (F20) and
$\Omega_{0}=3$ rad s-1 (F30), and their corresponding spectrograms. The GW
emissions in both models exhibit similar bounce, ring-down, $300$ Hz, and kHz
signals as discussed above for the SPHYNX simulations. Although the bounce and
ring-down signals in models F20 and F30 are consistent with the S20H and S30H
models, we can see that the $300$ Hz and kHz signals evolve differently,
especially with respect to the occurrence time and peak frequency of the kHz
signal. As discussed in Section 3.1, FLASH and SPHYNX show slightly different
dynamical evolutions of the PNS after core bounce (see Figure 1). In Section
3.3, we will show that the $300$ Hz and kHz signals are correlated with the
outer and inner structures of the PNS, respectively, and therefore causes the
differences between FLASH and SPHYNX models.

To investigate the origin of the $300$ Hz and kHz signals, we calculate the
contributions of GW emissions from different density regions by post-
processing the simulation data, using the formulae described in Section 2.3.
Figure 5 shows the GW contributions from regions within density
$\rho<10^{11}$, $10^{11}\leq\rho<10^{13}$, and $\rho\geq 10^{13}$ g cm-3 for
models S30H and F30, seen along the pole to eliminate bounce and ring-down
signals. In both models, we can see that the $300$ Hz and kHz signals emanate
from separate regions. The $300$ Hz signal is mainly from the region of
$10^{11}\leq\rho<10^{13}$ g cm-3, while the kHz signal is mainly from the PNS
inner core where $\rho\geq 10^{13}$ g cm-3.

Figure 6: Normalized mode amplitude at a radius of $15$ km for models S20H and
S30H (top), and F20 and F30 (bottom).

In Figure 5, we also evaluate the spherical harmonic components $a_{11}$ and
$a_{22}$ in the regions of $10^{11}\leq\rho<10^{13}$ and $\rho\geq 10^{13}$ g
cm-3 (defined in Section 2.4), and overlap their spectrograms in white and red
lines, respectively. The components $a_{11}$ and $a_{22}$ represent the $m=1$
and $m=2$ deformations in the specific density region, and their contours
indicate the mode frequency of the corresponding deformation,
$f_{\mathrm{mode},m}$. The mode frequency is related to the pattern frequency
by $f_{\mathrm{pat},m}=f_{\mathrm{mode},m}/m$ (Watts et al., 2005). Note that
the mode frequency of $a_{11}$ is doubled in Figure 5 to facilitate comparison
between the components $a_{11}$ and $a_{22}$. Comparing the spherical harmonic
modes with the $300$ Hz and kHz GW signals reveals that the $a_{22}$ component
coincides with both. This is because the dominant quadrupole component of GW
emissions stems from the $l=2,m=2$ mode, making it a natural source for both
signals. On the other hand, the pattern frequencies of $a_{11}$ and $a_{22}$
satisfy the relation $f_{\mathrm{pat},1}\simeq f_{\mathrm{pat},2}$, indicating
that the $a_{22}$ component is a daughter mode of the $a_{11}$ component. This
infers that both the $300$ Hz and kHz signals are associated with the $m=1$
spiral deformation induced by the low-$T/|W|$ instability. Shibagaki et al.
(2020) conducted a full-GR CCSN simulation of a rapid-rotating $70M_{\odot}$
progenitor. They found a transient quasi-periodic time modulation at $450$ Hz
from the $m=1$ spiral deformation in $50-100$ km. The $300$ Hz signal in our
models is similar to the $450$ Hz signal in Shibagaki et al. (2020). On the
other hand, the kHz signal resembles the $\sim 930$ Hz signal found by Ott et
al. (2007) in CCSN simulations of a $20M_{\odot}$ progenitor, which is
correlated with the $m=1$ mode at $10-15$ km, but the peak frequency is higher
in our cases. In addition to the $300$ Hz and kHz signals, some higher-order
modes of GW emissions at around $800$ Hz can be seen in Figure 5 as well. In
model S30H, the $800$ Hz signal is correlated to the kHz signal and the
$a_{22}$ component but is much weaker than the $300$ Hz and kHz signals. In
model F30, similar higher-order modes also exist between the $300$ Hz and kHz
signals, but the interactions among these GW features are more complex than
that in model S30H.

Figure 7: Radial profiles of rotational (blue), Brunt-Väisälä (orange), and
Lamb (green) frequencies at $t_{\mathrm{pb}}=30$ ms for models F30 (left) and
S30H (right). The black solid lines denote the half-peak frequencies of the
$300$ Hz and kHz gravitational-wave signals. Figure 8: Spectrograms of the
normalized spherical harmonic coefficients, $a_{11}$ (white, doubled
frequency) and $a_{22}$ (red), of the mean energy of trapped electron anti-
neutrinos for models S30H (left) and F30 (right). The GW spectrograms of the
plus mode, seen along the pole at $10$ kpc, are shown in filled contour to
ease comparison.

To investigate the potential impact of $m=4$ perturbations induced by the
Cartesian grid discretization in FLASH simulations, we also perform an
analysis of the azimuthal density modes in the equatorial plane by computing
the Fourier amplitude (Centrella et al., 2001)

$C_{m}(\varpi)=\frac{1}{2\pi}\int^{2\pi}_{0}\rho(\varpi,z=0)e^{im\phi}d\phi,$
(13)

where $\varpi$ is the cylindrical radius relative to the PNS center. Figure 6
shows the normalized mode amplitude, $|C_{m}|/C_{0}$, evaluated at a radius of
$\varpi=15$ km. We note that the relative difference in the mean density,
$C_{0}$, is below $30\%$ before $t_{\mathrm{pb}}\sim 35$ ms between SPHYNX and
FLASH simulations and then increases to $80\%$ due to the PNS kicks in FLASH
simulations, making the normalized mode amplitudes in FLASH simulations higher
than the SPHYNX counterpart simulations. Furthermore, we can see that in FLASH
simulations, the $m=4$ mode is the dominant mode in the early postbounce phase
($t_{\mathrm{pb}}\lesssim 10$ ms) due to grid effects, while it is always
subdominant in SPHYNX simulations, as expected for a meshless method. However,
there is no clear relation between the $m=4$ grid mode and the other
$m=\\{1,2,3\\}$ modes in FLASH simulations, and thus both F20 and F30 models
remain dynamically stable to grid perturbations. This is also consistent with
the earlier work done with the full GR simulations via the Whisky code in Ott
et al. (2007).

Recently, Takiwaki et al. (2021) proposed that the low-$T/|W|$ instability in
CCSN environments could be triggered by the Rossby waves growing near the
convective zone. In their explanation, such instability requires having a
corotation radius to coincide with the convective layer in the PNS. Therefore,
it is interesting to investigate whether the $300$ Hz and kHz low-$T/|W|$
signals in our models are satisfied with the same criteria. Figure 7 shows the
radial profiles of the rotational, Brunt-Väisälä, and Lamb frequencies,
evaluated in the equatorial plane, for models F30 and S30H at
$t_{\mathrm{pb}}=30$ ms. We evaluate the Lamb frequency via

$f_{\mathrm{Lamb}}=\frac{1}{2\pi}\frac{\sqrt{l(l+1)}c_{s}}{r},$ (14)

and the Brunt-Väisälä frequency using the Ledoux criterion (Buras et al.,
2006; Ott et al., 2013)

$f_{\mathrm{BV}}=\frac{\mathrm{sign}(C_{L})}{2\pi}\sqrt{\left|\frac{C_{L}}{\rho}\frac{d\Phi}{dr}\right|},$
(15)

where $l$ is taken to be $1$, $c_{s}$ is the local speed of sound, $d\Phi$ is
the local gravitational potential where the approximation $d\Phi/dr\sim-
GM(r)/r^{2}$ is adopted. The Ledoux criterion reads (Ledoux, 1947)

$C_{L}=-\left(\frac{\partial\rho}{\partial
P}\right)_{s,Y_{l}}\left[\left(\frac{\partial P}{\partial
s}\right)_{\rho,Y_{l}}\left(\frac{ds}{dr}\right)+\left(\frac{\partial
P}{\partial Y_{l}}\right)_{\rho,s}\left(\frac{dY_{l}}{dr}\right)\right],$ (16)

where we approximate the lepton fraction by the electron fraction, $Y_{l}\sim
Y_{e}$, for simplicity. In model F30, two corotation radii at $5-10$ km and
$30$ km, which are described by the intersection of the pattern frequencies of
the $300$ Hz and kHz signals (black lines) and the rotational frequency (blue
line), are in convective regions with negative Brunt-Väisälä frequency. This
is consistent with the statement proposed by Takiwaki et al. (2021). On the
other hand, in the case of model S30H, only one corotation radius of the kHz
signal at $5-10$ km coincides with the edge of a convective layer. However, we
note that the Brunt-Väisälä frequency oscillates rapidly between $30$ and $40$
km in S30H, suggesting that a convective zone could be developing in that
region as well.

In addition to the Rossby wave scenario, it is worth mentioning that the area
where the kHz signal originated, approximately at around the isodensity radius
$R_{13}$ ($r\sim 20$ km), aligns with the region where electron anti-neutrinos
are predominantly generated and remain coupled with the matter. Therefore, the
density variation driven by the $m=1$ deformation can affect the production of
electron anti-neutrinos in the PNS inner core. In the IDSA neutrino treatment,
electron-type neutrinos are decomposed into trapped and free-streaming
neutrinos. Among those, only trapped neutrinos are coupled with the fluid.
Since in our implementation of IDSA, the $\mathcal{O}(v/c)$ terms of the
neutrino pressure are included (see Equation 24 in Liebendörfer et al. 2009),
the neutrino pressure gradient could contribute to the development the $m=1$
deformation. To establish this, we apply the spherical harmonic decomposition
to the mean energy of trapped electron anti-neutrinos through a variant of
Equation (10):

$a_{lm}=\frac{(-1)^{|m|}}{\sqrt{4\pi(2l+1)}}\frac{\sum
w\left<E_{\bar{\nu}_{e}}\right>(\theta,\phi)Y_{l}^{m}(\theta,\phi)}{\sum w},$
(17)

where $\left<E_{\bar{\nu}_{e}}\right>=Z_{\bar{\nu}_{e}}/Y_{\bar{\nu}_{e}}$ is
the mean energy of trapped electron anti-neutrinos. We use the mean energy
$\left<E_{\bar{\nu}_{e}}\right>$, instead of $Y_{\bar{\nu}_{e}}$ or
$Z_{\bar{\nu}_{e}}$, to avoid the sharp decline in electron anti-neutrino
fractions, which could introduce strong numerical noise. We evaluate the
spherical harmonic components in the shell $R_{13}\pm 10$ km. Figure 8 shows
the spectrograms of the $a_{11}$ and $a_{22}$ components for models S30H and
F30, revealing that the distributions of trapped electron anti-neutrinos are
in good agreement with not only the kHz signal but also the $300$ Hz signal.
Albeit the asymmetric electron anti-neutrino distribution requires the $m=1$
deformation from the low-$T/|W|$ instability to grow. In the following
section, we show that there are distinct phase differences between the $m=1$
deformation and the asymmetric neutrino distribution in the PNS inner core.
This suggests that neutrino pressure could play a role in promoting the
development of $m=1$ deformation and exhibiting unique GW signatures in the
kHz window. We also find that the distribution of electron neutrinos has a
similar asymmetric effect but is less pronounced than electron anti-neutrinos.
This is because electron neutrinos have been produced since the collapse.

###  3.3 GW Dependence on Initial Angular Velocity

Figure 9: Time evolution of the PNS mass with density $\rho\geq 10^{13}$ g
cm-3 (top), the corresponding averaged isodensity radius (middle), and the
doubled dynamical frequency (bottom) in the SPHYNX simulations with $675$k
particles (solid lines). Different colors represent simulations with different
initial angular velocities. Dashed lines show the counterpart models (S20H and
S30H) with $1600$k particles for comparison. Figure 10: GW strains and
spectrograms of the plus mode seen along the pole at a distance of $10$ kpc.
Different panels represent the SPHYNX $675$k particle simulations with
different initial angular velocity. The white lines represent the doubled
dynamical frequency for the domain with $\rho\geq 10^{13}$ g cm-3 in each
model.

To investigate the dependence of the $300$ Hz and kHz low-$T/|W|$ signals on
the initial angular velocity, we perform a series of SPHYNX simulations with
$675$k particles (See Table 1). Figure 9 shows the time evolution of the
enclosed mass ($M_{13}$), the corresponding isodensity radius ($R_{13}$), and
the doubled dynamical frequency, $2f_{\rm dyn}\sim
2\sqrt{R_{13}^{3}/GM_{13}}/2\pi$, for the domain with $\rho\geq 10^{13}$ g
cm-3 for models with $\Omega_{0}$ ranging from $0.0$ to $4.0$ rad s-1, in
steps of $0.5$ rad s-1. In addition, we also plot the same quantities for
models S20H and S30H in dashed lines for comparison. First, we can see that
the PNS inner core structures are consistent with different particle numbers,
and the trend basically follows the description provided in Section 3.1 for
models with $\Omega_{0}=2$ and $3$ rad s-1: higher initial angular velocities
will result in a decrease of the enclosed mass and isodensity radius, while
the doubled dynamical frequency remains a relatively unchanged range of kHz
along the simulations.

Figure 10 shows the corresponding GW strain and spectrogram for these models,
assuming viewing from the pole and at a distance of $10$ kpc. We note that in
Figure 10, the color ranges are fixed at the same values for all panels to
facilitate comparison between different initial angular velocities. We first
focus on the models S20 and S30 to ensure that our simulations with $675$k
particles accurately capture the main GW features observed in our simulations
with $1600$k particles (models S20H and S30H). By comparing spectrograms in
Figure 4 and Figure 10, we find that both the $300$ Hz and kHz signals in
models S20 and S30 are similar to those of models S20H and S30H. Since the
main low-$T/|W|$ features discussed in our high-resolution runs are also
captured in the simulations with $675$k particles, we can conclude that we are
able to conduct a parameter study of the initial angular velocity at a
considerably lower computational cost using this resolution.

Figure 11: Normalized density fluctuation in the equatorial plane at
$t_{\mathrm{pb}}=50$ ms for models using SPHYNX with $675$k particles and
different initial angular velocities. The average density, $\bar{\rho}$, is
taken over the azimuthal direction in the plane. The black dashed curves
denote the isodensity contours at $\rho=10^{11}$, $10^{12}$, $10^{13}$, and
$10^{14}$ g cm-3. Figure 12: Mean energy of trapped electron anti-neutrino
distribution on the XZ-plane at $t_{\mathrm{pb}}=50$ ms for models using
SPHYNX with $675$k particles and different initial angular velocities. The red
and blue curves respectively illustrate positive and negative density
fluctuations centered around zero. The white dashed curves denote the
isodensity contours at $\rho=10^{12}$, $10^{13}$, and $10^{14}$ g cm-3.

As a result, from Figure 10 we find that the $300$ Hz signal starts to appear
when the initial angular velocity $\Omega_{0}\geq 1$ rad s-1 and it persists
consistently across different angular velocities. However, the kHz signal
appears only in a range within $1.5\leq\Omega_{0}\leq 3.5$ rad s-1, though the
kHz signal in model S35 occurred in a short duration between $30-60$ ms
postbounce. In addition, we overlay the doubled dynamical frequency on the
spectrograms in Figure 10. We find that the kHz signal, once it appears, has a
peak frequency described by the doubled dynamical frequency at the density of
$10^{13}$ g cm-3. As we have discussed in Section 3.2, the $m=1$ deformation
associated with the kHz signal originates from the PNS inner core where
$\rho>10^{13}$ g cm-3, and its pattern frequency ($a_{11}$ in Figure 5) should
be related to the characteristic frequency, the dynamical frequency here, in
this density region. Therefore, it is not surprising that the kHz signal can
be described by the doubled dynamical frequency at the density of $10^{13}$ g
cm-3.

To visualize the $m=1$ deformation, we follow Takiwaki et al. (2016) to show
the density fluctuation in the equatorial plane,
$(\rho-\bar{\rho})/\bar{\rho}$, at $t_{\mathrm{pb}}=50$ ms in Figure 11, where
$\bar{\rho}$ is the azimuthally averaged density. Different isodensity
contours are plotted as dashed lines as well. We can see that the $m=1$
deformation can develop in several density regions. In the density region of
$10^{11}\lesssim\rho\lesssim 10^{13}$ g cm-3, all models with $\Omega_{0}\geq
1.0$ rad s-1 show the spiral structure in the density fluctuation plot, but
not in models with $\Omega_{0}\leq 0.5$ rad s-1. This is consistent with the
presence of the $300$ Hz signal in Figure 10 and also confirms that the $300$
Hz signal is indeed emitted mainly from this density region in Figure 5. On
the other hand, we can see that in the region of $\rho>10^{13}$ g cm-3, where
the kHz signal emanated, the $m=1$ deformation develops only in the models
with $1.5\leq\Omega_{0}\leq 3.5$ rad s-1. This supports our hypothesis that
the kHz signal is associated with the low-$T/|W|$ instability. As a
consequence, the strength of the kHz signal is correlated with the presence
and amplitude of the $m=1$ deformation in the high-density region of
$\rho>10^{13}$ g cm-3.

Figure 12 shows the mean energy of trapped electron anti-neutrinos for models
with different initial angular velocities. We can see that the mean energy of
electron anti-neutrinos also exhibits prominent asymmetric distributions at
$R_{13}\sim 20$ km in the models with $1.5\leq\Omega_{0}\leq 3.5$ rad s-1,
which is consistent with the models for which the $m=1$ density deformation
and the kHz signal are present. This alignment once again highlights the
correlation between the $m=1$ density deformation and the neutrino
distribution. For neutrino treatments that include the $\mathcal{O}(v/c)$
terms, the density asymmetries can lead to an asymmetric neutrino pressure,
potentially contributing to the development of the $m=1$ density deformation
in the PNS inner core. Comparing Figure 11 and Figure 12, we find that the
phase of $m=1$ density deformation leads the phase of the mean energy
variation around $R_{13}$ (the system rotates counter-clockwise). This
suggests that the high-density region produces more energetic neutrinos,
subsequently heating the surrounding matter and facilitating expansion, which
eventually results in a low-density region later.

A similar effect but with a lower modulation frequency has been proposed by
Takiwaki & Kotake (2018). They suggest that the neutrino emissions have a time
modulation similar to the GW frequency from the low-$T/|W|$ instability ($\sim
100-300$ Hz) and could be detectable by the Hyper-Kamiokande and the IceCube
detectors. It is also expected that the asymmetric neutrino distributions
associated with the kHz GW signal should produce noticeable time modulations
on the neutrino emissions. However, in our current implementation of the IDSA,
the free-streaming neutrinos are averaged over angles and therefore cannot be
directly evaluated without additional approximations.

###  3.4 Detectability of GW Signals

In this section, we discuss the detectability of the GW features from rapid-
rotating CCSNe using the current ground-based GW detectors. Figure 13 shows
the amplitude spectral density (ASD) of the plus mode of GW emissions between
$t_{\mathrm{pb}}=-10$ and $100$ ms in the SPHYNX simulations with $675$k
particles, seen along the pole and assumed at a distance $10$ kpc. Note that
the bounce and ring-down signals are not present in this viewing angle which
makes it easier to focus on the GW signals related to the low-$T/|W|$
instability. The sensitivity curves of the GW detectors Advanced LIGO,
Advanced Virgo, and KAGRA (Abbott et al., 2020) are plotted as black lines for
reference. It is clear from Figure 13 that both the $300$ Hz and kHz signals
are detectable by the current ground-based GW detectors at our assumed source
distance, and they are among the strongest signals in this time period. The
peak frequencies of the $300$ Hz and kHz signals fall within narrow ranges
between $210-300$ Hz and $1280-1350$ Hz, respectively, and are not sensitive
to the initial angular velocity whenever these signals are present.

Figure 13: Amplitude spectral density of the plus mode of GW emissions between
$t_{\mathrm{pb}}=-10$ and $100$ ms, for models using SPHYNX with $675$k
particles and various initial angular velocities, seen along the pole at a
distance of $10$ kpc. The black solid, dashed, and dotted lines denote the
sensitive curves of advanced LIGO (aLIGO), advanced Virgo (AdV), and KAGRA,
respectively. Figure 14: Similar to Figure 13, but for the $\Omega_{0}=3.0$
models using different codes and resolutions. For the comparison, only the GW
emissions between $t_{\mathrm{pb}}=-10$ and $100$ ms are used for analysis.

In addition, it is worth noting that there is one weaker but noticeable GW
signal associated with the higher-order mode at around $800$ Hz, as discussed
in Section 3.2. By examining Figures 10 and 13, we can see that the peak
frequency of this higher-order mode signal increases gradually after
$t_{\mathrm{pb}}=70-100$ ms, which leads to a secondary peak around $900-1000$
Hz in the ASD. Among these GW features, the $300$ Hz signal is the most robust
signal from the low-$T/|W|$ instability and can be excited when the initial
angular velocity is higher than $1.0$ rad s-1. The kHz signal will appear when
the initial angular velocity is within the range $1.5\leq\Omega_{0}\leq 3.5$
rad s-1. The cases with $\Omega_{0}=2.5$ and $3.0$ rad s-1 have the strongest
emissions, suggesting a resonance frequency at around $2.5-3.0$ rad s-1. The
higher-order mode signal, which peaked around $800$ Hz, appears in a similar
range as the kHz signal ($1.5\leq\Omega_{0}\leq 3.5$ rad s-1) but does not
show a clear resonance frequency. Therefore, if we could detect the higher-
order mode signal in addition to the $300$ Hz and kHz signals, these would
provide additional constraints on narrowing down the initial angular velocity
of a collapsar.

In Figure 14, we compare the ASD of GW emissions of SPHYNX simulations with
different numbers of particles, using $\Omega_{0}=3.0$ rad s-1. We also plot
the ASD of model F30 for comparison. The ASD is evaluated within the time
range between $t_{\mathrm{pb}}=-10$ and $100$ ms. We can see that the GW
emissions have a qualitatively similar spectral density distribution between
the SPHYNX simulations with $200$k (S30L), $675$k (S30), and $1600$k particles
(S30H). However, model S30L shows a kHz signal that is one order of magnitude
weaker in amplitude compared to models S30 and S30H, and the peak frequency is
also $30-50$ Hz lower. We consider that these differences are due to an
underresolved PNS inner core in the lowest resolution model (S30L). Even with
its low resolution, the S30L model displays an ASD similar to those of the S30
and S30H models, which have converged results. This is because the primary
contributions to the GW emissions originate from the densest regions of the
PNS, which have the highest spatial resolution in SPHYNX. In addition, as long
as the overall evolution is accurately described, the GW emissions can be
adequately followed with integration over the entire domain, as discussed in
Section 2.3.

Comparing codes, the peak frequencies of the $300$ Hz and kHz low-$T/|W|$
signals in the FLASH simulation (F30) are qualitatively similar to those of
the models S30 and S30H. Nevertheless, the GW emission in the F30 model is
approximately one order of magnitude stronger in the signal amplitude than the
results from SPHYNX. As pointed out in Andresen et al. (2021), the grid
resolution in the post-shock region will affect the numerical damping of the
convective cells and the activities of the g-mode around the PNS, which could
explain the variation in the magnitudes of the GW emissions between SPHYNX and
FLASH.

##  4 SUMMARY & CONCLUSIONS

We have performed an analysis of the development of the low-$T/|W|$
instability and associated GW emissions in the early postbounce phase of
rotating CCSNe. To this end, we performed 3D hydrodynamical core-collapse
simulations of a $20M_{\odot}$ progenitor with different initial angular
velocities ($\Omega_{0}$), which are parametrically added to the progenitor
model. In this work, we computed models with $0\leq\Omega_{0}\leq 4$ rad s-1
using the smoothed particle hydrodynamics code SPHYNX, and models with
$\Omega_{0}=2$ and $3$ rad s-1 using the grid-based hydrodynamics code FLASH
for comparison.

Among our rotating models, the GW emissions exhibit two strong low-$T/|W|$
signals after $20$ ms postbounce in both SPHYNX and FLASH simulations. Both
signals are correlated with the $m=1$ deformation induced by the low-$T/|W|$
instability, with peak frequencies of about $300$ Hz and $1.3$ kHz (called kHz
here). The $300$ Hz signal is present in models with $\Omega_{0}\geq 1.0$ rad
s-1, which mainly emanated from the region of $10^{11}\leq\rho<10^{13}$ g
cm-3. The peak frequency of the $300$ Hz signal is not sensitive to the
initial angular velocity, the code used, or the spatial resolution. On the
other hand, the kHz signal is present only in models with a narrower range of
initial angular velocity, $1.5\leq\Omega_{0}\leq 3.5$ rad s-1, originates
mainly in the region of $\rho\geq 10^{13}$ g cm-3, and is highly associated
with the asymmetric distribution of electron anti-neutrinos. The peak
frequency of the kHz signal, once present, is not sensitive to the initial
angular velocity, but is moderately affected by the dynamical evolution of the
PNS inner core.

In addition to the $300$ Hz and kHz signals, there is an additional weaker
higher-order mode of GW emissions at around $800$ Hz, emanating mainly from
the region of $10^{11}\leq\rho<10^{13}$ g cm-3, in those of our models that
developed the low-$T/|W|$ instability. This higher-order signal is also
correlated with the $m=1$ deformation in the PNS, but its occurrence time is
different from those of the $300$ Hz and kHz signals. The range of initial
angular velocity where the higher-order signal exists is similar to that of
the kHz signal, around $1.5\leq\Omega_{0}\leq 3.5$ rad s-1. However, the peak
frequency of this higher-order signal gradually increases to $900-1000$ Hz
after $70-100$ ms postbounce, which leads to a secondary peak in the amplitude
spectral density. Therefore, the initial angular velocity of the CCSN
progenitors can be inferred from the detection of the higher-order signal, and
the $300$ Hz and kHz low-$T/|W|$ signals.

We note that the GW features and their peak frequencies presented in this work
can be dependent on the numerical methods and the physical models used,
especially the kHz signal. Furthermore, the kHz features associated with the
low-$T/|W|$ instability are typically originated from the asymmetric density
distribution in the PNS inner core, where electron anti-neutrinos are largely
produced and still coupled with the matter. The density asymmetries in this
region can induce an asymmetric neutrino distribution and consequently result
in an asymmetric neutrino pressure. This, in turn, could facilitate the
development of density deformation. A more in-depth stability analysis is
necessary and will be one of our future works.

In this work, we have not considered the effects of the magnetic field, which
can affect the dynamical evolution of the PNS and the explosion mechanism
(e.g., Kuroda et al., 2020; Matsumoto et al., 2020; Müller & Varma, 2020;
Obergaulinger & Aloy, 2020, 2021; Raynaud et al., 2022; Powell et al., 2023).
In the presence of a strong magnetic field, the low-$T/|W|$ instability could
be suppressed (Fu & Lai, 2011; Muhlberger et al., 2014). The angular momentum
transport driven by the magnetorotational instability can redistribute the
rotational profile (Bugli et al., 2018), which in turn affects the development
of the low-$T/|W|$ instability and weakens the associated GW signals by an
order of magnitude (Bugli et al., 2023). To obtain a more complete diagnostic
of the angular momentum profile from the low-$T/|W|$ signals, 3D
magnetohydrodynamics simulations should be performed.

We are grateful to the referee for a thoughtful report and useful suggestions
that helped us improve the manuscript. This work is supported by the National
Center for Theoretical Sciences of Taiwan, the Ministry of Education (Higher
Education Sprout Project NTU-112L104022), the National Science and Technology
Council of Taiwan through grant NSTC 111-2112-M-007-037 and
112-2811-M-002-113, the Center for Informatics and Computation in Astronomy
(CICA) at National Tsing Hua University through a grant from the Ministry of
Education of Taiwan, and the Swiss Platform for Advanced Scientific Computing
(PASC) projects SPH-EXA and SPH-EXA2: Optimizing Smooth Particle Hydrodynamics
for Exascale Computing. KCP is supported by the NSTC grant NSTC
111-2112-M-007-037 and 112-2112-M-007-040. This work has also been carried out
as part of the SKACH consortium through funding from SERI. Simulations and
data analysis have been carried out on the Taiwania supercomputer at the
National Center for High-Performance Computing (NCHC) in Taiwan, on the CICA
Cluster at the National Tsing Hua University, and on the scientific computing
core facility sciCORE (http://scicore.unibas.ch/) at the University of Basel.
Analysis and visualization of the simulation data were completed using the
analysis toolkit yt (Turk et al., 2011).

## References

  * Abbott et al. (2020) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2020, Living Reviews in Relativity, 23, 3, doi: 10.1007/s41114-020-00026-9
  * Abdikamalov et al. (2014) Abdikamalov, E., Gossan, S., DeMaio, A. M., & Ott, C. D. 2014, Phys. Rev. D, 90, 044001, doi: 10.1103/PhysRevD.90.044001
  * Abdikamalov et al. (2022) Abdikamalov, E., Pagliaroli, G., & Radice, D. 2022, in Handbook of Gravitational Wave Astronomy. Edited by C. Bambi, 21, doi: 10.1007/978-981-15-4702-7_21-1
  * Andresen et al. (2021) Andresen, H., Glas, R., & Janka, H. T. 2021, MNRAS, 503, 3552, doi: 10.1093/mnras/stab675
  * Andresen et al. (2019) Andresen, H., Müller, E., Janka, H. T., et al. 2019, MNRAS, 486, 2238, doi: 10.1093/mnras/stz990
  * Blondin et al. (2003) Blondin, J. M., Mezzacappa, A., & DeMarino, C. 2003, ApJ, 584, 971, doi: 10.1086/345812
  * Bugli et al. (2023) Bugli, M., Guilet, J., Foglizzo, T., & Obergaulinger, M. 2023, MNRAS, 520, 5622, doi: 10.1093/mnras/stad496
  * Bugli et al. (2018) Bugli, M., Guilet, J., Müller, E., et al. 2018, MNRAS, 475, 108, doi: 10.1093/mnras/stx3158
  * Buras et al. (2006) Buras, R., Rampp, M., Janka, H. T., & Kifonidis, K. 2006, A&A, 447, 1049, doi: 10.1051/0004-6361:20053783
  * Burrows et al. (2012) Burrows, A., Dolence, J. C., & Murphy, J. W. 2012, ApJ, 759, 5, doi: 10.1088/0004-637X/759/1/5
  * Burrows et al. (2019) Burrows, A., Radice, D., & Vartanyan, D. 2019, MNRAS, 485, 3153, doi: 10.1093/mnras/stz543
  * Burrows & Vartanyan (2021) Burrows, A., & Vartanyan, D. 2021, Nature, 589, 29, doi: 10.1038/s41586-020-03059-w
  * Cabezón et al. (2017) Cabezón, R. M., García-Senz, D., & Figueira, J. 2017, A&A, 606, A78, doi: 10.1051/0004-6361/201630208
  * Cabezón et al. (2008) Cabezón, R. M., García-Senz, D., & Relaño, A. 2008, Journal of Computational Physics, 227, 8523, doi: 10.1016/j.jcp.2008.06.014
  * Cabezón et al. (2018) Cabezón, R. M., Pan, K.-C., Liebendörfer, M., et al. 2018, A&A, 619, A118, doi: 10.1051/0004-6361/201833705
  * Centrella & McMillan (1993) Centrella, J. M., & McMillan, S. L. W. 1993, ApJ, 416, 719, doi: 10.1086/173272
  * Centrella et al. (2001) Centrella, J. M., New, K. C. B., Lowe, L. L., & Brown, J. D. 2001, ApJ, 550, L193, doi: 10.1086/319634
  * Dubey et al. (2008) Dubey, A., Reid, L. B., & Fisher, R. 2008, Physica Scripta Volume T, 132, 014046, doi: 10.1088/0031-8949/2008/T132/014046
  * Finn & Evans (1990) Finn, L. S., & Evans, C. R. 1990, ApJ, 351, 588, doi: 10.1086/168497
  * Fryxell et al. (2000) Fryxell, B., Olson, K., Ricker, P., et al. 2000, ApJS, 131, 273, doi: 10.1086/317361
  * Fu & Lai (2011) Fu, W., & Lai, D. 2011, MNRAS, 413, 2207, doi: 10.1111/j.1365-2966.2011.18296.x
  * García-Senz et al. (2012) García-Senz, D., Cabezón, R. M., & Escartín, J. A. 2012, A&A, 538, A9, doi: 10.1051/0004-6361/201117939
  * García-Senz et al. (2022) —. 2022, A&A, 659, A175, doi: 10.1051/0004-6361/202141877
  * Hopkins (2013) Hopkins, P. F. 2013, MNRAS, 428, 2840, doi: 10.1093/mnras/sts210
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Ivanov & Fernández (2021) Ivanov, M., & Fernández, R. 2021, ApJ, 911, 6, doi: 10.3847/1538-4357/abe59e
  * Jardine et al. (2022) Jardine, R., Powell, J., & Müller, B. 2022, MNRAS, 510, 5535, doi: 10.1093/mnras/stab3763
  * Kuroda et al. (2020) Kuroda, T., Arcones, A., Takiwaki, T., & Kotake, K. 2020, ApJ, 896, 102, doi: 10.3847/1538-4357/ab9308
  * Kuroda et al. (2016) Kuroda, T., Kotake, K., & Takiwaki, T. 2016, ApJ, 829, L14, doi: 10.3847/2041-8205/829/1/L14
  * Kuroda et al. (2014) Kuroda, T., Takiwaki, T., & Kotake, K. 2014, Phys. Rev. D, 89, 044011, doi: 10.1103/PhysRevD.89.044011
  * Lattimer & Swesty (1991) Lattimer, J. M., & Swesty, D. F. 1991, Nucl. Phys. A, 535, 331, doi: 10.1016/0375-9474(91)90452-C
  * Ledoux (1947) Ledoux, P. 1947, ApJ, 105, 305, doi: 10.1086/144905
  * Liebendörfer (2005) Liebendörfer, M. 2005, ApJ, 633, 1042, doi: 10.1086/466517
  * Liebendörfer et al. (2009) Liebendörfer, M., Whitehouse, S. C., & Fischer, T. 2009, ApJ, 698, 1174, doi: 10.1088/0004-637X/698/2/1174
  * Liu et al. (2007) Liu, Y., San Liang, X., & Weisberg, R. H. 2007, Journal of Atmospheric and Oceanic Technology, 24, 2093, doi: 10.1175/2007JTECHO511.1
  * Marek et al. (2006) Marek, A., Dimmelmeier, H., Janka, H. T., Müller, E., & Buras, R. 2006, A&A, 445, 273, doi: 10.1051/0004-6361:20052840
  * Matsumoto et al. (2020) Matsumoto, J., Takiwaki, T., Kotake, K., Asahina, Y., & Takahashi, H. R. 2020, MNRAS, 499, 4174, doi: 10.1093/mnras/staa3095
  * Melson et al. (2015) Melson, T., Janka, H.-T., & Marek, A. 2015, ApJ, 801, L24, doi: 10.1088/2041-8205/801/2/L24
  * Mezzacappa et al. (2020) Mezzacappa, A., Marronetti, P., Landfield, R. E., et al. 2020, Phys. Rev. D, 102, 023027, doi: 10.1103/PhysRevD.102.023027
  * Mezzacappa et al. (2023) —. 2023, Phys. Rev. D, 107, 043008, doi: 10.1103/PhysRevD.107.043008
  * Monaghan (2005) Monaghan, J. J. 2005, Reports on Progress in Physics, 68, 1703, doi: 10.1088/0034-4885/68/8/R01
  * Morozova et al. (2018) Morozova, V., Radice, D., Burrows, A., & Vartanyan, D. 2018, ApJ, 861, 10, doi: 10.3847/1538-4357/aac5f1
  * Muhlberger et al. (2014) Muhlberger, C. D., Nouri, F. H., Duez, M. D., et al. 2014, Phys. Rev. D, 90, 104014, doi: 10.1103/PhysRevD.90.104014
  * Müller & Varma (2020) Müller, B., & Varma, V. 2020, MNRAS, 498, L109, doi: 10.1093/mnrasl/slaa137
  * Obergaulinger & Aloy (2020) Obergaulinger, M., & Aloy, M. Á. 2020, MNRAS, 492, 4613, doi: 10.1093/mnras/staa096
  * Obergaulinger & Aloy (2021) —. 2021, MNRAS, 503, 4942, doi: 10.1093/mnras/stab295
  * O’Connor & Ott (2010) O’Connor, E., & Ott, C. D. 2010, Classical and Quantum Gravity, 27, 114103, doi: 10.1088/0264-9381/27/11/114103
  * O’Connor & Couch (2018) O’Connor, E. P., & Couch, S. M. 2018, ApJ, 865, 81, doi: 10.3847/1538-4357/aadcf7
  * Oohara et al. (1997) Oohara, K.-i., Nakamura, T., & Shibata, M. 1997, Progress of Theoretical Physics Supplement, 128, 183, doi: 10.1143/PTPS.128.183
  * Ott et al. (2007) Ott, C. D., Dimmelmeier, H., Marek, A., et al. 2007, Phys. Rev. Lett., 98, 261101, doi: 10.1103/PhysRevLett.98.261101
  * Ott et al. (2005) Ott, C. D., Ou, S., Tohline, J. E., & Burrows, A. 2005, ApJ, 625, L119, doi: 10.1086/431305
  * Ott et al. (2013) Ott, C. D., Abdikamalov, E., Mösta, P., et al. 2013, ApJ, 768, 115, doi: 10.1088/0004-637X/768/2/115
  * Pajkos et al. (2019) Pajkos, M. A., Couch, S. M., Pan, K.-C., & O’Connor, E. P. 2019, ApJ, 878, 13, doi: 10.3847/1538-4357/ab1de2
  * Pajkos et al. (2021) Pajkos, M. A., Warren, M. L., Couch, S. M., O’Connor, E. P., & Pan, K.-C. 2021, ApJ, 914, 80, doi: 10.3847/1538-4357/abfb65
  * Pan et al. (2018) Pan, K.-C., Liebendörfer, M., Couch, S. M., & Thielemann, F.-K. 2018, ApJ, 857, 13, doi: 10.3847/1538-4357/aab71d
  * Pan et al. (2021) —. 2021, ApJ, 914, 140, doi: 10.3847/1538-4357/abfb05
  * Pan et al. (2016) Pan, K.-C., Liebendörfer, M., Hempel, M., & Thielemann, F.-K. 2016, ApJ, 817, 72, doi: 10.3847/0004-637X/817/1/72
  * Pan et al. (2019) Pan, K.-C., Mattes, C., O’Connor, E. P., et al. 2019, Journal of Physics G Nuclear Physics, 46, 014001, doi: 10.1088/1361-6471/aaed51
  * Perego et al. (2014) Perego, A., Gafton, E., Cabezón, R., Rosswog, S., & Liebendörfer, M. 2014, A&A, 568, A11, doi: 10.1051/0004-6361/201423755
  * Powell & Müller (2020) Powell, J., & Müller, B. 2020, MNRAS, 494, 4665, doi: 10.1093/mnras/staa1048
  * Powell et al. (2023) Powell, J., Müller, B., Aguilera-Dena, D. R., & Langer, N. 2023, MNRAS, 522, 6070, doi: 10.1093/mnras/stad1292
  * Powell et al. (2021) Powell, J., Müller, B., & Heger, A. 2021, MNRAS, 503, 2108, doi: 10.1093/mnras/stab614
  * Rahman et al. (2022) Rahman, N., Janka, H. T., Stockinger, G., & Woosley, S. E. 2022, MNRAS, 512, 4503, doi: 10.1093/mnras/stac758
  * Raynaud et al. (2022) Raynaud, R., Cerdá-Durán, P., & Guilet, J. 2022, MNRAS, 509, 3410, doi: 10.1093/mnras/stab3109
  * Read et al. (2010) Read, J. I., Hayfield, T., & Agertz, O. 2010, MNRAS, 405, 1513, doi: 10.1111/j.1365-2966.2010.16577.x
  * Richers et al. (2017) Richers, S., Ott, C. D., Abdikamalov, E., O’Connor, E., & Sullivan, C. 2017, Phys. Rev. D, 95, 063019, doi: 10.1103/PhysRevD.95.063019
  * Rosswog (2015) Rosswog, S. 2015, MNRAS, 448, 3628, doi: 10.1093/mnras/stv225
  * Saijo et al. (2003) Saijo, M., Baumgarte, T. W., & Shapiro, S. L. 2003, ApJ, 595, 352, doi: 10.1086/377334
  * Saijo & Yoshida (2006) Saijo, M., & Yoshida, S. 2006, MNRAS, 368, 1429, doi: 10.1111/j.1365-2966.2006.10229.x
  * Saitoh & Makino (2013) Saitoh, T. R., & Makino, J. 2013, ApJ, 768, 44, doi: 10.1088/0004-637X/768/1/44
  * Scheidegger et al. (2008) Scheidegger, S., Fischer, T., Whitehouse, S. C., & Liebendörfer, M. 2008, A&A, 490, 231, doi: 10.1051/0004-6361:20078577
  * Scheidegger et al. (2010) Scheidegger, S., Käppeli, R., Whitehouse, S. C., Fischer, T., & Liebendörfer, M. 2010, A&A, 514, A51, doi: 10.1051/0004-6361/200913220
  * Shibagaki et al. (2020) Shibagaki, S., Kuroda, T., Kotake, K., & Takiwaki, T. 2020, MNRAS, 493, L138, doi: 10.1093/mnrasl/slaa021
  * Shibagaki et al. (2021) —. 2021, MNRAS, 502, 3066, doi: 10.1093/mnras/stab228
  * Sotani et al. (2021) Sotani, H., Takiwaki, T., & Togashi, H. 2021, Phys. Rev. D, 104, 123009, doi: 10.1103/PhysRevD.104.123009
  * Szczepańczyk et al. (2021) Szczepańczyk, M. J., Antelis, J. M., Benjamin, M., et al. 2021, Phys. Rev. D, 104, 102002, doi: 10.1103/PhysRevD.104.102002
  * Takiwaki & Kotake (2018) Takiwaki, T., & Kotake, K. 2018, MNRAS, 475, L91, doi: 10.1093/mnrasl/sly008
  * Takiwaki et al. (2021) Takiwaki, T., Kotake, K., & Foglizzo, T. 2021, MNRAS, 508, 966, doi: 10.1093/mnras/stab2607
  * Takiwaki et al. (2016) Takiwaki, T., Kotake, K., & Suwa, Y. 2016, MNRAS, 461, L112, doi: 10.1093/mnrasl/slw105
  * Torrence & Compo (1998) Torrence, C., & Compo, G. P. 1998, Bulletin of the American Meteorological Society, 79, 61, doi: 10.1175/1520-0477(1998)079<0061:APGTWA>2.0.CO;2
  * Torres-Forné et al. (2018) Torres-Forné, A., Cerdá-Durán, P., Passamonti, A., & Font, J. A. 2018, MNRAS, 474, 5272, doi: 10.1093/mnras/stx3067
  * Turk et al. (2011) Turk, M. J., Smith, B. D., Oishi, J. S., et al. 2011, ApJS, 192, 9, doi: 10.1088/0067-0049/192/1/9
  * van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science and Engineering, 13, 22, doi: 10.1109/MCSE.2011.37
  * Vartanyan & Burrows (2020) Vartanyan, D., & Burrows, A. 2020, ApJ, 901, 108, doi: 10.3847/1538-4357/abafac
  * Virtanen et al. (2019) Virtanen, P., Gommers, R., Burovski, E., et al. 2019, scipy/scipy: SciPy 1.3.2, v1.3.2, Zenodo, Zenodo, doi: 10.5281/zenodo.3533894
  * Watts et al. (2005) Watts, A. L., Andersson, N., & Jones, D. I. 2005, ApJ, 618, L37, doi: 10.1086/427653
  * Woosley & Heger (2007) Woosley, S. E., & Heger, A. 2007, Phys. Rep., 442, 269, doi: 10.1016/j.physrep.2007.02.009